Maxima code based on equation 10 from "Calculation of the Magnetic Field Created by a Thick Coil"

/* 

####### 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.

Output is in form of vektor: [HR, Hz, H].

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)])

)$