楕円体の反磁界係数

if(a==b){a=a*1.001;}
if(a==c){a=a*1.001;}
if(b==c){a=a*1.0001;b=b*1.0001;}
k = sqrt((1-(b/a)^2)/(1-(c/a)^2));
θ = acos(c/a);
x = sin(θ);
L = (b/a)*(c/a)/((sin(θ)^3)*(k^2))*(ellipticF(x,k)-ellipticE(x,k));
M = (b/a)*(c/a)/((sin(θ)^3)*(k^2)*(1-k^2))*(ellipticE(x,k)-(1-k^2)*ellipticF(x,k)-(k^2)*(sin(θ))*(c/a)/(b/a));
N = (b/a)*(c/a)/((sin(θ)^3)*(1-k^2))*((sin(θ)*(b/a)/(c/a))-ellipticE(x,k));
μeff = 1/(L+(1-L)/μr);
μeffmax = 1/L;
print(a);
print(b);
print(c);
print(k);
print(θ);
print(L);
print(M);
print(N);
print(μr);
print(μeff);
print(μeffmax);