%
%% 6.4 Radiation of the Hertzian dipole (Computational example)
% A Hertzian dipole consists of a negative point charge at the coordinate origin and a harmonic-moving positive point charge at $\vec a\left( t \right)$.
% Trajectory of the positive charge: \[\vec a\left( t \right) = {a_0}\cos \left( {\omega t} \right){{\hat e}_z}\]
% Charge density: \[{\rho \left( {\vec r,t} \right) = q\left( {\delta \left( {\vec r - \vec a\left( t \right)} \right) - \delta \left( {\vec r} \right)} \right)}\]
% Current: \[{\vec j\left( {\vec r,t} \right) = \vec v\left( {\vec r,t} \right)\rho \left( {\vec r,t} \right) = q\dot \vec a\left( t \right)\delta \left( {\vec r - \vec a\left( t \right)} \right){\kern 1pt} {\kern 1pt} }\]
% Vector potential: \[{\vec A\left( {\vec r,t} \right) = \frac{{{\mu _0}}}{{4\pi }}\int {dt'{d^3}\vec r'\frac{{\vec j\left( {\vec r',t'} \right)\delta \left( {t - {c^{ - 1}}\left| {\vec r - \vec r'} \right| - t'} \right)}}{{\left| {\vec r - \vec r'} \right|}}} }\]
% Parameter: \[{a_0} = q = \mu  = \omega  = c = 1\]

clear all
syms r t tb a0 q av positive
syms theta phi omega mu c x z real

Par=[a0==1 q==1 mu==1 omega==1 c==1]

%
%% 1 Spherical coordinates
% Spherical coordinates: \[\vec r = \left( {\begin{array}{*{20}{c}}
%   {r\sin \theta \cos \varphi } \\
%   {r\sin \theta \sin \varphi } \\
%   {r\cos \theta }
% \end{array}} \right)\]
% Jacobian: \[J = \left( {\frac{{\partial {x_i}}}{{\partial {q_j}}}} \right) = \left( {{{\vec g}_r},{{\vec g}_\theta },{{\vec g}_\varphi }} \right),\quad  \to \quad {\hat e_r} = \frac{{{{\vec g}_r}}}{{\left| {{{\vec g}_r}} \right|}},\;{\hat e_\theta } = \frac{{{{\vec g}_\theta }}}{{\left| {{{\vec g}_\theta }} \right|}},\;{\hat e_\varphi } = \frac{{{{\vec g}_\varphi }}}{{\left| {{{\vec g}_\varphi }} \right|}}\]
% z-direction: \[{\hat e_z} = \left( {\begin{array}{*{20}{c}}
%   0 \\
%   0 \\
%   1
% \end{array}} \right)\]

r_v=[r*sin(theta)*cos(phi);r*sin(theta)*sin(phi);r*cos(theta)]
J=jacobian(r_v,[r,theta,phi]);
e_r=J(:,1)
e_theta=J(:,2)/r
e_phi=J(:,3)/(r*sin(theta))
E_z_xA=[0;0;1]
XPlane=[phi==0 theta==pi/2 r==x]
XZPlane=[phi==0 sin(theta)==x/r cos(theta)==z/r r==sqrt(x^2+z^2)]

%
%% 2 Vector potential in far field approximation
% Current of the Hertzian dipole: \[\vec j\left( {\vec r,t} \right) = q\dot \vec a\left( t \right)\delta \left( {\vec r - \vec a\left( t \right)} \right)\mathop  \to \limits_{{\text{far field approximation}}} q\dot \vec a\left( t \right)\delta \left( {\vec r} \right)\]
% Vector potential: \[{\vec A\left( {\vec r,t} \right) = \frac{{{\mu _0}}}{{4\pi }}\int {dt'{d^3}\vec r'\frac{{\vec j\left( {\vec r',t'} \right)\delta \left( {t - {c^{ - 1}}\left| {\vec r - \vec r'} \right| - t'} \right)}}{{\left| {\vec r - \vec r'} \right|}}} }\]
%                                      \[ \to \quad \frac{{{\mu _0}}}{{4\pi }}\int {dt'{d^3}\vec r'\frac{{q\dot \vec a\left( {t'} \right)\delta \left( {\vec r'} \right)\delta \left( {t - {c^{ - 1}}\left| {\vec r - \vec r'} \right| - t'} \right)}}{{\left| {\vec r - \vec r'} \right|}}} \]
%                                         \[ = \frac{{{\mu _0}}}{{4\pi r}}q\dot \vec a\left( {t - r/c} \right) =  - \frac{{{\mu _0}}}{{4\pi r}}q\omega \sin \left( {t - r/c} \right){{\hat e}_z}\]

A_v=(-mu*omega*a0*q*sin(omega*(t-r/c))/(4*pi*r))*E_z_xA

% in spherical coordinates basis:  ${\vec A^{\left( {Ku} \right)}} = \left( {\matrix{
%    {\vec A \cdot {{\hat e}_r}}  \cr
%    {\vec A \cdot {{\hat e}_\vartheta }}  \cr
%    {\vec A \cdot {{\hat e}_\varphi }}  \cr
%
%  } } \right)$

A_Ku=[dot(A_v,e_r);dot(A_v,e_theta);dot(A_v,e_phi)]

%
%% 3 B and E field in spherical coordinates basis
% Magnetic field: \[{{{\vec B}^{\left( {Ku} \right)}} = rot\;{{\vec A}^{\left( {Ku} \right)}}}\]
% Electric field: \[{{{\vec E}^{\left( {Ku} \right)}} = \int\limits_{}^{} {rot\;{{\vec B}^{\left( {Ku} \right)}}dt} }\]
% The curl in spherical coordinates is calculated with the function Curl3d at the end of this script.

B_Ku=Curl3d(A_Ku,[r,theta,phi],[1,r,r*sin(theta)])
E_Ku=int(Curl3d(B_Ku,[r,theta,phi],[1,r,r*sin(theta)]),t)

%
%% 4 B and E on the x-axis
% x-axis: r=x,\ \theta=\frac{\pi}{2},\ \varphi=0
% E field: \vec E=
% E_r^{(Ku)}\hat e_r+E_\theta^{(Ku)}
% B field: \vec B=B_\varphi^{(Ku)}\hat e_\varphi

E_v=E_Ku(1)*e_r+E_Ku(2)*e_theta

% on the x axis

sube(E_v,XPlane)
E_z_xA=ans(3)

% y component of the B field on the x axis

B_v=B_Ku(3)*e_phi
sube(B_v,XPlane)
B_y_xA=ans(2)

%
%% 5 Graphics

R=linspace(0,20);
subplot(2,2,1)
hold on
plot3([0 20],[0 0],[0 0],':k')
B_y_xAM=matlabFunction(sube(B_y_xA,Par))
plot3(R,B_y_xAM(0,R),0*R)
E_z_xAM=matlabFunction(sube(E_z_xA,Par))
plot3(R,0*R,E_z_xAM(0,R))
axis([0 20 -.05 .05 -.05 .05])
view([-74 23])
axis square
xlabel('x axis')
ylabel('B_y axis')
zlabel('E_z axis')
hold off
subplot(2,2,2)
plot(R,B_y_xAM(0,R),R,-E_z_xAM(0,R),[0 20],[0 0],':k')
axis([0 20 -0.015 0.05])
xlabel('x axis')
axis square
legend('B_y','-E_z')

subplot(2,2,1)
view([-38.673 15.726])
subplot(2,2,2)
view([0.00 90.00])

%
%% 6 Field lines of the E field in the x-z plane
% The electric field lines can be plotted in the x-z plane as contour lines of the magnetic field integrated over time:
%

B_v
sube(B_v,XZPlane)
B_yt_xzE=int(ans(2),t)
B_yt_xzEM=matlabFunction(sube(B_yt_xzE,Par))

N=30
[R,The]=meshgrid(linspace(0,40,N),linspace(0,2*pi,N));
X=R.*sin(The);
Z=R.*cos(The);
subplot(1,1,1)
surfc(X,Z,B_yt_xzEM(0,X,Z))
axis square
xlabel('x axis')
ylabel('z axis')
zlabel('\int{ B_z(x,0,z) dt}')

%
%% 7 Animation of the field lines in the x-z plane (-> Hertz_Dipol.gif)

subplot(1,1,1);
%syms x y
fun=@(Ta)fcontour(B_yt_xzEM(Ta,x,z),40*[-1 1],'LevelList',0.002*linspace(-1,1,4));

%

fanimator(fun,'AnimationRange',[0 2*pi],'FrameRate',5);
hold on
fanimator(@(Ta)text(0,45,sprintf('Node lines of the E field in x-z plane, T=%2.1f',Ta),'HorizontalAlignment','center'),...
    'AnimationRange',[0 2*pi],'FrameRate',5);
axis square
xlabel('x')
ylabel('z')

%

writeAnimation('Hertz_Dipol.gif','Loopcount',inf)

%
%% Curl in generalized coordinates (function)

function Res=Curl3d(Ex,Ko,Me)    
% Curl in generalized coordinates
% (Expression F,coordinates q, metric g)

% \[{{\hat e}_1}\left( {\vec \nabla  \times \vec F} \right) = \frac{1}{{{g_2}{g_3}}}\left( {\frac{\partial }{{\partial {q^2}}}{g_3}{F_3} - \frac{\partial }{{\partial {q^3}}}{g_2}{F_2}} \right)\quad \left( {cyclical} \right),    {g_i} 	\equiv  \left| {\frac{{\partial \vec r}}{{\partial {q^i}}}} \right|, \quad     F_i	\equiv\vec F \cdot \hat e_i

    Res=[...
        1/(Me(2)*Me(3))*(diff(Me(3)*Ex(3,1),Ko(2))-diff(Me(2)*Ex(2,1),Ko(3)));...
        1/(Me(3)*Me(1))*(diff(Me(1)*Ex(1,1),Ko(3))-diff(Me(3)*Ex(3,1),Ko(1)));...
        1/(Me(1)*Me(2))*(diff(Me(2)*Ex(2,1),Ko(1))-diff(Me(1)*Ex(1,1),Ko(2)))...
        ];
    Res=simplify(Res);
end