Looking for bug in a graphics program for plotting dipole fields
14 次查看(过去 30 天)
显示 更早的评论
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:
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)
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 个评论
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 个)
Milan Bansal
2024-7-25
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!
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Number Theory 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!