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)

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

[dpb]

Format the code first, then show/tell us explicitly what is the problem.

[David Gabriel]

The following plot is the output of the function below. The electric field lines are perfect. The little red squares are the plotted electric potential. In the code, I had to make them red squares or I would not see them at all. The code plots each point as a dot and the increment is small enough so that the field lines plot as a line. But, the electric porential are not plotting. Just a few points, which are the squares.

Thank you for your patience while I navigate this.

Best Wishes,

David.

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

if(ckEq)
    JJ=1;
    DELTA = 0.2;
    ANGLE = 45
    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 (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

[dpb]

plotit ( [1 -1], [-1 0; 1 0], 0, 1, 0.01, 0.001, 20, 200, 1)

   200

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)
    disp(NLV)
    JJ=1;
    DELTA = 0.2;
    ANGLE = 45;
    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 (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

The difference is in the logic used to compute the points -- I don't know precisely what it is being calculated without more looking into than have time for, but in the first section

     for I =1:NLE       
            THETA = 2*pi*(I-1)/(NLE);
            XS=XQ(K)+0.1*cos(THETA);
            YS=YQ(K)+0.1*sin(THETA);
            ...

the x,y points are dependent upon the number NLE items and the theta is generated for those number of points and so XS, YS is dependent upon those theta.

In the other section,

    ANGLE = 45;
    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
        end
     ...        

the ANGLE is constant and given there are only two charges NLV == 2 and so only two sets of four points each are within the bounds regardless of what is given for NLV or the DLV input parameters.

You'll need to consider just what it is that is being calculated and how to refine the position calculation that makes sense, but the basic logic appears to be faulty in that section.

I set the NLV and DLV parameters in the above and turned off the first to 200 for NLV and you can see it just recalculated the same four points for each charge 200 times given there is no change in the angle being used based on that number of iterations.

[David Gabriel]

Thank you very much for your time. I can take a look into that section more deeply.

Best Wishes,

David

[dpb]

        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.

[David Gabriel]

Thank you for noting that.

dg

Answer 1

[David Gabriel]

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

Answer 2

[Milan Bansal]

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!

[David Gabriel]

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