C set constants
      
      dtime=0.001
      q=0.9

      print *, 'Input initial x,y,z'
      read *, x,y,z
      print *, 'Input initial vx,vy,vz'
      read *, vx,vy,vz

      open (unit=10, file='orbit.dat')
      
      time=0

      do while (time.LT.50) 

C        er2 is the square of the ellipsoidal radius
C        acc=-grad(phi)

         er2=x*x+y*y+z*z/(q*q)
         phi=0.5*log(er2)
         ax0=-x/er2
         ay0=-y/er2
         az0=-z/(q*q)/er2

         vxt=vx+ax0*dtime
         vyt=vy+ay0*dtime
         vzt=vz+az0*dtime

         xt=x+vxt*dtime+0.5*ax0*dtime*dtime
         yt=y+vyt*dtime+0.5*ay0*dtime*dtime
         zt=z+vzt*dtime+0.5*az0*dtime*dtime

         ert2=xt*xt+yt*yt+zt*zt/q*q
         axt=-xt/ert2
         ayt=-yt/ert2
         azt=-zt/(q*q)/ert2

         ax=0.5*(ax0+axt)
         ay=0.5*(ay0+ayt)
         az=0.5*(az0+azt)

         x=x+vx*dtime+0.5*ax*dtime*dtime
         y=y+vy*dtime+0.5*ay*dtime*dtime
         z=z+vz*dtime+0.5*az*dtime*dtime

         vx=vx+ax*dtime
         vy=vy+ay*dtime
         vz=vz+az*dtime

         time=time+dtime

         write (10,66) time,x,y,z,vx,vy,vz,phi
 66      format (f10.4,2x,7(e10.4,2x))

      enddo

      close (10)

      end