%% Method of Images
% Here we find the electrostatic fields for a point charge above an infinite PEC plane.
% This is easy once we find the expression for two charges mirrored in the plane.
% We then derive the E-field as -\nabla V


close all

h=1;

x = -1:0.05:1;
z = (0:0.05:0.8)';

[X,Z] = meshgrid(x,z);
V = 1./(sqrt(X.^2+(Z-h).^2)) - 1./(sqrt(X.^2+(Z+h).^2));


v = 1./(sqrt(x.^2+(z-h).^2)) - 1./(sqrt(x.^2+(z+h).^2));
subplot(211)
contour(x,z,v)
[Ex,Ez] = gradient(v);

hold on

quiver(x,z,Ex,Ez)
axis equal

for xx = -1:0.2:1 % Generate field lines for E field.
   streamline(x,z,Ex,Ez,xx,0) 
end

xlabel('x')
ylabel('z | V and E')
title('Point charge above PEC plane, point charge omitted due to singularity, \rho_s(x) overlaid')
plot(x,Ez(end,:)./5) % Overlay surface charge density for prettification.
xlim([-1 1])
ylim([0 0.8])

% Using the boundary condition for dielectric media \hat{n_2 } \cdot \left(\overrightarrow{\;D_1 } -\overrightarrow{\;D_2 } \right)=\rho_s
% we get \rho_s =-\varepsilon_0 \frac{\partial \;V}{\partial \;z} evaluated at z=0.
% This is simply E_z =-\hat{\;z} \cdot \nabla V which we have already computed.

subplot(212)
plot(x,Ez(end,:))
xlabel('x')
ylabel('\rho_s')
title('Surface charge density magnitude')