function M = coaxial_circle(Rp,Rs,rw,distance)
a=Rp;
b=Rs;
d=distance;
mu_not=4*pi*10^-7;
% 
if d==0 && (Rp-Rs)==0
  k=sqrt(4*a*(a-rw)/(2*a-rw)^2);
  [K,E]=ellipke(k^2);
    M=mu_not*sqrt(a*(a-rw))*((2/k-k)*K-2/k*E);
else
k=2*sqrt(a*b/((a+b)^2+d^2));
[K,E]=ellipke(k^2);
M=2*mu_not/k*sqrt(a*b)*((1-k^2/2)*K-E);
end
end