Up: Introduction to some Basic Previous: 3 Conclusion

# 4 Example Code

Here is a program that implements the calculation of d over an ellipsoid for an electric field from a unit point charge.
      program surf
implicit double precision (a-h,o-z)
parameter (zero=0.d0,one=1.d0,two=2.d0,three=3.d0,four=4.d0)
parameter (five=5.d0,six=6.d0,seven=7.d0,eight=8.d0,anine=9.d0)
parameter (ten=10.d0,tenth=.1d0,half=.5d0,third=1.d0/3.d0)
dimension r0(3),s00(3),sp0(3),sm0(3),s0p(3),s0m(3),da(3)
dimension ds1(3),ds2(3),dx(3),efield(3)
write (*,*) ' x y z of unit charge'
write (*,*) ' semiaxes lengths: a b c'
write (*,*) ' Nu Nv '
pi=four*atan(one)
hu=pi/nu
hv=two*pi/nv
gauss=zero
c
c here I recalculate all the trig functions -- store them in  a table
c if you want to make this faster
c
c loop over u
c
do 10 i=1,nu
u=(i-half)*hu
su=sin(u)
cu=cos(u)
sup=sin(u+half*hu)
cup=cos(u+half*hu)
sum=sin(u-half*hu)
cum=cos(u-half*hu)
c
c loop over v
c
do 20 j=1,nv
v=(j-half)*hv
sv=sin(v)
cv=cos(v)
svp=sin(v+half*hv)
cvp=cos(v+half*hv)
svm=sin(v-half*hv)
cvm=cos(v-half*hv)
c
c calculate the central and 4 edge surface points
c
s00(1)=a*su*cv
s00(2)=b*su*sv
s00(3)=c*cu
sp0(1)=a*sup*cv
sp0(2)=b*sup*sv
sp0(3)=c*cup
sm0(1)=a*sum*cv
sm0(2)=b*sum*sv
sm0(3)=c*cum
s0p(1)=a*su*cvp
s0p(2)=b*su*svp
s0p(3)=c*cu
s0m(1)=a*su*cvm
s0m(2)=b*su*svm
s0m(3)=c*cu
c
c evaluate the approximate elemental area vector
c
ds1(1)=sp0(1)-sm0(1)
ds1(2)=sp0(2)-sm0(2)
ds1(3)=sp0(3)-sm0(3)
ds2(1)=s0p(1)-s0m(1)
ds2(2)=s0p(2)-s0m(2)
ds2(3)=s0p(3)-s0m(3)
da(1)=ds1(2)*ds2(3)-ds1(3)*ds2(2)
da(2)=ds1(3)*ds2(1)-ds1(1)*ds2(3)
da(3)=ds1(1)*ds2(2)-ds1(2)*ds2(1)
c
c evaluate e field of point charge
c
dx(1)=s00(1)-r0(1)
dx(2)=s00(2)-r0(2)
dx(3)=s00(3)-r0(3)
r=sqrt(dx(1)**2+dx(2)**2+dx(3)**2)
r3i=one/r**3
efield(1)=dx(1)*r3i
efield(2)=dx(2)*r3i
efield(3)=dx(3)*r3i
c
c calculate normal e field times elemental area
c
dot=efield(1)*da(1)+efield(2)*da(2)+efield(3)*da(3)
gauss=gauss+dot
20 continue
10 continue
write (*,*) 'integral over ellipsoid/4 pi = ',gauss/(four*pi)
end