%
%% 6.2 Helmholtz coil (Computational example)
% A Helmholtz coil arrangement consists of two current-carrying conductor rings in the same direction on a common axis. For an optimal homogeneous magnetic field between the coils, we must choose the coil radius and distance equal$({R_0} = D)$.
%
% We calculate the magnetic field lines in the x-z plane, plot them, and compare them to the experiment.
% Vector potential: \[{\vec A\left( {\vec R} \right) = \frac{1}{c}\int\limits_{}^{} {\vec j\left( {\vec r} \right)\frac{1}{{\left| {\vec R - \vec r} \right|}}{d^3}r} }\]
% Current: \[{\vec j\left( {\vec r} \right) = J{{\hat e}_\varphi }\delta \left( {r - {R_0}} \right)\left( {\delta \left( {z - \frac{{{R_0}}}{2}} \right) + \delta \left( {z + \frac{{{R_0}}}{2}} \right)} \right)}\]
% Parameter: R_0=J=c=1
%

clear all
syms z  J G jv
syms r R R0 Z phi Phi  positive
Par=[R0==1 J==1]

%
%% 1 Basic vectors in cylinder coordinates
% \[\vec r = \left( {\begin{array}{*{20}{c}}
% {r\cos \varphi }\\
% {r\sin \varphi }\\
% z
% \end{array}} \right)\], \[\vec R = \left( {\begin{array}{*{20}{c}}
% {R\cos \Phi }\\
% {R\sin \Phi }\\
% Z
% \end{array}} \right)\], \[{\hat e_\varphi } = \frac{1}{r}\frac{\partial }{{\partial \varphi }}\vec r\]

rv=[r*cos(phi);r*sin(phi);z]
Rv=[R*cos(Phi);R*sin(Phi);Z]
e_phi=diff(rv,phi)/r

%
%% 2 Vector potential
% Vector potential:  \[{\vec A\left( {\vec R} \right) = \frac{1}{c}\int\limits_{}^{} {\vec j(\vec r)\frac{1}{{\left| {\vec R - \vec r} \right|}} rdrd\varphi dz
% Current   \[\vec j\left( {\vec r} \right) = J{{\hat e}_\varphi }\delta \left( {r - {R_0}} \right)\left( {\delta \left( {z - \frac{{{R_0}}}{2}} \right) + \delta \left( {z + \frac{{{R_0}}}{2}} \right)} \right)\]

jv=J*e_phi*dirac(r-R0)*(dirac(z-R0/2)+dirac(z+R0/2))

% Propagator G=\frac{1}{|\vec r-\vec R|}

G=1/simplify(norm(rv-Rv))

% Potential without \varphi integration \vec A\!^{(\varphi )}=\int\vec j G rdrdz,\ \vec A=\int \vec A\!^{(\varphi )}d\varphi

Av_phi=int(int(r*jv*G,r,0,inf),z,-inf,inf)

% in the x-z plane (\Phi=0) and with parameters:

Av_phi=sube(Av_phi,[Phi==0 Par])

%
%% 3   x and z component of the vector potential
% \[{A_x} = \int\limits_0^{2\pi } {{{{A_1}^{\left( \varphi  \right)}}d\varphi } \],  \[{A_z} = \int\limits_0^{2\pi } {{{{A_3}^{\left( \varphi  \right)}}d\varphi } \]

A_x=int(Av_phi(1),phi,0,2*pi)
A_z=int(Av_phi(3),phi,0,2*pi)

%% 4    y component A_y with field lines (see also appendix below)
% \[{A_y} = \int\limits_0^{2\pi } {{{{A_2}^{\left( \varphi  \right)}}d\varphi } \]

[Rn,Zn,Ayn]=GetAy(Av_phi(2),24);
subplot(1,2,1)
hold on
surfc(Rn,Zn,Ayn);
axis([-4 4 -4 4 -6 6])
view([-120.2 25.0])

[Rn,Zn,Ayn]=GetAy(Av_phi(2),100);
subplot(1,2,2)
contour(Zn,Rn,Ayn,6*linspace(-1,1,25))
axis equal
axis off

%
%% Comparison with the experiment:
% Helmholtz coil and iron filings on a glass plate
%
% _________________________________________________________________________
%
%% Appendix: Field lines for Helmholtz coil in the z-x plane
% Vector potential: \[{\vec A\left( {\vec R} \right) = \frac{1}{c}\int\limits_{}^{} {\vec j\left( {\vec r} \right)\frac{1}{{\left| {\vec R - \vec r} \right|}}{d^3}r} }\]
% Current: \[{\vec j\left( {\vec r} \right) = J{{\hat e}_\varphi }\delta \left( {r - {R_0}} \right)\left( {\delta \left( {z - \frac{{{R_0}}}{2}} \right) + \delta \left( {z + \frac{{{R_0}}}{2}} \right)} \right)}\]
% Helmholtz coil and current components j:
%
% From symmetry considerations follows {A_x} = {A_z} = 0:
%
% The magnetic field lines in the z-x plane can be represented as contour lines of the vector potential {A_y}:
%
% _________________________________________________________________________
%
%% Function for \varphi-integration (numerical)

function [Rn,Zn,Ayn]=GetAy(Ay_phi,N)
[Rn,Zn]=meshgrid(linspace(-4,4,N));
Ay_phim=matlabFunction(Ay_phi);
Ayn=integral(@(phi)Ay_phim(Rn,Zn,phi),0,2*pi,'ArrayValued',true);
end