Application Center - Maplesoft
  Submit your work
Application Center Applications 3D Equipotential and electric field lines due to point charges Preview

App Preview:

3D Equipotential and electric field lines due to point charges

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application

 

Image 

3D Equipotential and electric field lines due to point charges 

Takao Takeuchi, Ph.D.
Department of Mathematics and Physics
State University of New York at Alfred
Alfred, New York 14802 USA
takeuct@alfredstate.edu http://web.alfredstate.edu/takeuchi/index.htm      ? 2006 Takao Takeuchi 

Introduction 

This worksheet demonstrates the use of Maple for calculating and displaying three dimensional equipotential surfaces and electric field lines due to point charges of your choice. It takes quite a long time to plot electric field lines if your computer's cpu is slow. Please modify at your own risk. Try with only two charges first to learn how to use it.Be sure to input reasonable coordinates and charges. 

Ver. 9.5.1 

Initialization 

> restart:
 

> with(plots):with(plottools):
 

Input1 

The procedure below obtains plotting area size and calls subroutines InputXYZQ and CheckInput. 

> Input1:=proc(xxmin,xxmax,yymin,yymax,zzmin,zzmax)
  global ncontours,xmin,xmax,ymin,ymax,zmin,zmax;
  if (xxmin>=xxmax) then
  print ("Error! xxmin should be smaller that xxmax");
  break;
  end if;
  if (yymin>=yymax) then
     print ("Error! yymin should be smaller that yymax");
     break;
  end if;
  if (zzmin>=zzmax) then
     print ("Error! zzmin should be smaller that zmax");
     break;
  end if;
  xmin:=xxmin;
  xmax:=xxmax;
  ymin:=yymin;
  ymax:=yymax;
  zmin:=zzmin;
  zmax:=zzmax;
  
  InputXYZQ(xxmin,xxmax,yymin,yymax,zzmin,zzmax);
  CheckInput();
end proc:
 

InputXYZQ 

The procedure below obtains coordinates and charge of each charge from the user. 

> InputXYZQ:=proc(xxmin,xxmax,yymin,yymax,zzmin,zzmax)
global Nocharge,xx,yy,zz,Q,Chargeall;
local i,j,charge;
Nocharge:=-10;
while Nocharge<=0 do
  print("Be sure that the number of charges is positive. ");
  Nocharge:=readstat("How many charges? ");
end do;

for i from 1 to Nocharge do
  printf("For charge %d\n",i);
  xx[i]:=readstat(" XX of charge: ");
  if (xx[i]<xxmin) or (xx[i]>xxmax) then
     printf("Warning! The input value XX not between xxmin= %d and xxmax= %d.\n", xxmin, xxmax);
  end if;
  yy[i]:=readstat(" YY of charge: ");
  if (yy[i]<yymin) or (yy[i]>yymax) then
     printf("Warning! The input value YY not between yymin= %d and yymax= %d.\n", yymin, yymax);
  end if;
  zz[i]:=readstat(" ZZ of charge: ");
  if (zz[i]<zzmin) or (zz[i]>zzmax) then
     printf("Warning! The input value ZZ not between zzmin= %d and zzmax= %d.\n", zzmin, zzmax);
  end if;

  Q[i]:=readstat(" Q of charge: ");
  if (Q[i]=0.0) then
     printf("Warning! The charge is zero.");
  end if;
end do;
    
 

> end proc:
 

DrawCharge 

The procedure below creates a plot structure of the charges. 

> DrawCharge:=proc(xmin,xmax,ymin,ymax,zmin,zmax)
  global Nocharge,xx,zz,Q,Chargeall;
  local i,charge,Chargesizex,Chargesizey,Chargesizez,Chargesize;
 
  Chargesizex:=abs(xmax-xmin)*0.5/32;
  Chargesizey:=abs(ymax-ymin)*0.5/32;
  Chargesizez:=abs(zmax-zmin)*0.5/32;
  Chargesize:=min(Chargesizex,Chargesizey,Chargesizez);
 

    for i from 1 to Nocharge do
       if Q[i]>0.0 then  
          charge[i]:=(sphere([xx[i],yy[i],zz[i]],Chargesize,color=blue));
       elif Q[i]<0.0 then  
          charge[i]:=(sphere([xx[i],yy[i],zz[i]],Chargesize,color=red));
       else
          charge[i]:=(sphere([xx[i],yy[i],zz[i]],Chargesize,color=yellow));
       end if;
    end do;
    Chargeall:=seq(charge[i],i=1..Nocharge); 

> end proc:
 

CheckInput 

The procedure below checks for duplicate charges. 

> CheckInput:=proc()
  
global Nocharge,xx,yy,zz,Q;
local i,j;
    
for i from 1 to Nocharge do
 for j from 1 to i-1 do
      if xx[i]=xx[j] and yy[i]=yy[j] and zz[i]=zz[j] and Q[i]=Q[j] then
         printf("Warning! Charge %d identical to charge %d\n",i,j);
         printf("Input data for charge %d again!\n",i);
            xx[i]:=readstat(" XX of charge: ");
            yy[i]:=readstat(" YY of charge: ");

            zz[i]:=readstat(" ZZ of charge: ");

            Q[i]:=readstat(" Q of charge: ");
         CheckInput();
      end if;
 end do;     # i
end do;       # j
 

> end proc:
 

CalculateV 

The procedure below computes potential V at (x,y,z) due to the charges. 

> CalculateV:=proc()
  global V,Nocharge,xx,yy,zz,Q;
  local r,i;

  V:=0;
  for i from 1 to Nocharge do
     r:=sqrt((x-xx[i])^2+(y-yy[i])^2+(z-zz[i])^2);
     V:=V+Q[i]/r;
  end do;
 

> end proc:
 

MinMaxV 

The procedure below computes maxmum and  minimum potential due to the charges. 

> MinMaxV:=proc(ncontours,xmin,xmax,ymin,ymax,zmin,zmax)
  global Nocharge,xx,yy,zz,Q,ContourValues;
  local i,ii,jj,kk,x,y,z,r,NgridX,NgridY,NgridZ,DeltaX,DeltaY,DeltaZ,V,Vfac,Vmin,Vmax;
  Vmax:=10^(-30);Vmin:=10^30;
  V:=0;
  NgridX:=10;NgridY:=10;NgridZ:=10;
  DeltaX:=(xmax-xmin)/(NgridX-1);
  DeltaY:=(ymax-ymin)/(NgridY-1);
  DeltaZ:=(zmax-zmin)/(NgridZ-1);
   
  for ii from 1 to NgridX do
     for jj from 1 to NgridY do
        for kk from 1 to NgridZ do
           x:=xmin+DeltaX*(ii-1);
           y:=ymin+DeltaY*(jj-1);
           z:=zmin+DeltaZ*(kk-1);
           V:=0;
           for i from 1 to Nocharge do
              r:=sqrt((x-xx[i])^2+(y-yy[i])^2+(z-zz[i])^2);
              if evalf(r)>=0.01 then V:=V+Q[i]/r end if;
           end do;
           if evalf(V)<=Vmin then Vmin:=evalf(V) end if;
           if evalf(V)>=Vmax then Vmax:=evalf(V) end if;
        end do;
     end do;
   end do;

   printf(" Vmin = %E, Vmax = %E\n ",Vmin, Vmax);
   
  Vfac:=0.85;
  if Vmin>0 and Vmax>0 then
     Vmin:=Vmin/Vfac;
  elif Vmin<0 and Vmax<1.0e-8 then
     Vmax:=Vmin*(1-Vfac);
  end if;

#  printf(" Vmin = %E, Vmax = %E\n ",Vmin, Vmax);
  for i from 1 to ncontours do
     ContourValues[i]:=Vmin+(i-1)*(Vmax-Vmin)/(ncontours-1);
  end do;
 

> end proc:
 

CalculateE 

> CalculateE:=proc()
  global Ex,Ey,Ez,Nocharge,xx,yy,zz,Q;
  local r,i;

  Ex:=0.0;Ey:=0.0;Ez:=0.0;
  for i from 1 to Nocharge do
     r:=sqrt((x-xx[i])^2+(y-yy[i])^2+(z-zz[i])^2);
     if r<>0.0 then  
        Ex:=Ex+Q[i]*(x-xx[i])/r^3;
        Ey:=Ey+Q[i]*(y-yy[i])/r^3;
        Ez:=Ez+Q[i]*(z-zz[i])/r^3;
     end if;
  end do;
 

> end proc:
 

Display3dContoursOfV 

The procedure below displays equipotential surfaces or electric field lines. 

> Display3dContoursOfV:=proc(potflag::nonnegint)
 
local i,c,pcurve,Ecurve,grad3d,graphtitle,ncontours;

if potflag=0 then
  graphtitle:="Electric field lines around all charges";
elif potflag=2 then
  graphtitle:="Electric field lines around one charge";
#elif potflag=3 then
#   graphtitle:="Efieldplot3d";
#elif potflag=4 then
#   graphtitle:="Gradplot3d";
else
  graphtitle:="Equipotential surfaces";
end if;  

setoptions3d(labels=["x","y","z"],view=[xmin..xmax,ymin..ymax,zmin..zmax],style=patchnogrid,axes=boxed,orientation=[-55,53],title=graphtitle);

DrawCharge(xmin,xmax,ymin,ymax,zmin,zmax);
if potflag=1 then
  ncontours:=-1;
  while ncontours<=0 do
     print("Be sure that the number of contours is larger than 1. ");
     ncontours:=readstat("Number of equipotential surfaces(>1): ");
  end do;
  if ncontours=1 then ncontours:=2;end if;
  MinMaxV(ncontours,xmin,xmax,ymin,ymax,zmin,zmax);
  CalculateV();
end if;

if potflag=0 or potflag=2 then
  PlotEfieldLines(potflag,xmin,xmax,ymin,ymax,zmin,zmax);

#elif potflag=3 then
#   CalculateE();
   #Ecurve:=fieldplot3d([Ex,Ey,Ez],x=xmin..xmax,y=ymin..ymax,z=zmin..zmax,grid=[10,10,10],orientatio#n=[-120,60],scaling=constrained,axes=frame,labels=["x","y","z"],arrows=LINE);

#   display(Chargeall,Ecurve);

#elif potflag=4 then
#   grad3d:=gradplot3d(V,x=xmin..xmax,y=ymin..ymax,z=zmin..zmax);
#   display(Chargeall,grad3d);

else
  for i to ncontours do
     c[i]:=V=ContourValues[i];
  end do;
  c[ncontours+1]:=V=0.0;
pcurve:=(implicitplot3d({seq(c[k],k=1..ncontours+1)},x=xmin..xmax,y=ymin..ymax,z=zmin..zmax,grid=[15,15,15],orientation=[-120,60],scaling=constrained,axes=frame,labels=["x","y","z"],transparency=0.0,shading=ZHUE));

  display(Chargeall,pcurve);

end if;
 

> end proc:
 

PlotEfieldLines 

The procedure below generates electric field lines, plot structures and display them.. 

> PlotEfieldLines:=proc(potflag,xmin,xmax,ymin,ymax,zmin,zmax)
global Nocharge,xx,yy,zz,Q,EfieldLineSegment,Chargeall;
local i,j,k,nn,NoFieldLines,s,sx,sy,sz,WhichCharge,EfieldOneLine,EfieldLineAllOneCharge,EfieldLineAll,NNocharge,x0,y0,z0;

sx:=abs(xmax-xmin)*0.1/32;
sy:=abs(ymax-ymin)*0.1/32;
sz:=abs(zmax-zmin)*0.1/32;
s:=min(sx,sy,sz)*2;

NoFieldLines:=1;

if potflag=2 then
  WhichCharge:=0;
  while (WhichCharge<=0 or WhichCharge>Nocharge) do
       WhichCharge:=readstat("Around which charge to draw E-field lines? ");
  end do;
end if;

if potflag=2 then
  NNocharge:=1;
elif potflag=0 then
  NNocharge:=Nocharge;
end if;

for nn to NNocharge do
  if potflag=0 then
     WhichCharge:=nn;
  end if;

  NoFieldLines:=1;

for i in [-1,0,1] do
  x0:=xx[WhichCharge]+i*sx*5;
  for j in [-1,0,1] do
      y0:=yy[WhichCharge]+j*sy*5;
      for k in [-1,0,1] do
        z0:=zz[WhichCharge]+k*sz*5;
#        if not((i=0 and j=0) or (i=0 and k=0) or (j=0 and k=0)) then
          s:=abs(s);
          GetOneEfieldLineSegment(1,x0,y0,z0,s,xmin,xmax,ymin,ymax,zmin,zmax);
  
          s:=-s;
          GetOneEfieldLineSegment(2,x0,y0,z0,s,xmin,xmax,ymin,ymax,zmin,zmax);
  
          EfieldOneLine[NoFieldLines]:=seq(EfieldLineSegment[j],j=1..2);
          NoFieldLines:=NoFieldLines+1
#         end if;
     end do;   # k
  end do;      # j
end do;         # i


EfieldLineAllOneCharge[nn]:=seq(EfieldOneLine[i],i=1..NoFieldLines-1);

end do;         # nn

EfieldLineAll:=seq(EfieldLineAllOneCharge[i],i=1..NNocharge);

display(Chargeall,EfieldLineAll);   

end proc:
 

GetOneEfieldLineSegment 

The procedure below creates one electric field line. 

> GetOneEfieldLineSegment:=proc(segno,x0,y0,z0,s,xmin,xmax,ymin,ymax,zmin,zmax)
global Nocharge,xx,yy,zz,Q,EfieldLineSegment;
local n,i,j,k,smallinc,H1,H2,H3,H1N,H2N,H3N,Rx,Ry,Rz,R,R3,Ex,Ey,Ez,E,x1,y1,z1,xN,yN,zN,LineSegment,Ntime;

Ntime:=1;
smallinc:=0.00001;
LineSegment:=0;

x1:=x0;y1:=y0;z1:=z0;
xN:=x0:yN:=y0;zN:=z0;
H1:=0;H2:=0;H3:=0;
H1N:=0;H2N:=0;H3N:=0;

for k do

  Ex:=0;Ey:=0;Ez:=0;

  for i from 1 to Nocharge do
     Rx:=x1+H1/2-xx[i];
     Ry:=y1+H2/2-yy[i];
     Rz:=z1+H3/2-zz[i];
     R:=sqrt(Rx^2+Ry^2+Rz^2);
     R3:=R^3;
 

       if R<>0.0 then  
          Ex:=Ex+Q[i]*Rx/R3;
          Ey:=Ey+Q[i]*Ry/R3;
          Ez:=Ez+Q[i]*Rz/R3;
       end if;
     end do;

     E:=sqrt(Ex^2+Ey^2+Ez^2);
     if E<>0.0 then
        H1:=s*Ex/E;
        H2:=s*Ey/E;
        H3:=s*Ez/E;
     end if;
     x1:=x1+H1;
     y1:=y1+H2;
     z1:=z1+H3;

     if ((abs(H1+H1N)<smallinc and abs(H2+H2N)<smallinc and abs(H3+H3N)<smallinc)) or Ntime=600  then
       EfieldLineSegment[segno]:=seq(LineSegment[n],n=1..Ntime-1);
       break;
     end if;
  
     for j from 1 to Nocharge do
        if ((abs(x1-xx[j])+abs(y1-yy[j])+abs(z1-zz[j]))<abs(0.9*s)) then
           EfieldLineSegment[segno]:=seq(LineSegment[n],n=1..Ntime-1);
           break;
        end if;
     end do;

     if (x1>=xmax or x1<=xmin or y1>=ymax or y1<=ymin or z1>=zmax or z1<=zmin) then
        EfieldLineSegment[segno]:=seq(LineSegment[n],n=1..Ntime-1);
        break;
     end if;

     LineSegment[Ntime]:=line([xN,yN,zN],[x1,y1,z1]);
     Ntime:=Ntime+1;
#     printf ("k= %d and Ntime= %d", k, Ntime);
     xN:=x1:yN:=y1;zN:=z1;
     H1N:=H1;H2N:=H2;H3N:=H3


  end do;   # k loop
#     printf ("k final= %d and Ntime= %d", k, Ntime);print (xmin);
end proc:
 

Main Procedure 

Important! You have to run both 

1. Input1  

2. Display3dContoursOfV 

 

Input1(xmin,xmax,ymin,ymax,zmin,zmax)
 

xmin, xmax : Window size in x direction   Use number like -1.5, 1.5
ymin, ymax : Window size in y direction   Use number like -2, 2 

zmin, zmax : Window size in z direction   Use number like -2, 2

Input charges
  Number of charges: number of charges you want to consider. Use 1, 2, 3, etc.
  XX(i): x coordinate of the i-th charge    Use numbers like 1, -1  etc.
  YY(i): y coordinate of the i-th charge    Use numbers like 1, -1 etc.
  ZZ(i): z coordinate of the i-th charge    Use numbers like 1, -1, etc.
  Q(i): charge of i-th charge            Use numbers like 1, -2, 3, etc. 

 

Important! Be sure to type a semicolon ";" after typing in an input number in ver. 9.5.1. But in ver.10, type a number and press return.
 

Diplay3dContoursOfV(potflag)
potflag=0 : Electric field lines will be plotted around all charges
           1 : Equipotential surfaces
           2 : Electric field lines will be plotted around one charge of your choice
***********************************************************************  
Reference: J.R. Merrill, "Using Computers in Physics", p.45, 1976, Houghton Mifflin Comp., Boston
*********************************************************************** 

Example 

> Input1(-1.5,1.5,-1.5,1.5,-1,0.1);
 

Be sure that the number of charges is positive.  

For charge 1 

For charge 2
 

> Display3dContoursOfV(1);
 

Be sure that the number of contours is larger than 1.  

Vmin = -2.220633E+00, Vmax = 1.000000E-30
 

Plot
 

>
 

Example output 

Image 

Image 

Image 

Image 

Legal Notice: The copyright for this application is owned by the author(s). Neither Maplesoft nor the author are responsible for any errors contained within and are not liable for any damages resulting from the use of this material. This application is intended for non-commercial, non-profit use only. Contact the author for permission if you wish to use this application in for-profit activities.
 

Image