Looking for bug in a graphics program for plotting dipole fields

3 次查看（过去 30 天）
David Gabriel
David Gabriel 2024-7-15
回答： David Gabriel 2024-7-25
Ran in:
Greetings,
First, I apologize for postiing in the "General" area. This is the fisrt time I am asking for help in the MATLAB community and navagting this tyupe of forum. In looking of a MATAB program to plot both electric and potential fields of a dipole I came across a Book Chapter under the Academia profile of Darvin Messi on Numerical Methods. See the following link:
https://www.academia.edu/7995677/NUMERICAL_METHODS?email_work_card=view-paper&li=0
I typed in the code and worked through a majority of bugs, a couple due to typos in the text. One had to do with adding a symbol to the plot function. I got the electric fields portion to work perfectly. However, the electric potential plots still do not work. I have tried to the best of my ability to error trap. I was able to get a couple of points to plot but nothing more. I really like this approach which does not make use of MATLAB's mesh or gradient functions because of the application I have in working with students. On the other hand, I do not know why this portion is not working. Any help would be greatly appreciated as I would not trouble the MATLAB community without exhausting the combination/permutations of what could be wrong.
I would be great to then keep the corrected version on this community as through my searches, a program like this has been requested by students very frequently.
Best wishes,
David.
%. Program below
plotit ( [-1 1], [-1.5 0; 1.5 0], 1, 1, 0.01, 0.01, 20, 20, 5)
FACTOR = 0.5000
FACTOR = 0.5000
function plotit(charges, location, ckEField, ckEq, DLE, DLV, NLE, NLV, PTS)
figure;
hold on
% Program for plotting the electric field lines
% and equipotential lines due to coplanar point charges
% the plot is to be within the range -5<x,y<5
%
% This is the correct usage:
% function plotit(charges, location,ckEField,ckEq,DLE,DLV,NLE,NLV,PTS)
%
% where,
% charges = a vector containing the charges
% location = a matrix where each row is a charge location
% ckEField = Flag set to 1 plots the Efield lines
% ckEq = Flag set to 1 plots the Equipotential lines
% DLE or DLV = the increment along E & V lines
% NLE = No. of E-Field lines per charge
% NLV = No. of Equipotential lines per charge
% PTS => Plots every PTS point (i.e. if PTS = 5 then plot every 5th point)
% note that constant Q/4*Pie*ErR is set equal to 1.0
% Determine the E-Field Lines
% For convenience, the starting points( XS,YS) are radially distributed about charge locations
Q=charges;
XQ = location(:,1);
YQ = location(:,2);
JJ=1;
NQ = length(charges);
if (ckEField)
for K=1:NQ
for I =1:NLE
THETA = 2*pi*(I-1)/(NLE);
XS=XQ(K)+0.1*cos(THETA);
YS=YQ(K)+0.1*sin(THETA);
XE=XS;
YE=YS;
JJ=JJ+1;
if (~mod(JJ,PTS))
plot(XE, YE, 'k.')
end
while (1)
% Find increment and new point (X,Y)
EX=0;
EY=0;
for J=1:NQ;
R =sqrt((XE-XQ(J))^2 + (YE - YQ(J))^2 );
EX = EX +Q(J)*(XE-XQ(J))/(R^3);
EY = EY +Q(J)*(YE-YQ(J))/(R^3);
end
E = sqrt(EX^2 + EY^2);
% CHECK FOR A SINGULAR POINT
if (E <=0.00005)
break;
end
DX = DLE*EX/E;
DY = DLE*EY/E;
% FOR NEGATIVE CHARGE, NEGATE DX & DY SO THAT INCREMENT IS AWAY FROM THE CHARGE
if (Q(K) < 0)
DX = -DX;
DY = -DY;
end
XE = XE + DX;
YE = YE + DY;
% CHECK WHETHER NEW POINT IS WITHIN THE GIVEN RANGE OR TOO
% CLOSE TO ANY OF THE POINT CHARGES - TO AVOID SINGULAR POINT
if ((abs(XE) >= 5) | (abs(YE) >= 5))
break;
end
if (sum(abs(XE-XQ) < 0.05 & abs(YE-YQ) < 0.05) > 0)
break;
end
JJ=JJ+1;
if (~mod(JJ,PTS))
plot(XE,YE,'k.')
end
end % while loop
end % I =1:NLE
end % K = 1:NQ
end % if
% NEXT, DETERMINE THE EQUIPOTENTIAL LINES FOR CONVENIENCE, THE STARTING POINTS (XS,YS) ARE
% CHOSEN LIKE THOSE FOR THE E-FIELD LINES
if(ckEq)
JJ=1;
DELTA = 0.2;
ANGLE = 45*pi/180;
for K =1:NQ
FACTOR = 0.5
for KK = 1:NLV
XS = XQ(K) + FACTOR*cos(ANGLE);
YS = YQ(K) + FACTOR*sin(ANGLE);
if ( abs(XS) >= 5 | abs(YS) >=5 )
break;
end
DIR = 1;
XV = XS;
YV = YS;
JJ=JJ+1;
if (~mod(JJ,PTS))
plot(XV,YV, 'rs')
end
% FIND INCREMENT AND NEW POINT (XV,YV)
N=1;
while(1)
EX = 0;
EY = 0;
for J = 1:NQ
R = sqrt((XV-XQ(J))^2 + (YV-YQ(J))^2);
EX = EX + Q(J)*(XV-XQ(J))/(R^3);
EY = EY + Q(J)*(YV-YQ(J))/(R^3);
end
E=sqrt(EX^2 + EY^2);
if (E <= 0.00005)
FACTOR = 2*FACTOR;
break;
end;
DX = -(DLV*EX)/E;
DY = (DLV*EY)/E;
XV = XV + DIR*DX;
YV = YV + DIR*DY;
% CHECK IF THE EQUIPOTENTIAL LINE LOOPS BACK TO (X,YS)
R0 = sqrt((XV - XS)^2 + (YV - YS)^2);
if (R0 < DELTA & N < 50)
FACTOR = 2*FACTOR;
break;
end
% CHECK WHETHER NEW POINT IS WITHIN THE GIVEN RANGE IF FOUND OUT OF RANGE, GO BACK TO THE STARTING POINT
% (S,YS) BUT INCREMENT IN THE OPPOSITE DIRECTION
if (abs(XV) > 5 | abs(YV) > 5)
DIR = DIR - 2;
XV = XS;
YV = YS;
end
if (abs(DIR) > 1)
FACTOR = 2*FACTOR;
break;
end
if ( sum( abs(XV-XQ) < 0.005 & abs(YV-YQ) < 0.005) > 0 )
break;
end
end
JJ=JJ+1;
if (~mod(JJ,PTS))
N=N+1;
plot(XV,YV,'rs')
end
end % WHILE loop
end % KK
end % K
end % if
  6 个评论
显示 4更早的评论
dpb
dpb 2024-7-15
for KK = 1:NLV
XS = XQ(K) + FACTOR*cosd(ANGLE);
YS = YQ(K) + FACTOR*sind(ANGLE);
if ( abs(XS) >= 5 | abs(YS) >=5 )
break;
end
end
Nota Bene: The loop over KK NLV times has nothing inside it that is dependent upon KK...there's no point in the loop at all as written.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

David Gabriel
David Gabriel 2024-7-25
OK. This is as good as a day of trial and error could get it with altering the other parameters as mentioned it seems that the equipotentential lines are plotted like spirals eminating from the charge, plus the lines double back on themselves. I will include the altered code, if anyone can tell me why it does not plot as circles.
dg
plotit ( [1 -1], [-1 0; 1 0], 1, 1, 0.01, 0.03, 20, 2, 1)
% plotit ( [1 -1], [-1 0; 1 0], 1, 1, 0.01, 0.03, 20, 1, 1)
function plotit(charges,location,ckEField,ckEq,DLE,DLV,NLE,NLV,PTS)
figure;
hold on
% Program for plotting the electric field lines
% and equipotential lines due to coplanar point charges
% the plot is to be within the range -5<x,y<5
%
% This is the correct usage:
% function plotit(charges, location,ckEField,ckEq,DLE,DLV,NLE,NLV,PTS)
%
% where,
% 1 = charges = a vector containing the charges
% 2 = location = a matrix where each row is a charge location
% 3 = ckEField = Flag set to 1 plots the Efield lines
% 4 = ckEq = Flag set to 1 plots the Equipotential lines
% 5 = DLE = the increment along E lines *** resolution of E lines: <=0.01
% 6 = DLV = the increment along V lines
% 7 = NLE = No. of E-Field lines per charge
% 8 = NLV = No. of Equipotential lines per charge
% 9 = PTS => Plots every PTS point (i.e. if PTS = 5 then plot every 5th point)
% note that constant Q/4*Pie*ErR is set equal to 1.0
% Determine the E-Field Lines
% For convenience, the starting points( XS,YS) are radially distributed about charge locations
Q=charges;
XQ = location(:,1);
YQ = location(:,2);
JJ=1;
NQ = length(charges);
if (ckEField)
for K=1:NQ
for I =1:NLE
THETA = 2*pi*(I-1)/(NLE);
XS=XQ(K)+0.1*cos(THETA);
YS=YQ(K)+0.1*sin(THETA);
XE=XS;
YE=YS;
JJ=JJ+1;
if (~mod(JJ,PTS))
plot(XE, YE, 'k.')
end
while (1)
% Find increment and new point (X,Y)
EX=0;
EY=0;
for J=1:NQ;
R =sqrt((XE-XQ(J))^2 + (YE - YQ(J))^2 );
EX = EX +Q(J)*(XE-XQ(J))/(R^3);
EY = EY +Q(J)*(YE-YQ(J))/(R^3);
end
E = sqrt(EX^2 + EY^2);
% % CHECK FOR A SINGULAR POINT
if (E <=0.0005)
break;
end
DX = DLE*EX/E;
DY = DLE*EY/E;
% FOR NEGATIVE CHARGE, NEGATE DX & DY SO THAT INCREMENT IS AWAY FROM THE CHARGE
if (Q(K) < 0)
DX = -DX;
DY = -DY;
end
XE = XE + DX;
YE = YE + DY;
% CHECK WHETHER NEW POINT IS WITHIN THE GIVEN RANGE OR TOO
% CLOSE TO ANY OF THE POINT CHARGES - TO AVOID SINGULAR POINT
if ((abs(XE) >= 5) | (abs(YE) >= 5))
break;
end
if (sum(abs(XE-XQ) < 0.05 & abs(YE-YQ) < 0.05) > 0) % 0.05 affect appearance of field lines
break;
end
JJ=JJ+1;
if (~mod(JJ,PTS))
plot(XE,YE,'k.')
end
end % while loop
end % I =1:NLE
end % K = 1:NQ
end % if
% NEXT, DETERMINE THE EQUIPOTENTIAL LINES FOR CONVENIENCE, THE STARTING POINTS (XS,YS) ARE
% CHOSEN LIKE THOSE FOR THE E-FIELD LINES
if(ckEq)
JJ=1;
DELTA = 0.01; % resolution of dotted lines
ANGLE = 45;
for K =1:NQ
FACTOR = 0.2 % must be less than 0.001
for KK = 1:NLV
XS = XQ(K) + FACTOR*cosd(ANGLE);
YS = YQ(K) + FACTOR*sind(ANGLE);
if ( abs(XS) >= 5 | abs(YS) >=5 )
break;
end
DIR = 1;
XV = XS;
YV = YS;
JJ=JJ+1;
if ~(mod(JJ,PTS))
plot(XV,YV, 'r.')
end
% FIND INCREMENT AND NEW POINT (XV,YV)
N=1;
while(1)
EX = 0;
EY = 0;
for J = 1:NQ
R = sqrt((XV-XQ(J))^2 + (YV-YQ(J))^2);
EX = EX + Q(J)*(XV-XQ(J))/(R^3);
EY = EY + Q(J)*(YV-YQ(J))/(R^3);
end
E=sqrt(EX^2 + EY^2);
if (E <= 0.0005)
FACTOR = 2*FACTOR;
break;
end
DX = -DLV*EY/E;
DY = DLV*EX/E;
XV = XV + DIR*DX;
YV = YV + DIR*DY;
% CHECK IF THE EQUIPOTENTIAL LINE LOOPS BACK TO (X,YS)
R0 = sqrt((XV - XS)^2 + (YV - YS)^2);
if (R0 < DELTA && N < 50)
FACTOR = 2*FACTOR;
break;
end
% CHECK WHETHER NEW POINT IS WITHIN THE GIVEN RANGE IF FOUND OUT OF RANGE, GO BACK TO THE STARTING POINT (S,YS) BUT INCREMENT IN THE OPPOSITE DIRECTION
if (abs(XV) > 5 || abs(YV) > 5)
DIR = DIR - 2;
XV = XS;
YV = YS;
end
if (abs(DIR) > 1)
FACTOR = 2*FACTOR;
break;
end
if ( sum( abs(XV-XQ) < 0.005 & abs(YV-YQ) < 0.005) > 0 )
break;
end
%%--Moved this code in the while loop--%%
JJ=JJ+1;
if ~(mod(JJ,PTS))
N=N+1;
plot(XV,YV,'r.')
end
end % WHILE loop
end % KK
end % K
end % if ckEq
end % if
Milan Bansal
Milan Bansal 2024-7-25
Ran in:
Hi David Gabriel,
I see that you are trying to plot electrict field lines and equipotential points but you are only able to plot a few "red squares" for the equipotential points due to a possible in the code.
In the attached code you have put the logic of plotting the squares (mentioned below) outside the "while(1)" loop but it seems it should be within the loop.
% JJ = JJ + 1;
% if ~(mod(JJ, PTS))
% N = N + 1;
% plot(XV, YV, 'rs')
% end
Here is the complete code after making this correction:
plotit ( [1 -1], [-1 0; 1 0], 1, 1, 0.01, 1, 20, 20, 1)
function plotit(charges,location,ckEField,ckEq,DLE,DLV,NLE,NLV,PTS)
figure;
hold on
% Program for plotting the electric field lines
% and equipotential lines due to coplanar point charges
% the plot is to be within the range -5<x,y<5
%
% This is the correct usage:
% function plotit(charges, location,ckEField,ckEq,DLE,DLV,NLE,NLV,PTS)
%
% where,
% charges = a vector containing the charges
% location = a matrix where each row is a charge location
% ckEField = Flag set to 1 plots the Efield lines
% ckEq = Flag set to 1 plots the Equipotential lines
% DLE or DLV = the increment along E & V lines
% NLE = No. of E-Field lines per charge
% NLV = No. of Equipotential lines per charge
% PTS => Plots every PTS point (i.e. if PTS = 5 then plot every 5th point)
% note that constant Q/4*Pie*ErR is set equal to 1.0
% Determine the E-Field Lines
% For convenience, the starting points( XS,YS) are radially distributed about charge locations
Q=charges;
XQ = location(:,1);
YQ = location(:,2);
JJ=1;
NQ = length(charges);
if (ckEField)
for K=1:NQ
for I =1:NLE
THETA = 2*pi*(I-1)/(NLE);
XS=XQ(K)+0.1*cos(THETA);
YS=YQ(K)+0.1*sin(THETA);
XE=XS;
YE=YS;
JJ=JJ+1;
if (~mod(JJ,PTS))
plot(XE, YE, 'k.')
end
while (1)
% Find increment and new point (X,Y)
EX=0;
EY=0;
for J=1:NQ
R =sqrt((XE-XQ(J))^2 + (YE - YQ(J))^2 );
EX = EX +Q(J)*(XE-XQ(J))/(R^3);
EY = EY +Q(J)*(YE-YQ(J))/(R^3);
end
E = sqrt(EX^2 + EY^2);
% % CHECK FOR A SINGULAR POINT
if (E <=0.00005)
break;
end
DX = DLE*EX/E;
DY = DLE*EY/E;
% FOR NEGATIVE CHARGE, NEGATE DX & DY SO THAT INCREMENT IS AWAY FROM THE CHARGE
if (Q(K) < 0)
DX = -DX;
DY = -DY;
end
XE = XE + DX;
YE = YE + DY;
% CHECK WHETHER NEW POINT IS WITHIN THE GIVEN RANGE OR TOO
% CLOSE TO ANY OF THE POINT CHARGES - TO AVOID SINGULAR POINT
if ((abs(XE) >= 5) || (abs(YE) >= 5))
break;
end
if (sum(abs(XE-XQ) < 0.05 & abs(YE-YQ) < 0.05) > 0)
break;
end
JJ=JJ+1;
if (~mod(JJ,PTS))
plot(XE,YE,'k.')
end
end % while loop
end % I =1:NLE
end % K = 1:NQ
end % if
% NEXT, DETERMINE THE EQUIPOTENTIAL LINES FOR CONVENIENCE, THE STARTING POINTS (XS,YS) ARE
% CHOSEN LIKE THOSE FOR THE E-FIELD LINES
% Determine the Equipotential Lines
if (ckEq)
JJ = 1;
DELTA = 0.2;
ANGLE = 45; % Added missing semicolon
for K = 1:NQ
FACTOR = 0.5;
for KK = 1:NLV
XS = XQ(K) + FACTOR * cosd(ANGLE);
YS = YQ(K) + FACTOR * sind(ANGLE);
if (abs(XS) >= 5 || abs(YS) >= 5)
break;
end
DIR = 1;
XV = XS;
YV = YS;
JJ = JJ + 1;
if ~(mod(JJ, PTS))
plot(XV, YV, 'rs')
end
% Find increment and new point (XV, YV)
N = 1;
while (1)
EX = 0;
EY = 0;
for J = 1:NQ
R = sqrt((XV - XQ(J))^2 + (YV - YQ(J))^2);
EX = EX + Q(J) * (XV - XQ(J)) / (R^3);
EY = EY + Q(J) * (YV - YQ(J)) / (R^3);
end
E = sqrt(EX^2 + EY^2);
if (E <= 0.00005)
FACTOR = 2 * FACTOR;
break;
end
DX = -(DLV * EY) / E;
DY = (DLV * EX) / E;
XV = XV + DIR * DX;
YV = YV + DIR * DY;
% Check if the equipotential line loops back to (XS, YS)
R0 = sqrt((XV - XS)^2 + (YV - YS)^2);
if (R0 < DELTA && N < 50)
FACTOR = 2 * FACTOR;
break;
end
% Check whether the new point is within the given range
if (abs(XV) > 5 || abs(YV) > 5)
DIR = DIR - 2;
XV = XS;
YV = YS;
end
if (abs(DIR) > 1)
FACTOR = 2 * FACTOR;
break;
end
if (sum(abs(XV - XQ) < 0.005 & abs(YV - YQ) < 0.005) > 0)
break;
end
%%--Moved this code in the while loop--%%
JJ = JJ + 1;
if ~(mod(JJ, PTS))
N = N + 1;
plot(XV, YV, 'rs')
end
end % while loop
end % KK loop
end % K loop
end % if ckEq
end % if
In the output plot as shown above, a lot of red squares are plotted in the middle of the field lines which seems to be the correct location of equipotential points.
Hope this helps!
  1 个评论
David Gabriel
David Gabriel 2024-7-25
Thank you so much. You have been an enormous help as you can see below. I used the square symbol only to see the plot when only few to none were being graphed. When I cange the symbol to '.' to and alter the different paramters which affect the resolution, the lines start to appear. A you can see this is only the beginning. I can now play through trial and error how to get them to show properly. And something is affecting the symmetruy and amongst other things to figure out. I will post when I finally come up an answer unless someone who more knowledgable and can see it without the many trials and errors can do it first. But, we are close, thanks to your careful eye.
Bets Wishes,
David

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Symbolic Math Toolbox 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by