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'
      read (*,*) r0
      write (*,*) ' semiaxes lengths: a b c'
      read (*,*) a,b,c
      write (*,*) ' Nu Nv '
      read (*,*) 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 directed area
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