Calculating Electric dipole using Coulomb's law with Quiver

Below is the code block that calculates Electric Fields due to dipoles in a 2 dimensional plane using Coulomb's Law. I need 'Quivers' to be drawn on this code.

%-------------------------------------------------------------------------% 
%
%     Computing a two-dimensional sketch of the equipotential lines and 
%     electric field lines of an electric dipole consisting of a pair of 
%     opposite charges of magnitude Q using Coulomb's Law.
%
%-------------------------------------------------------------------------%
clear; 
close all; 
clc;
%-------------------------------------------------------------------------%
%                             INITIALIZATION                                  
%-------------------------------------------------------------------------%
% Constant 1/(4*pi*epsilon_0) = 9*10^9
k = 9*10^9;
% eps_r = Relative permittivity
eps_r = 1;
charge_order = 10^-9; % milli, micro, nano etc..
const = k*charge_order/eps_r;
% Nx = Number of grid points in X- direction
% Ny = Number of grid points in Y-Direction
Nx = 101; % For 1 meter
Ny = 101; % For 1 meter
% n = Number of charges
n = 2;
% Electric fields Initialization
% E = Total electric field
E_f = zeros(Nx,Ny);
% Ex = X-Component of Electric-Field
% Ey = Y-Component of Electric-Field
Ex = E_f;
Ey = E_f;
% ex = unit vector for x-component electric field
% ey = unit vector for y-component electric field
ex = E_f;
ey = E_f;
% r = distance between a selected point and the location of charge
r = E_f;
r_square = E_f;
% Array of charges
% Q = All the 'n' charges are stored here
Q = [1,-1];
% Array of locations
X = [10,-10];
Y = [0,0];

%-------------------------------------------------------------------------%
%                      COMPUTATION OF ELECTRIC FIELDS
%-------------------------------------------------------------------------%
%  Repeat for all the 'n' charges
for k = 1:n
    q = Q(k);
    
    % Compute the unit vectors
    for i=1:Nx
        for j=1:Ny
            r_square(i,j) = (i-51-X(k))^2+(j-51-Y(k))^2;
            r(i,j) = sqrt(r_square(i,j));
            ex(i,j) = ex(i,j)+(i-51-X(k))./r(i,j);
            ey(i,j) = ey(i,j)+(j-51-Y(k))./r(i,j);
        end
    end
    
    E_f = E_f + q.*const./r_square;
    Ex = Ex + E_f.*ex.*const;
    Ey = Ex + E_f.*ey.*const;
end
%-------------------------------------------------------------------------%
%                           PLOT THE RESULTS
%-------------------------------------------------------------------------%

x_range = (1:Nx)-51;
y_range = (1:Ny)-51;
contour_range = -10:0.01:10;
contour(x_range,y_range,E_f',contour_range,'linewidth',0.5);
axis([-15 15 -15 15]);
xlabel('x','fontsize',14);
ylabel('y','fontsize',14);
title('Electric Dipole Distribution','fontsize',14);

Ultimately, I should get a result like the images above. I looked at how it was written about 'Quiver' from this site, but I couldn't implement it to this code.

Thank you in advance for your help.

[William Rose]

@Kutlu Yigitturk, What part does not work - the quivrs plot, or the contour plot, or both? Ar the plots above taken from another source? If they are from your program, what is wrong with them?

[William Rose]

@Kutlu Yigitturk The contour plot looks good. Quiver plot - wait a sec.

Answer 1

[William Rose]

@Kutlu Yigitturk,

I modified the program a bit and added the calculation of voltage as well as electric field strength.

I think part of the problem before was that the total E field value was being updated inside a loop, in a place where adding up vector (Ex, Ey) was reasonable, but adding up magnitudes was not. Now I calculate Ex and Ey inside the k=1:n loop, then I calculate E_f outside that loop, using the array of Ex and the array of Ey.

I also added a surface plot of voltage because it is kind of fun to see a 3D plot.

[Kutlu Yigitturk]

Thank you for your help. These very well for me.

[William Rose]

You're welcome, @Kutlu Yigitturk.

Answer 2

[William Rose]

@Kutlu Yigitturk

[Added this later: the code I attached to this message has one or more mistakes. Please ignore it. Sorry.]

mistake that party accounts for the bad plotThe attached code makes a quiver plot as well as the contour plot. (Plot attached.) The contour plot looks fine but the quver plot does not. Are you sure your equations for ex and ey are correct?

I notice you had to transpose the E_f array when making the contour plot. Why? That suggests to me that the data in E_f was not oriented the way you expected.

Is E_f supposd to be the potential or the magnitude of the electric field at each point?

Try making Q(2)=0. Then the contour lines and arrows should like like a monopole field. The contour plot does look correct in this case, but the quiver plot does not look correct. Maybe that approach will help you figure out f there is an error in the equations for ex and ey.