# Help with creating contour plot

24 views (last 30 days)
Andrew on 4 Jun 2019
Edited: Andrew on 8 Jun 2019
I am trying to graph the electric potential lines in a plane perpendicular to a ring of charge in the x-y plane. Deriving the formula for off axes points is tedious so I set up a the integrand and then summed it up over a discretized interval in the form of a reimann sum using a for loop. I know my expression for V (Vtot) is correct because it approximately matches analytical solutions for V and -grad(Vtot) appears as it should in the quiver plot to the right below. The problem is probably somewhere in my substitution for the symbols x,y and z (perhaps the meshgrid).
%% Computing a symbolic expression for V for anywhere in space
syms x y z % phiprime is angle that an elemental dq of the circular charge is located at, x,y and z are arbitrary
% points in space outside the charge distribution
N = 200; % number of increments to sum
R = 2; % radius of circle is 2 meters
dphi = 2*pi/N; % discretizing the circular line of charge which spans 2pi
integrand = 0;
for phiprime = 0:dphi:2*pi %phiprime ranges from 0 to 2pi in increments of dphi
integrand = integrand + dphi./(sqrt(((x - R.*cos(phiprime) )).^2 + ((y - R.*sin(phiprime) ).^2) + z.^2));
end
integral = sum(integrand); % uncessary but harmless step that I leave to show that I am using the
% sum of the above expression for each dphi
eps0 = 8.854e-12;
kC = 1/(4*pi*eps0);
rhol = 1*10^-9; % linear charge density
Vtot = kC*rhol*R.*integral; % symbolic expression for Vtot
%% Graphing V & E in plane perpedicular to ring & passing through center
[Y1, Z1] = meshgrid(-4:.5:4, -4:.5:4);
Vcont1 = subs(Vtot, [x,y,z], {0,Y1,Z1}); % Vcont1 stands for V contour; 1 is becuase I do the plane of the ring next
contour(Y1,Z1,Vcont1)
xlabel('y - axis [m]')
ylabel('z - axis [m]')
title('V in a plane perpedicular to a ring of charge (N = 200)')
str = {'Red line is side view', 'of ring of charge'};
text(-1,2,str)
hold on
circle = rectangle('Position',[-2 0 4 .1],'Curvature',[1,1]); % visually displaying line of charge on plot
set(circle,'FaceColor',[1, 0, 0],'EdgeColor',[1, 0, 0]);
g = gradient(-1.*(kC*rhol*R.*integral),[x,y,z]); % taking negative gradient of V and finding symbolic equations
% for Ex, Ey and Ez
% now substituting all the values of the 2D coordinate system for the symbolic x and y variables to get
% numeric values for Ex and Ey because the gradient command can't differentiate a scalar
Ey1 = subs(g(2), [x y z], {0,Y1,Z1});
Ez1 = subs(g(3), [x y z], {0,Y1,Z1});
E1 = sqrt(Ey1.^2 + Ez1.^2); % full numeric magnitude of E in y-z plane
Eynorm1 = Ey1./E1; % This normalizes the electric field lines
Eznorm1 = Ez1./E1;
quiver(Y1,Z1,Eynorm1,Eznorm1);
hold off
***NOTE FOR BELOW FIGURES: The axes are labeled wrong. My mistake! ***** The x axis is (y) and the y axis is (z) ******  As you can see, my contour plot is not filling space. When I lower the value of N, suddenly my contour plot begins to fill space as seen below which makes no sense. N should be totally unrelated; it only determines the accuracy of Vtot (NOTE: I changed the increments of the meshgrid of the plot below, that's why the vectors are less dense). I've used the method of symbols and a for loop and then substitution for a LINE of charge and everything worked perfectly. I don't know why this contour plot is being uncooperative for a ring of charge. Thanks in advance!

Andrew on 8 Jun 2019
OOOOOOHHHHH. I understand what you mean now. I didn't realize you could pick out which lines were plotted by using another vector defined as a linspace and plug that vector into the 'levels' spot of the contour() command.
Now that I re-read the section on 'contour' I realize what "To draw the contour lines at specific heights, specify levels as a vector of monotonically increasing values" means.
Wow you don't how much you've helped me. I've spent hours trying to figure out ways around this and the solution was right under my nose. Thank you Walter. And thanks KSSV. Really appreciate and wish there was some way to repay you besides accepting the answer (which I don't know how to do by the way, sorry I'm new).
Not that you care but here is my beautiful contour plot now: That's exactly what I needed.
I don't know if I should ask this or create a new question but,
is there a way I could make a contour plot of a riemann sum WITHOUT using syms and the 'subs' command to sub values for the variables which are acting as constants in the riemann sum?
Essentially, what I already did here but a method which avoids the 'subs' command since it's my understanding that 'subs' takes a long time to execute for many values.
Walter Roberson on 8 Jun 2019
Not at the moment, but you can improve performance with
Vcont1 = double( subs( vpa(Vtot), [x,y,z], {0,Y1,Z1}) );
Andrew on 8 Jun 2019
I've never heard of Variable-precision arithmetic. After some cursory research, it seems very interesting and useful; and you expressed it in terms of my code. Wow. THANK YOU. Thank you so much. Your responses were prompt (not that they needed to be) and very helpful. Thanks and have a nice day!