c===========================================================
c     twobody:  Integrates restricted gravitational 2-body 
c     problem using LSODA.
c
c     usage: twobody <x0> <y0> <vx0> <vy0> <tmax> <dt> [<tol>]
c
c     Output to standard output 
c
c        0.0    x(0.0)   y(0.0)   dEtot(0.0)  dJtot(0.0)
c        dt     x(dt)    y(dt)    dEtot(dt)   dJtot(dt)  
c       2*dt    x(2*dt)  y(2*dt)  dEtot(2*dt) dJtot(2*dt)
c                             .                           
c                             .                           
c                             .                           
c       tmax    x(tmax)  y(tmax)  dEtot(tmax) dJtot(tmax)
c===========================================================
      program        twobody

      implicit       none

      integer        iargc,      i4arg
      real*8         r8arg

      real*8         r8_never
      parameter    ( r8_never = -1.0d-60 )

c-----------------------------------------------------------
c     Command line arguments (initial position and velocity 
c     components will be read directly into y() array).
c-----------------------------------------------------------
      real*8         tmax,       dt,         tol

c-----------------------------------------------------------
c     LSODA Variables.
c-----------------------------------------------------------
      integer        neq
      parameter    ( neq = 4)

      external       fcn,        jac

      real*8         y(neq)
      real*8         tbgn,       tend
      integer        itol
      real*8         rtol,       atol
      integer        itask,      istate,     iopt
      integer        lrw

      parameter    ( lrw = 22 + neq * 16 )
      real*8         rwork(lrw)

      integer        liw
      parameter    ( liw = 20 + neq )
      integer        iwork(liw)
      integer        jt

      real*8         default_tol
      parameter    ( default_tol = 1.0d-6 )

c-----------------------------------------------------------
c     Locals 
c
c     Etot:   Instantaneous total mechanical energy 
c     Jtot:   Instantaneous total angular momentum
c     Etot0:  Initial total mechanical energy 
c     Jtot0:  Initial total angular momentum
c-----------------------------------------------------------
      real*8         t,           ts,         tf
      integer        ieq
      real*8         Etot,        Jtot,
     &               Etot0,       Jtot0

c-----------------------------------------------------------
c     Common communication with routine 'fcn' in 'fcn.f' ...
c-----------------------------------------------------------
      include       'fcn.inc'

c-----------------------------------------------------------
c     Initialize parameters defined  in common block ...
c-----------------------------------------------------------
      G = 1.0d0
      M = 1.0d0

c-----------------------------------------------------------
c     Parse command line arguments (initial values) ...
c-----------------------------------------------------------
      if( iargc() .lt. 6 ) go to 900

      do ieq = 1 , 4
         y(ieq) = r8arg(ieq,r8_never)
         if( y(ieq) .eq. r8_never ) go to 900
      end do
      tmax = r8arg(5,r8_never)
      if( tmax .eq. r8_never ) go to 900
      dt   = r8arg(6,r8_never)
      if( dt   .eq. r8_never ) go to 900
      tol  = r8arg(7,default_tol)

c-----------------------------------------------------------
c     Set LSODA parameters ...
c-----------------------------------------------------------
      itol = 1
      rtol = tol
      atol = tol
      itask = 1 
      iopt = 0
      jt = 2
      
c-----------------------------------------------------------
c     Compute initial energy, angular momentum, then output 
c     initial time, particle coordinates, 
c     Etot - Etot0 and Jtot - Jtot0
c
c     Note use of format statement to ensure that all five
c     numbers are output on a single line, 'write(*,*)'
c     will break lines, inhibiting further processing with 
c     Unix utilities.  The format statement is good for 
c     up to 10 numbers per line.
c-----------------------------------------------------------
      t = 0.0d0
      call calc_ej(y,Etot0,Jtot0)
      call calc_ej(y,Etot,Jtot)

      write(*,1000) t, y(1), y(2),
     &              Etot - Etot0, Jtot - Jtot0
1000  format(1P, 10 E25.16, 0P)

c-----------------------------------------------------------
c     Do the integration ...
c-----------------------------------------------------------
      istate = 1
      do while( t .le. tmax ) 
         ts = t
         tf = t + dt
c-----------------------------------------------------------
c        Integrate EOM from t=ts to t=tf ...
c-----------------------------------------------------------
         call lsoda(fcn,neq,y,ts,tf,
     &              itol,rtol,atol,itask,    
     &              istate,iopt,rwork,lrw,iwork,liw,jac,jt)

c-----------------------------------------------------------
c        Check return code; bail-out with an error message
c        if routine was not successful ...
c-----------------------------------------------------------
         if( istate .lt. 0 ) then
            write(0,*) 'twobody:  Error return ', istate, 
     &                 ' from LSODA '
            write(0,*) 'twobody: Current interval ', t, t + dt
            stop
         end if

         t = t + dt
c-----------------------------------------------------------
c        Compute new energy and angular momentum, output
c        as previously (i.e. use the same format statement)
c-----------------------------------------------------------
         call calc_ej(y,Etot,Jtot)
         write(*,1000) t, y(1), y(2),
     &                 Etot - Etot0, Jtot - Jtot0
      end do

      stop

900   continue
         write(0,*) 'usage: twobody <x0> <y0> <vx0> <vy0> '//
     &              '<tmax> <dt> [<tol>]'
      stop
      end
