c=======================================================================
c     twsplt:  This program generates a strobe plot of the driven 
c     triangle well system.  This system has a Hamiltonian
c     H = p**2/(2*m) - V0*x - V*x*cos(w*t+theta)    x < 600
c     The system is strobed at the period of the driving field.
c     Because of the presence of the infinite wall the program
c     do2cje from the NAG library is used.  This program integrates
c     the equations of motion until the particle hits the wall.
c     The momentum is then reversed and the equations of motion
c     are again integrated until another collision, etc.
c=======================================================================
      program twsplt

      implicit none

c-----------------------------------------------------------------------
c     n: number of odes
c     iwork: size of work array used by d02ejf
c-----------------------------------------------------------------------
      integer n, iwork
      parameter (n = 2, iwork = 28 + 21*n)

c-----------------------------------------------------------------------
c     tbgn: initial time of integration (overwritten by do2ejf)
c     tend: final time
c     tol: tolerance for d02cje
c-----------------------------------------------------------------------
      real   tbgn, tend, tdt, dt
      parameter (dt = 1.0e-4)
      real   tol, default_tol
      parameter (default_tol = 1.0e-8)

c-----------------------------------------------------------------------
c     i: index used to increment output times
c     ifail: error return from do2cje
c-----------------------------------------------------------------------
      integer ifail

c-----------------------------------------------------------------------
c     work: work array used by d02ejf
c     y: solution vector
c     g: function that monitors collisions with the wall
c-----------------------------------------------------------------------
      real   work(iwork), y(n)
      real   g

c-----------------------------------------------------------------------
c     pmax: maximum value of p used as initial condition
c     xmin: minimum value of x "
c     np: number of initial values of p
c     nx: "                           x
c     niter: number of strobe times (approx)
c-----------------------------------------------------------------------
      real pmax, xmin
      integer np, nx, niter
      integer ip, ix

c-----------------------------------------------------------------------
c     argument parsing variables
c-----------------------------------------------------------------------
      real r8arg, r8_never
      parameter (r8_never = -1.0e-60)
      integer iargc

c-----------------------------------------------------------------------
c     fcn: subroutine that calculates deriviatives
c     pederv: "   "     "  Jacobian
c     out: subroutine that outputs solution at specified times
c-----------------------------------------------------------------------
      external g, d02cje, fcn, out, cjwd02, cjxd02

c-----------------------------------------------------------------------
c     Common communication with subroutines
c-----------------------------------------------------------------------
      include 'fcn.inc'

      real   pi

c-----------------------------------------------------------------------
c     parse command line arguments
c-----------------------------------------------------------------------
      if (iargc() .lt. 2) go to 900

      w = r8arg(1, r8_never)
      v = r8arg(2, r8_never)
      if (w .eq. r8_never .or. v .eq. r8_never) go to 900

      tol = r8arg(3, default_tol)
      if (tol .le. 0.0e0) tol = default_tol

c-----------------------------------------------------------------------
c     Set values of parameters
c-----------------------------------------------------------------------
      v0 = 4.0e-1
      m = 1.0e0
      theta = 0.0e0
      niter = 200
      pi = 4.0e0*atan(1.0e0)
      per = 2.0e0*pi/w
      tend = 2.0e0*niter*per
      l = 6.0e2

c-----------------------------------------------------------------------
c     set iteration parameters
c-----------------------------------------------------------------------
      pmax = 2.1e1
      xmin = 5.4e2
      np = 1
      nx = 40

c-----------------------------------------------------------------------
c     iterate over initial conditions
c-----------------------------------------------------------------------
      do ip = 1, np
         
         do ix = 1, nx

c-----------------------------------------------------------------------
c     Initialize solution at t = 0
c-----------------------------------------------------------------------
            tbgn = 0.0e0
            y(1) = xmin + real(ix-2)*(l-xmin)/real(nx-1)
            y(2) = 0.0e0

            i = 1

 100        continue

c-----------------------------------------------------------------------
c     Call d02ejf to integrate y until collsion with wall
c-----------------------------------------------------------------------
            ifail = 0

            call d02cje(tbgn,tend,n,y,fcn,tol,'D',out,g,work,ifail)

c-----------------------------------------------------------------------
c     Error message if d02cje fails
c-----------------------------------------------------------------------
            if (ifail .gt. 0) then
               write(0,*) 'd02cje returns ifail = ',ifail
               stop
            end if

c-----------------------------------------------------------------------
c     Reverse momentum and return to d02cje unless at the end
c-----------------------------------------------------------------------
            if (tbgn .lt. niter*per) then
               y(2) = -y(2)
               ifail = 0
               i = i - 1
               tdt = tbgn + dt
               call d02cje(tbgn,tdt,n,y,fcn,tol,'D',cjxd02,cjwd02,
     &              work,ifail)
               go to 100
            end if
         end do
      end do

      stop

c-----------------------------------------------------------------------
c     usage message
c-----------------------------------------------------------------------
 900  continue
      
      write(0,*) 'usage: twsplt <w> <V> [<tol>]'
      write(0,*)
      write(0,*) 'w: frequency of driving field'
      write(0,*) 'V: strength of driving field'
      write(0,*) 'tol: optional tolerance for d02cje'

      stop

      end