 Application Center - Maplesoft

# 3D Equipotential and electric field lines due to point charges

You can switch back to the summary page by clicking here. 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]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]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]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); For charge 1

 For charge 2

 > Display3dContoursOfV(1); Vmin = -2.220633E+00, Vmax = 1.000000E-30 >

Example output    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. 