Demagnetization coefficient of an ellipsoid

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)) * (ellipticF(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);