% Imaging in metal cylinder
% cylinder and line charges in z-direction (in/out of screen).

close all
a = 0.1e-3; % radius or cylindrical PEC.
plot(a*cos(linspace(0,2*pi,100)), a*sin(linspace(0,2*pi,100)))
axis equal
hold on
d = 4*a; % Original line charge at 

% \overrightarrow{\;r^{\prime } } =\hat{\;x} d

plot(d, 0, 'r*')

di = a^2/d;
plot(di, 0, 'b*')
grid on

% Solve fields
x = 0.5*a:a/10:d-a/10;
y = -2*a:a/10:2*a;
[X,Y] = meshgrid(x,y);

theta = atan2(Y,X);
r = sqrt(X.^2 + Y.^2);

V = log(sqrt(  ((r.*cos(theta)-di).^2 + (r.*sin(theta)).^2)   ./ (( r.*cos(theta)-d).^2 + (r.*sin(theta)).^2)  ));
[Ex, Ey] = gradient(-V);

%%V(V<1) = -1
%V(V>1) = 1

quiver(x,y,Ex,Ey)

contour(X,Y,V)
for theta = linspace(pi/2+0.01, 3*pi/2-0.01, 7)
    streamline(x,y,Ex,Ey, d-1*a/10+a/10*cos(theta), a/10*sin(theta))
end
