/*
####### Maxima package file (hthickcoil.mac) #######
This code is based on equation 10 from "Calculation of the Magnetic Field Created by a Thick Coil". DOI:10.1163/156939310791958653.
Note: Equ. 10 has two errors. These are fixed in the code.
Errata:
1) In Hz(r,z) plus sign is replaced by minus sign.
2) In EtaZ term z[j] is replaced with z[k].
---------------------------------------------------------
Save this code in hthickcoil.mac file,then copy the file in one of the folders returned by 'file_search_maxima' command. Then execute load("hthickcoil.mac")$
For automatic load of Hring function one can put this line into maxima-init.mac file: setup_autoload("hthickcoil.mac",Hring)$
Usage: Hring(r,z,r1,r2,z1,z2,current).
r and z are cylindrical coordinates of point where you seek the H field.
r1 - inner radius, r2 - outer radius, z1 - lower position, z2 - upper position.
For using built-in elliptic integrals in MAXIMA, see comments FIX1 and FIX2.
Maxima 5.31.2 version has issue with elliptic_ec function of negative argument.
*/
ellipticK(m):=block(
[qtemp1:0.0,qtemp2:0.0],qtemp1:float(qtemp1),
[qtemp1,qtemp2,qtemp2,qtemp2]:quad_qags(1/sqrt(1-m*sin(x)^2),x,0,%pi/2),
return (qtemp1)
)$
ellipticE(m):=block(
[qtemp1:0.0,qtemp2:0.0],qtemp1:float(qtemp1),
[qtemp1,qtemp2,qtemp2,qtemp2]:quad_qags(sqrt(1-m*sin(x)^2),x,0,%pi/2),
return (qtemp1)
)$
Hring(r,z,r_1,r_2,z_1,z_2,I):=
block
(
[Hr,Hz,H,j,theta,v,a,b,fi,kappa,d,c,f,eps:10^(-15),Eta_r,Eta_z,bracket1:0,bracket2:0,qtemp1,qtemp2],
r:float(r),z:float(z),r_1:float(r_1),r_2:float(r_2),z_1:float(z_1),z_2:float(z_2),I:float(I),j:float(j),
H:float(H),Hr:float(Hr),Hz:float(Hz),v:float(v),a:float(a),b:float(b),fi:float(fi),kappa:float(kappa),d:float(d),
c:float(c),f:float(f),bracket1:float(bracket1),bracket2:float(bracket2), qtemp1:float(qtemp1),
/* FIX1: uncomment line below if elliptic_kc/ec(nagative_argument) works. */
/* local (K,E), K(m) := elliptic_kc(m),E(m):=elliptic_ec(m), */
/* FIX2: if elliptic_kc/ec(nagative_argument) does not work use line below.Otherwise if you use FIX1 above then comment line below. */
local (K,E), K(m) := ellipticK(m),E(m):=ellipticE(m),
array(r,flonum,2),array(z,flonum,2),
r[1]:r_1,r[2]:r_2,z[1]:z_1,z[2]:z_2,
j:I/((r_2-r_1)*(z_2-z_1)),
v:r*sin(theta),
for i:1 thru 2 do (
for k:1 thru 2 do (
a: r^2+r[i]^2+(z-z[k])^2,
b:2*r*r[i],
fi:(4/3)*sqrt(a-b)/(b+eps),
kappa: 2*b/(b-a),
d:r[i]-r*cos(theta),
c:sqrt(a-b*cos(theta)),
f:z-z[k]+eps,
Eta_r:r*cos(theta)^2*log(d+c),
Eta_z:r*cos(theta)*log(f+c)+v*atan((d*f*csc(theta)/(r*c+eps)))-z[k]+v*atan(v/f),
[qtemp1,qtemp2,qtemp2,qtemp2]:( quad_qags(Eta_r,theta,0+eps,%pi-eps) + quad_qags(Eta_r,theta,%pi+eps,2*%pi-eps) ),
bracket1: bracket1 + ((-1)^(i+k))*(qtemp1+fi*((a+b)*K(kappa)-a*E(kappa))),
[qtemp1,qtemp2,qtemp2,qtemp2]:( quad_qags(Eta_z-f*log(d+c),theta,0+eps,%pi-eps)+quad_qags(Eta_z-f*log(d+c),theta,%pi+eps,2*%pi-eps) ),
bracket2: bracket2 + ((-1)^(i+k))*qtemp1
) /* end_for_k */
), /* end_for_i */
Hr:ev((j/(4*%pi))*bracket1,numer), Hz:ev((j/(4*%pi))*bracket2,numer),
/* OUTPUT: [HR, Hz, H] */
return ([Hr,Hz,sqrt(Hr^2+Hz^2)])
)$
