program ex061
c
c     Kepler Motion in the x and y coordinate system
c     by using the improved Euler method
c     subroutines:
c        onestep(x,y,vx,vy,dt)
c        force(x,y,fx,fy)
c
      implicit real*8 (a-h,o-z)
c
c
c Default Parameters
      n   = 2000
      ti  = 0.0d0
      tf  = 50.0d0
      dt  = (tf-ti)/n
      xi  = 4.0d0
      yi  = 0.0d0
      vxi = 0.0d0
      vyi = 0.4d0
c
      write(*,810) '# Vy = ?'
      read(*,*) vyi
810   format(a)
c
      t  = ti
      x  = xi
      y  = yi
      vx = vxi
      vy = vyi
c
         call engang(x,y,vx,vy,eng,ang)
         write(*,800) x,y,t,eng,ang
      do i = 1,n
         call onestep(x,y,vx,vy,dt)
         t = t + dt
         call engang(x,y,vx,vy,eng,ang)
         write(*,800) x,y,t,eng,ang
      end do
c
      stop
 800  format(5(1x,f15.7))
      end
c **********************************************************************
      subroutine onestep(x,y,vx,vy,dt)
c **********************************************************************
      implicit real*8 (a-h,o-z)
c
c Force at x,y
      call force(x,y,fx,fy)

c Trial Propagation and Force at x,y
      x1  = x + vx*dt
      y1  = y + vy*dt
      vx1 = vx + fx*dt
      vy1 = vy + fy*dt
      call force(x1,y1,fx1,fy1)
c
      dx =  ( vx + vx1 )*dt/2   ! Improved Euler
      dy =  ( vy + vy1 )*dt/2   ! Improved Euler
      dvx = ( fx + fx1 )*dt/2   ! Improved Euler
      dvy = ( fy + fy1 )*dt/2   ! Improved Euler

c Real Propagation in Improved Euler method
      x = x + dx
      y = y + dy
      vx = vx + dvx
      vy = vy + dvy
c
      return
      end
c **********************************************************************
      subroutine force(x,y,fx,fy)
c **********************************************************************
      implicit real*8 (a-h,o-z)
c
      r = sqrt(x*x+y*y)
      fx = -x/r/r/r
      fy = -y/r/r/r
      return
      end
c **********************************************************************
      subroutine engang(x,y,vx,vy,eng,ang)
c **********************************************************************
      implicit real*8 (a-h,o-z)
      eng=(vx*vx+vy*vy)/2-1.0d0/sqrt(x**2+y**2)
      ang=x*vy-y*vx
      end