% T-3.25
% E and V in the xy-plane for a cylinder with radius a=1 with surface charge Q, oriented along the z-axis.
% V defined as 0 when r=b>a.
% Cylinder is a PEC, so E=0 within it, and V=constant
% The V constant comes from:
% V\left(\overrightarrow{\;r} \right)=-\int_{L_{\textrm{ground}} }^{\overrightarrow{\;r} } \overrightarrow{E} \cdot {\textrm{dl}}^{\prime }

EM_constants

e0 = eps0;
q = q0;
L = 1;

a = 1;
b = 10;
r = 0:0.01:b;
mask = [zeros(1,a/0.01) ones(1,(b-(a-0.01))/0.01)] % Conductivity mask

% \overrightarrow{E} =\frac{Q}{2\pi \epsilon_0 L\;r}\hat{r}
% V=\frac{Q}{2\pi \epsilon_0 L}\ln \left(\frac{b}{r}\right)\;\;\;\;,r>a

E =  q./(2*pi*e0*L*r(r>a));
V1 = q./(2*pi*e0*L)*(log(b/a));
V2 = q./(2*pi*e0*L).*(log(b./r(r>a)));


figure
hold on
plot(r(r<=a), repmat(0, [1, length(r(r<=a))]), 'r') 

idx = find(r>a, 1,'first');
plot([a r(idx)], [0 q./(2*pi*e0*L*r(idx))], 'r') % just ignore this shit
plot(r(r>a), E, 'r')
title('E-field')

figure 
hold on
plot(r(r<=a), repmat(V1, [1, length(r(r<=a))]), 'k')
plot(r(r>a), V2, 'k')
title('Electrostatic Potential V')