% Plots the potential and electric field created by two spherical 
% and uniform charge distributions with the following total
% charges, radii, and positions (in MKS units):
%   Q1=-2.e-9   R1=0.08   [x1,y1,z1]=[-1,0,0]
%   Q2=+1.e-9   R2=0.05   [x2,y2,z2]=[+1,0,0]
% Notice that the potential of a uniform spherical charge is
%   V(r) = K*Q*(3*R^2-r^2)/(2*R^3)  inside  (r<R)
%        = K*Q/r                    outside (r>R)
% J.M.Soler. Sept.2010 & Sep.2021

% Set parameters, in MKS units
clear all
K_coulomb=9.e9;                     % Coulomb constant in J*m^2/C^2
Q1=-2.e-9; R1=0.08; x1=-0.15; y1=0; z1=0; % 1st charge, radius, and position
Q2=+1.e-9; R2=0.05; x2=+0.15; y2=0; z2=0; % 2nd charge, radius, and position
dx=0.02; dy=dx; dz=dx;              % Interval between plot points
xmin=-0.35; xmax=0.35;              % Plot range in x axis
ymin=-0.25; ymax=0.25;              % Plot range in y axis
zmin=-2*dz; zmax=2*dz;              % Plot range in z axis

% Set plot points along x and y axes
x = xmin : dx : xmax;
y = ymin : dx : ymax;
z = zmin : dz : zmax;
nx = numel(x);
ny = numel(y);
nz = numel(z);

% Find index at which z=0
iz0 = find(z==0);

% Convert plot points into matrices
% [x,y,z] = meshgrid(x,y,z);        % Produces arrays x,y,z with order x(iy,ix,iz)
[x,y,z] = ndgrid(x,y,z);            % Produces arrays x,y,z with order x(ix,iy,iz)
nPoints = numel(x);                 % Total number of points

% Find electrostatic potential, using for loops
% V = zeros(size(x));                            % Create matrix for V values
% for i = 1:nPoints                              % For every plot point
%    d1 = sqrt( (x(i)-x1)^2 + (y(i)-y1)^2 + (z(i)-z1)^2 );  % Distance to 1st charge
%    d2 = sqrt( (x(i)-x2)^2 + (y(i)-y2)^2 + (z(i)-z2)^2 );  % Distance to 2nd charge
%    V1 = sphericalChargePotential(Q1,R1,d1);    % Potential of 1st charge
%    V2 = sphericalChargePotential(Q2,R2,d2);    % Potential of 2nd charge
%    V(i) = V1 + V2;                             % Total potential at point i
% end

% Alternative: find electrostatic potential for all points at once
d1 = sqrt( (x-x1).^2 + (y-y1).^2 + (z-z1).^2 );  % Distance to 1st charge
d2 = sqrt( (x-x2).^2 + (y-y2).^2 + (z-z2).^2 );  % Distance to 2nd charge
V1 = sphericalChargePotential(Q1,R1,d1);         % Potential of 1st charge
V2 = sphericalChargePotential(Q2,R2,d2);         % Potential of 2nd charge
V = V1 + V2;                                     % Total potential

% Surface plot of the potential
figure(1)
% surfc(x(:,:,iz0),y(:,:,iz0),V(:,:,iz0))                              % with meshgrid
surfc(x(:,:,iz0)',y(:,:,iz0)',V(:,:,iz0)')                             % with ndgrid
xlabel('x (m)'), ylabel('y (m)'), zlabel('V (V)')
title('Potential of two uniform spherical charges')

% Find the electric field E = -gradient(V)
[Ex,Ey,Ez] = gradient(-V,dx,dy,dz);

% Plot a contour map of V and a quiver (arrows) map of E
figure(2)
% contour(x(:,:,iz0),y(:,:,iz0),V(:,:,iz0))
contour(x(:,:,iz0)',y(:,:,iz0)',V(:,:,iz0)')
axis equal
xlabel('x (m)'), ylabel('y (m)')
hold on
arrowLength = 1.5;    % max. arrow length, in units of mesh interval
% quiver(x(:,:,iz0),y(:,:,iz0),Ex(:,:,iz0),Ey(:,:,iz0),arrowLength)    % with meshgrid
quiver(x(:,:,iz0)',y(:,:,iz0)',Ex(:,:,iz0)',Ey(:,:,iz0)',arrowLength)  % with ndgrid
title('Potential and electric field of two uniform spherical charges')
hold off

% Find rho = -laplacian(V)/(4*pi*Kc) = (dEx/dx + dEy/dy) / (4*pi*Kc)
[dExdx,dExdy,dExdz] = gradient(Ex,dx,dy,dz);
[dEydx,dEydy,dEydz] = gradient(Ey,dx,dy,dz);
[dEzdx,dEzdy,dEzdz] = gradient(Ez,dx,dy,dz);
rho = (dExdx + dEydy + dEzdz) / (4*pi*K_coulomb);

% Plot the density (in nC/m^3 to avoid an apparent bug in surfc)
figure(3)
% surfc(x(:,:,iz0),y(:,:,iz0),rho(:,:,iz0)/1.e-9)                      % with meshgrid
surfc(x(:,:,iz0)',y(:,:,iz0)',rho(:,:,iz0)'/1.e-9)                     % with ndgrid
xlabel('x (m)'), ylabel('y (m)'), zlabel('rho (nC/m^3)')
title('Density of two uniform spherical charges')

% Find the density from the laplacian of the potential using
% d2V/dx2 + d2V/dy2 + d2V/dz2 = -4*pi*Kc*rho
ix=2:nx-1; iy=2:ny-1; iz=2:nz-1;
% Notice that, with meshgrid, first (row) index refers to y, 
%   and second (column) index refers to x
% rho(iy,ix,iz) = ...                                                  % with meshgrid
%    -( (V(iy+1,ix,iz)-2*V(iy,ix,iz)+V(iy-1,ix,iz))/dy^2   ...
%      +(V(iy,ix+1,iz)-2*V(iy,ix,iz)+V(iy,ix-1,iz))/dx^2   ...
%      +(V(iy,ix,iz+1)-2*V(iy,ix,iz)+V(iy,ix,iz-1))/dz^2 ) ...
%     / (4*pi*K_coulomb);
% Natural order index using ndgrid:
rho(ix,iy,iz) = ...                                                    % with ndgrid
   -( (V(ix+1,iy,iz)-2*V(ix,iy,iz)+V(ix-1,iy,iz))/dx^2   ...
     +(V(ix,iy+1,iz)-2*V(ix,iy,iz)+V(ix,iy-1,iz))/dy^2   ...
     +(V(ix,iy,iz+1)-2*V(ix,iy,iz)+V(ix,iy,iz-1))/dz^2 ) ...
    / (4*pi*K_coulomb);

% Plot the density again, calculated with the Laplacian of the potential
figure(4)
% surfc(x(iy,ix,iz0),y(iy,ix,iz0),rho(iy,ix,iz0)/1.e-9)                % with meshgrid
surfc(x(ix,iy,iz0)',y(ix,iy,iz0)',rho(ix,iy,iz0)'/1.e-9)               % with ndgrid
xlabel('x (m)'), ylabel('y (m)'), zlabel('rho (nC/m^3)')
title('Density of two uniform spherical charges')
