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
| > |
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;
|
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);
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
|
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; |
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;
|
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;
|
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;
|
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
| > |
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.
