Matlab code to calculate the potential of line charges

Contents

• set up arrays of line charge locations
• strength of each line charge
• plot the Efield arrows and labels for the voltage contours

% Matlab code to calculate the potential due to a set of line charges.  The voltage
% reference is taken as the origin.  The code will have problems if any line
% charges are located at the origin.

% Clear the screen
clc;
clear all;
% close all windows
close all;

set up arrays of line charge locations
xcharge = [3.5 2 -2];  % x position of line charge in mm
ycharge = [-0.1 -1 3];  % y position of line charge in mm
xcharge = xcharge * .001; % convert to meters
ycharge = ycharge * .001;


strength of each line charge
rho_l_basic = 1e-7;   % units are C/m
rho_l = rho_l_basic * [1 3 -2];

% Set up grid for plotting
Ngrid = 33;
gridmax = .004;
% grid is the x and y values at which the calculation is done
grid = linspace(-gridmax, gridmax, Ngrid);

% Calculate voltage on grid due to each charge
% create 2 3D grids.  size of the arrays are Ngrid x Ngrid x ncharge
% The first array represents all combination of xgrid and ygrid
% locations with the xposition of the line charge.  The second
% represents all combinations with the y position of the line charges.
[xx1, yy1, xc1] = meshgrid(grid, grid, xcharge);
[xx2, yy2, yc2] = meshgrid(grid, grid, ycharge);
% Distance between the line charges and the grid points (3D array)
dist1 = sqrt( (xx1-xc1).^2 + (yy2 - yc2).^2);
% Grid points very close to a line charge will have a very high voltage that
%  will dominate contour plot.  Suppress this.
  dist1min =  .0002;
  dist1 = max(dist1, dist1min);
%  distance to reference voltage location (the origin).
dist2 = sqrt(xc1.^2 + yc2.^2);

% Calculate voltage at each grid point.  Vparts is the partial contribution to
% the voltage at each grid point (except for a charge strength factor).  It is
% an(Ngrid x Ngrid x ncharge) array.  The loop to calculate Vgrid then sums up
% the contributions of the various line charges.  Vgrid is an Ngrid x Ngrid
% array.
fac = 1./(2*pi*8.854e-12); % 1/(2*pi* epsilon0)
Vparts = fac.* log(dist2./dist1);
ncharge = length(xcharge); % get number of charges from length of array
Vgrid = zeros(Ngrid, Ngrid); % initialize Vgrid array
for ii = 1:ncharge
  Vgrid = Vgrid + Vparts(:,:,ii).*rho_l(ii);
end
[Ex, Ey] = gradient(-Vgrid);  % calculate electric field everywhere





plot the Efield arrows and labels for the voltage contours
% Plot results
% create position arrays in mm
grid_mm = grid*1000;
figure(1)
[cc, hh] = contour(grid_mm, grid_mm, Vgrid,15);

% adds contour value labels the voltage contours
clabel(cc,hh);

axis('equal');
axis([-4 4 -4 4]);
hold on;

% plot E field vectors
quiver(grid_mm, grid_mm, Ex, Ey);

title('Voltage contours and E field arrows');
ylabel('y (mm)');
xlabel('x (mm)');
hold off;

figure(2)
grid_cen = fix((Ngrid+1)/2);
plot(grid_mm, Vgrid(grid_cen,:));
title('Voltage on x axis');
ylabel('Voltage');
xlabel('x (mm)')