c===========================================================
c     dumb: Uses LSODA to integrate equations of motion for 
c     orbiting dumbbell.
c
c     Performs output via 'vsxynt' interface.  Currently
c     configured to write to .sedgat files which can 
c     then be sent to 'ser' or 'jser' visualization servers.
c
c     Output to standard output is 
c
c Column:    1  2  3  4  5 6     7     8    9    10   11
c Quantity:  t x1 y1 x2 y2 theta omega KE_t KE_r PE_g E_tot
c
c     at t=0, dtout, 2 dtout, ... , tmax
c===========================================================
      program        dumb

      implicit       none

      character*2    itoc
      real*8         r8arg
      integer        iargc,         indlnb

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

c-----------------------------------------------------------
c     Command-line arguments
c-----------------------------------------------------------
      real*8         tmax,       dtout

c-----------------------------------------------------------
c     LSODA Variables.
c-----------------------------------------------------------
      external       fcn,        jac

      integer        neq
      parameter    ( neq = 6 )

      real*8         y(neq),     yprime(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         tol
      real*8         default_tol
      parameter    ( default_tol = 1.0d-6 )

c-----------------------------------------------------------
c     Common communication with routine 'fcn' in 'fcn.f'.
c-----------------------------------------------------------
      include       'fcn.inc'
c-----------------------------------------------------------
c     Locals
c-----------------------------------------------------------
      real*8         t,          tout

c-----------------------------------------------------------
c     Parse command line arguments.
c-----------------------------------------------------------
      if( iargc() .ne. 4 ) go to 900

      y(4)  = r8arg(1,r8_never)
      tmax  = r8arg(2,r8_never)
      dtout = r8arg(3,r8_never)
      tol   = r8arg(4,r8_never)
      if( y(4)  .eq. r8_never  .or. tmax .eq. r8_never  .or. 
     &    dtout .eq. r8_never  .or. tol  .eq. r8_never ) 
     &  go to 900
 
c-----------------------------------------------------------
c     Hard-code the remainder of the problem parameters:
c
c     ( x_c(0),  y_c(0)  )   =  ( 1.0 , 0.0 )
c     ( vx_c(0), vy_c(0) )   =  ( 0.0 , vy0 )
c
c     theta(0) = 0
c     omega(0) = 0
c
c     m1/m2    = 0.5
c     d        = 0.1
c-----------------------------------------------------------
      G  = 1.0d0
      MM = 1.0d0

      y(1)  = 1.0d0
      y(2)  = 0.0d0
      y(3)  = 0.0d0

      y(5)  = 0.0d0
      y(6)  = 0.0d0

c     m1bym2 = 0.5d0
c     Energy not conserved for m1bym2 .ne. 1
      m1bym2 = 2.0d0
      mu     = 1.0d0 / (1.0d0 + m1bym2)
      d      = 0.1d0

      m1     = 1.0d0 
      m2     = m1 / m1bym2

      d1     = 1.0d0 / (1.0d0 + m1bym2) * d
      d2     = d - d1

c-----------------------------------------------------------
c     Set LSODA parameters.
c-----------------------------------------------------------
      itol = 1
      rtol = tol
      atol = tol
      itask = 1 
      iopt = 0
      jt = 2

c-----------------------------------------------------------
c     Call the RHS-evaluating routine to initialize the 
c     auxiliary quantities, and output initial values.
c-----------------------------------------------------------
      t = 0.0d0
      call fcn(neq,t,y,yprime)
      write(*,1100) t, x1, y1, x2, y2, th, om, 
     &              ketrans, kerot, pegrav, etot
1100  format(1P,12E24.16,0P)

c-----------------------------------------------------------
c     Do the integration.
c-----------------------------------------------------------
      do while( t .le. tmax )
         istate = 1
         tout   = t + dtout  
         call lsoda(fcn,neq,y,t,tout,
     &              itol,rtol,atol,itask,    
     &              istate,iopt,rwork,lrw,iwork,liw,jac,jt)
         if( istate .lt. 0 ) then
            write(0,*) 'dumb:  Error return ', istate, 
     &                 ' from LSODA '
            write(0,*) 'dumb: Current interval ', 
     &                 t, t + dtout
            stop
         end if
c-----------------------------------------------------------
c     Call the RHS-evaluating routine to compute the 
c     auxiliary quantities, and output them.
c-----------------------------------------------------------
         call fcn(neq,t,y,yprime)
         write(*,1100) t, x1, y1, x2, y2, th, om, 
     &                 ketrans, kerot, pegrav, etot
      end do

      stop

900   continue
         write(0,*) 'usage: dumb <y_0> <tmax> <dtout> <tol>'
         write(0,*) ' '
         write(0,*) 
     &      '       Specify <y_0> = 1.0 for circular orbit'
      stop
      end
