% Plot electric vector field and potential scalar field between two point
% charges.
close all

q0 = 1.60217662e-19; % Elementary charge
eps0 = 8.8541878128e-12

Q = [1 -1].*(1/(4*pi*eps0))*q0
q1 = [-0.5 0];
q2 = [0.5 0];
q = [q1;q2];

n = 7
h = 0.1;
[X Y] = meshgrid((-n:h:n), (-n:h:n));

rx = X % Generate coordinates
ry = Y

Ex = zeros(size(X));
Ey = zeros(size(Y));

V = zeros(size(X));

for i = 1:2
    rsourcex = q(i,1); % Coordinate of point charge.
    rsourcey = q(i,2);

    sourcetofieldx = rx - rsourcex % Vector from point charge to each field point.
    sourcetofieldy = ry - rsourcey

    sourcetofieldMAG = sqrt(sourcetofieldx.^2 + sourcetofieldy.^2); % Distance between point charge and each field point.


    Ex = Ex + Q(i)./(sourcetofieldMAG.^2).*(sourcetofieldx./sourcetofieldMAG); % Coulomb's law
    Ey = Ey + Q(i)./(sourcetofieldMAG.^2).*(sourcetofieldy./sourcetofieldMAG);
    
    V = V + Q(i)./sourcetofieldMAG.^2;
end


figure(1)
s=pcolor(X,Y,V);
        s.FaceColor = 'interp';
        s.EdgeColor = 'none';
        view(2);
        colorbar
        caxis([-5e-11 5e-11])
        axis square
        %title('Voltage')
hold on


figure(1)
%quiver(X,Y,Ex,Ey,2)
hold on
plot(q(1,1),q(1,2), 'r*')
plot(q(2,1),q(2,2), 'r*')
axis square
grid on

hold on

%plot(0.2*cos(0:2*pi/10:2*pi)+0.5, 0.2*sin(0:2*pi/10:2*pi), 'g*')

hlines = streamline(X,Y,Ex,Ey, 0.2*cos(0:2*pi/20:2*pi)+0.5, 0.2*sin(0:2*pi/20:2*pi))
set(hlines,'LineWidth',0.5,'Color','k')
hlines = streamline(X,Y,-Ex,-Ey, 0.2*cos(0:2*pi/20:2*pi)+0.5, 0.2*sin(0:2*pi/20:2*pi))
set(hlines,'LineWidth',0.5,'Color','k')

hlines = streamline(X,Y,Ex,Ey, 0.2*cos(0:2*pi/20:2*pi)-0.5, 0.2*sin(0:2*pi/20:2*pi))
set(hlines,'LineWidth',0.5,'Color','k')
hlines = streamline(X,Y,-Ex,-Ey, 0.2*cos(0:2*pi/20:2*pi)-0.5, 0.2*sin(0:2*pi/20:2*pi))
set(hlines,'LineWidth',0.5,'Color','k')

ax = gca;
exportgraphics(ax, 'dipole.pdf','ContentType','vector')


%streamlines in circle around perimeter of system
% 


%streamparticles(X,Y,Ex,Ey, 0.4*cos(0:2*pi/20:2*pi)+0.5, 0.4*sin(0:2*pi/20:2*pi))
