program prog062
c
c     Kepler Motion in the polar coordinate system
c     by using the improved Euler method
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
      t  = ti
      x  = xi
      y  = yi
      vx = vxi
      vy = vyi
c
      theta = atan2(y,x)
      r = sqrt(x*x+y*y)
      u = vx*cos(theta)+vy*sin(theta)
      w = (-vx*sin(theta)+vy*cos(theta))/r
c
         write(*,800) x,y,t
      do i = 1,n
c Trial Propagation
         call force(r,u,w,fu,fw)
         r1     = r + u*dt
         u1     = u + fu*dt
         w1     = w + fw*dt
c
         call force(r1,u1,w1,fu1,fw1)
c
         dr     = ( u  + u1  )*dt/2
         dtheta = ( w  + w1  )*dt/2
         du     = ( fu + fu1 )*dt/2
         dw     = ( fw + fw1 )*dt/2

c Real Propagation
         t      = t + dt
         r      = r + dr
         theta  = theta + dtheta
         u      = u + du
         w      = w + dw
c
         write(*,800) r*cos(theta),r*sin(theta),t
      end do
c
      stop
 800  format(3(1x,f15.7))
      end
c **********************************************************************
      subroutine force(r,u,w,fu,fw)
c **********************************************************************
      implicit real*8 (a-h,o-z)
c
      fu = r*w*w - 1.0d0/r/r
      fw = -2.0d0*u*w/r
      end