Looking for bug in a graphics program for plotting dipole fields
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)
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 % ifGreetings,
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)
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 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)
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 graphics and visualization pl MATLAB Answers — New Questions