%
%% 6.3 Boundary value problem: charge over grounded sphere (Computational example)
% A charge Q is located above a grounded metal sphere. The electrical potential can be calculated with the help of a mirror charge q (Method of Images).
%
%                             .
%
% Parameter: Q=1,\ P=3.5

syms r theta P p Q q x z real
Par=[Q==1,P==3.5]

%
%% 1 Coordinates
% Position vector in the x-z plane (\varphi=0): \[\vec r = \left( {\begin{array}{*{20}{c}}
% {r\sin \theta }\\
% 0\\
% {r\cos \theta }
% \end{array}} \right)\]

r_v=[r*sin(theta);0;r*cos(theta)]

% Positions of charge Q and mirror charge q
% with \[{\hat e_z} = \left( {\begin{array}{*{20}{c}}
% 0\\
% 0\\
% 1
% \end{array}} \right)\]:
% Position of charge Q: \vec P=P\hat e_z
% Position of mirror charge q: \vec p=p\hat e_z

e_z=[0;0;1];
P_v=P*e_z
p_v=p*e_z

%
%% 2 Contributions to the potential
% Potential of the charge: \varphi_Q(\vec r)=\frac{Q}{\left|\vec r-P\hat e_z\right|}

phi_Q=simplify(Q/norm(P_v-r_v))

% Potential of the mirror charge \varphi_q(\vec r)=\frac{q}{|\vec r-p\hat e_z|}

phi_q=simplify(q/norm(p_v-r_v))

% Total potential: \phi=\phi_Q+\phi_q\ (\forall r>1)

phi=phi_Q+phi_q

%
%% 3 Boundary conditions
% Boundary conditions:  \left \varphi \right|_{r=1}=0  \to\ \ p(Q,P),\ q(Q,P)

sube(phi,r==1)==0
isolate(ans,cos(theta))

% For this equation to hold for all \theta\in[0,\pi], the numerator and denominator of the right side must be zero!

[A,B]=numden(rhs(ans));
assume(Q>0 & q<0 & P>1 & p<1)
solve([A,B],[q,p]);
[p==ans.p,q==ans.q]
Phi=simplify(sube(phi,ans))

%
%% 4 Potential in Cartesian coordinates
% Transformation: \cos\theta ={z\over r},\ r=\sqrt{x^2+z^2}

Phi_xy=sube(Phi,[Par,cos(theta)==z/r r==sqrt(x^2+z^2)])

%
%% 5 Contour lines of the potential

N=100;
[XX,ZZ]=meshgrid(linspace(-7,7,N),linspace(-7,7,N));
F=matlabFunction(Phi_xy)
contour(XX,ZZ,F(XX,ZZ),linspace(0,0.75,100))
axis equal
set(0,'defaultaxesfontsize',14,'defaultaxeslinewidth',1)
xlabel('Position x')
ylabel('Position z')
title('Contour lines of \phi(x,0,z)')

%
%% 6 Electric field lines
% Electric field in the presentation plane  $\left( \matrix{
%   {E_x}  \cr
%   {E_y}  \cr}  \right) = -\left( {\matrix{
%    {{{\partial \phi } \over {\partial x}}}  \cr
%    {{{\partial \phi } \over {\partial y}}}  \cr
%
%  } } \right)$

Ex=matlabFunction(-diff(Phi_xy,x));
Ez=matlabFunction(-diff(Phi_xy,z));
r0=0.2;
phi0=linspace(0,2*pi,50);
h=streamline(XX,ZZ,Ex(XX,ZZ),Ez(XX,ZZ),...
    r0*sin(phi0),r0*cos(phi0)+double(sube(P,Par)));
set(h,'Color','r','LineWidth',1)
rectangle('Position',[-1 -1 2 2],'FaceColor',[.8 .8 1],'Curvature',[1 1])
title('Contour lines and field lines of \phi(x,0,z)')

%