c=======================================================================
c     History: sode.f (Matt Choptuik)
c
c     accel:  This program calulates the expectation value of the 
c     acceleration for a quantum system of the form:
c
c     H = H0 - V*x*cos(w*t + theta)
c
c     where H0 is a time-independent, integrable Hamiltonian.  The
c     calculation proceeds by computing decomposing the time-dep.
c     wavefunction into a basis of eigenstates of H0 and computing
c     the coefficients of this decomposition as a function of time.
c     
c     For a particular system the subroutines 'dipole' and 'getE'
c     must be set up to calculate or read the proper dipole matrix
c     elements and energy eigenvalues for that particular H0.  The
c     common block in 'fcn.inc' will also have to be modified to 
c     contain the correct parameters for the specific H0.
c
c     The output of this program is given as a list of time values
c     and the corresponding value of the (real part) of the expectation
c     value of the acceleration written to standard out.  The imaginary 
c     part of the accel. (should be 0) and the total probability 
c     (should be 1) are written to standard error.
c
c     A Fourier transform of the output of this program gives the 
c     radiation spectrum of the system.  The maximum harmonic that 
c     can be seen will be determined by the command-line 
c     argument, hmax, used with accel.
c=======================================================================
      program accel

      implicit none

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

c-----------------------------------------------------------------------
c     Argument parsing variables
c-----------------------------------------------------------------------
      integer iargc, i4arg
      real*8     r8arg
      real*8     r8_never
      parameter (r8_never = -1.0d-60)

c-----------------------------------------------------------------------
c     minst: minimum state in basis
c     maxst: maximum "
c     nst : number of states in basis
c     E: vector containing energy eigenvalues of H0
c     x: array containing dipole matrix elements of H0
c-----------------------------------------------------------------------
      integer minst, maxst
      parameter (minst = 1, maxst = 400)
      integer nst
      parameter (nst = maxst - minst + 1)
      real*8     e(minst:maxst)
      real*8     x(minst:maxst,minst:maxst)
c-----------------------------------------------------------------------
c     nc: Number of coefficients in the expansion of the wavefunction
c           = 2*Number of states in basis
c     ist,jst,kst: state indices for do loops
c     neq: vector containing system dimensions (used in fcn.f)
c     ic: coefficient index for do loops
c-----------------------------------------------------------------------
      integer nc, ist, jst, kst, ic, jc
      parameter (nc = 2*nst)
      integer neq(4)

c-----------------------------------------------------------------------
c     Initial state at t = 0
c-----------------------------------------------------------------------
      integer inst
      parameter (inst = 60)

c-----------------------------------------------------------------------
c     Maximum number of output times
c-----------------------------------------------------------------------
      integer maxout, p2max, p2
      parameter (p2max = 14)
      parameter (maxout = 2**p2max)

c-----------------------------------------------------------------------
c     ncyc: Number of field cycles over which to compute accel
c     hmax: Maximum harmonic seen in radiation spectrum
c     tout: Vector holding output times
c     ntout: Total number of output times
c-----------------------------------------------------------------------
      integer ncyc, hmax
      parameter (ncyc = 256)
      real*8     tout(maxout), work(maxout)
      integer ntout, itout

c-----------------------------------------------------------------------
c     lsoda tolerance
c-----------------------------------------------------------------------
      real*8     tol, default_tol
      parameter (default_tol = 1.0d-8)

c-----------------------------------------------------------------------
c     phi: phase angle of initial state
c     pi: 3.14...
c-----------------------------------------------------------------------
      real*8      phi, pi
      parameter (phi = 0.0d0)

c-----------------------------------------------------------------------
c     Logical variables to determine tracing
c-----------------------------------------------------------------------
      logical ltrace
      parameter (ltrace = .true.)

      logical lsodatrace
      parameter (lsodatrace = .false.)

      logical strace
      parameter (strace = .false.)

      logical ptrace
      parameter (ptrace = .true.)

      logical vartrace
      parameter (vartrace = .true.)

c-----------------------------------------------------------------------
c     Variables for output of wavefunction, if (ptrace)
c     
c     fname: Name of file in which the data for the wavefunction will
c            be stored (it will have a .hdf ending added)
c     pst: array storing occupation probabilities for energy states
c     npout: # of times for which pst will be written to output
c            (npout should always be a power of 2)
c     pint: pst written every pint'th step
c-----------------------------------------------------------------------
      character*(*) fname
      parameter (fname = 'wvfcn')

      integer pint, npout
      parameter (npout = 2 048)

      real*8     pst(nst)
      integer shape, rank
      real*8 bbox(2)

c-----------------------------------------------------------------------
c     varfile: file to which variances will be written
c     navg,n2avg: 1st and 2nd moments of distribution
c     var: variance of distribution in n (energy level index)
c-----------------------------------------------------------------------
      character*(*) varfile
      parameter (varfile = 'varn.dat')

      integer ufrom, getu
      real*8 navg, var

c-----------------------------------------------------------------------
c     Efile: file from which energy eigenvalues will be read
c-----------------------------------------------------------------------
      character*(*) Efile
      parameter (Efile = 'trianglE')

c-----------------------------------------------------------------------
c     lsoda variables
c-----------------------------------------------------------------------
      external fcn, jac

      real*8     y(nst*(nst + 3))
      real*8     yprime(nc)
      real*8     tbgn, tend
      integer itol
      real*8     rtol, atol
      integer itask, istate, iopt
      integer lrw
      parameter (lrw = 22 + 9*nc + nc**2)
      real*8     rwork(lrw)
      integer liw
      parameter (liw = 20 + nc)
      integer iwork(liw)
      integer jt

c-----------------------------------------------------------------------
c     tdacc(j,i): coefficient of time dependent term in acceleration
c                 'matrix element' (multiply by V*z(t)*cos(wt+theta))
c     accti(j,i): time independent term "
c     accre: real part of the expectation value of acceleration
c     accim: imaginary part (should be 0)
c     prob: total probability in the wavefunction (should be 1)
c-----------------------------------------------------------------------
      real*8     tdacc(minst:maxst,minst:maxst)
      real*8     accti(minst:maxst,minst:maxst)
      real*8     accre, accim
      real*8     prob

c-----------------------------------------------------------------------
c     Set parameters definied in common block
c     V0: strength of triangle well
c     L: width of triangle well
c     m: mass of quantum particle
c     hbar: obvious
c     theta: initial phase of radiation field
c-----------------------------------------------------------------------
      v0 = 4.0d-1
      l = 6.0d2
      m = 1.0d0
      hbar = 1.0d0
      theta = 0.0d0

c-----------------------------------------------------------------------
c     Set parameters for output of wavefunction if (ptrace)
c-----------------------------------------------------------------------
      bbox(1) = minst
      bbox(2) = maxst
      rank = 1
      shape = nst

c-----------------------------------------------------------------------
c     Argument parsing and definition of pi.
c     w: frequency of the radiation field
c     V: strength of the radiation field
c     hmax: maximum harmonic in radiation spectrum
c     tocyc: # of cycles in turn-on
c     tol: optional tolerance for lsoda
c-----------------------------------------------------------------------
      if (iargc() .lt. 4) 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

      hmax = i4arg(3, -1)
      tocyc = i4arg(4, -1)
      if (hmax .lt. 0 .or. tocyc .lt. 0) go to 900

      tol = r8arg(5, default_tol)

      pi = 4.0d0*atan(1.0d0)

c-----------------------------------------------------------------------
c     Set values for totime, ntout.
c     totime: total time for turn-on of laser field
c-----------------------------------------------------------------------
      totime = 2.0d0*pi*tocyc/w

      ntout = 2*hmax*ncyc

      if (ntout .gt. maxout) then
         write(0,*) 'Number of output times exceeds internal'//
     &        ' storage. Use lower hmax or change source code.'
         stop
      end if

c-----------------------------------------------------------------------
c     Check that ntout .gt. npout and set value of pint
c-----------------------------------------------------------------------
      if (ntout .lt. npout) then
         write(0,*) 'ntout .lt. npout! (modify npout in source code)'
         stop
      else
         pint = ntout/npout
      end if
      
c-----------------------------------------------------------------------
c     Check to see if ntout is a power of 2
c-----------------------------------------------------------------------
      do p2 = 1, p2max
         if (ntout .eq. 2**p2) go to 10
      end do

      write(0,*) 'Number of output times should be a power of 2!'
      write(0,*) 'Make sure that hmax and ncycle are powers of 2.'
      stop

 10   continue

c-----------------------------------------------------------------------
c     Set values in neq(4) for use with subroutine 'fcn'.
c-----------------------------------------------------------------------
      neq(1) = nc
      neq(2) = minst
      neq(3) = maxst
      neq(4) = nst

c-----------------------------------------------------------------------
c     Calculate and store output times
c-----------------------------------------------------------------------
      do itout = 1, ntout
         tout(itout) = 0.0d0 + (itout - 1)*pi/(w*hmax)
      end do

c-----------------------------------------------------------------------
c     Read in vector of energy eigenvalues using subroutine getE
c-----------------------------------------------------------------------
      call getE(Efile,minst,maxst,E)

c-----------------------------------------------------------------------
c     Read in dipole matrix elements using subroutine dipole
c-----------------------------------------------------------------------
      call dipole(minst,maxst,E,x)

c-----------------------------------------------------------------------
c     Tracing of getE and dipole subroutine calls
c-----------------------------------------------------------------------
      if (strace) then
         write(0,*) 'minst = ',minst,', maxst = ',maxst
         do ist = minst, maxst
            write(0,*) 'E(',ist,') = ',E(ist)
         end do
         write (0,*) 'dipole elements'
         do ist = minst, maxst
            write(0,450) (x(ist,jst), jst = minst, maxst)
 450        format(1p,5e15.6)
         end do
      end if

c-----------------------------------------------------------------------
c     Calculate accti(j,i) and tdacc(j,i)
c-----------------------------------------------------------------------
      do ist = minst, maxst
         do jst = minst, maxst
            accti(jst,ist) = -((E(jst)-E(ist))**2)*x(jst,ist)/
     &           hbar**2
            tdacc(jst,ist) = 0.0d0
            
            do kst = minst, maxst
               tdacc(jst,ist) = tdacc(jst,ist) + (2.0d0*e(kst)-
     &              E(ist)-E(jst))*x(jst,kst)*x(kst,ist)/hbar**2
            end do

         end do
      end do

c-----------------------------------------------------------------------
c     Initialize the coefficients of the wavefunction
c     Note: y(ic) = a(ic+maxst-nst)
c           y(ic+nst) = b(ic+maxst-nst)
c-----------------------------------------------------------------------
      do ic = 1, nc
         y(ic) = 0.0d0
      end do

      y(inst-maxst+nst) = cos(phi)
      y(inst-maxst+2*nst) = sin(phi)

c-----------------------------------------------------------------------
c     Store E and x as elements of the vector y.  This is necessary to
c     pass the values of E and x to the subroutine 'fcn'.
c     y(ic+2*nst) = E(ic+maxst-nst)
c     y(nst*ic+jc+2*nst) = x(jc+maxst-nst,ic+maxst-nst)
c-----------------------------------------------------------------------
      do ic = 1, nst
         y(ic+2*nst) = E(ic+maxst-nst)
      end do

      do ic = 1, nst
         do jc = 1, nst
            y(nst*ic+jc+2*nst) = x(jc+maxst-nst,ic+maxst-nst)
         end do
      end do

c-----------------------------------------------------------------------
c     Calculate accre, accim, and prob at 
c     t = tout(1) = 0. 
c-----------------------------------------------------------------------
      if (tout(1) .ge. totime) then
         accre = tdacc(inst,inst)*v*cos(theta)
      else
         accre = 0.0d0
      end if
    
      accim = 0.0d0
      prob = 1.0d0

      if (ltrace) then
         write(0,500) tout(1), prob, accim
 500     format(' ','At step 1, t = ',1pe10.3,', prob = ',1pe12.5,
     &        ', accim = ',1pe12.5)
      end if

      write(*,*) tout(1), accre

c-----------------------------------------------------------------------
c     Tracing output of probabilities in each state to file fname
c-----------------------------------------------------------------------
      if (ptrace) then
         do ic = 1, nst
            pst(ic+maxst-nst) = y(ic)*y(ic) + y(ic+nst)*y(ic+nst)
         end do
         call gft_write_bbox(fname,tout(1),shape,rank,bbox,pst)
      end if

c-----------------------------------------------------------------------
c     Calculate and write var at initial time
c-----------------------------------------------------------------------
      if (vartrace) then
         ufrom = getu()
         open(ufrom,file=varfile,status='new')
         navg = 0.0d0
         var = 0.0d0
         do ic = 1, nst
            navg = navg + (ic+maxst-nst)*(y(ic)*y(ic) + 
     &           y(ic+nst)*y(ic+nst))
         end do
         write(0,*) 'At t = ',tout(1),' navg = ',navg
         do ic = 1, nst
            var = var + (real(ic+maxst-nst)-navg)*
     &           (real(ic+maxst-nst)-navg)*(y(ic)*y(ic) +
     &           y(ic+nst)*y(ic+nst))
         end do
         write(ufrom,*) tout(1)*w/(2.0d0*pi), var
      end if

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

c=======================================================================
c     BEGIN DO LOOP OVER TIME VALUES
c=======================================================================
      do itout = 2, ntout
         istate = 1

c-----------------------------------------------------------------------
c     Need these temporaries since lsoda overwrites tend ...
c-----------------------------------------------------------------------
         tbgn = tout(itout - 1)
         tend = tout(itout)

c-----------------------------------------------------------------------
c     Call ODE integrator lsoda.  Returns solution at t = tend in
c     y(neq) (i.e. y is overwritten).
c-----------------------------------------------------------------------
         call lsoda(fcn,neq,y,tbgn,tend,itol,rtol,atol,itask,
     &        istate,iopt,rwork,lrw,iwork,liw,jac,jt)

         if (lsodatrace) then
            write(0,1000) itout, ntout, tout(itout-1), tout(itout)
 1000       format(' ','accel: Step ',i4,' of ',i4,' t = ',
     &           1pe10.3,' .. ', 1pe10.3)
            write(0,*) 'accel: lsoda returns ',istate
         end if

c-----------------------------------------------------------------------
c     Error message if lsoda fails
c-----------------------------------------------------------------------
         if (istate .lt. 0) then
            write(0,1500) istate, itout, ntout, tout(itout-1),
     &           tout(itout)
 1500       format(/' ','accel: Error retrun ',i2,
     &           ' from integrator lsoda.'/
     &           '   At output time ',i5,' of ',i5/
     &           '   Interval ',1pe10.3,' .. ',1pe10.3/)
            stop
         end if

c-----------------------------------------------------------------------
c     Tracing output of probabilities p(ist) to file fname
c-----------------------------------------------------------------------
         if (ptrace .and. (mod(ntout,pint) .eq. 0)) then
            do ic = 1, nst
               pst(ic+maxst-nst) = y(ic)*y(ic) + y(ic+nst)*y(ic+nst)
            end do
            call gft_write_bbox(fname,tout(itout),shape,rank,bbox,pst)
         end if

c-----------------------------------------------------------------------
c     Output of variance to file varfile
c-----------------------------------------------------------------------
         if (vartrace) then
            navg = 0.0d0
            var = 0.0d0
            do ic = 1, nst
               navg = navg + (ic+maxst-nst)*(y(ic)*y(ic) +
     &              y(ic+nst)*y(ic+nst))
            end do
            write(0,*) 'At t = ',tout(itout),' navg = ',navg
            do ic = 1, nst
               var = var + (real(ic+maxst-nst)-navg)*
     &              (real(ic+maxst-nst)-navg)*(y(ic)*y(ic) +
     &              y(ic+nst)*y(ic+nst))
            end do
            write(ufrom,*) tout(itout)*w/(2.0d0*pi), var
         end if

c-----------------------------------------------------------------------
c     Calculate total probability at time t = tout(itout)
c-----------------------------------------------------------------------
         if (ltrace) then
            prob = 0.0d0

            do ic = 1, nc
               prob = prob + y(ic)*y(ic)
            end do
         end if

c-----------------------------------------------------------------------
c     Calculate accre and accim at time t = tout(itout)
c-----------------------------------------------------------------------
         accre = 0.0d0
         accim = 0.0d0

         if  (tout(itout) .lt. totime) then
            do ic = 1, nst
               do jc = 1, nst
                  accre = accre + 
     &                 (y(jc)*y(ic) + y(jc+nst)*y(ic+nst))*
     &                 (accti(jc+maxst-nst,ic+maxst-nst) 
     &                 + tdacc(jc+maxst-nst,ic+maxst-nst)*v*
     &                 (sin(w*tout(itout)/(4.0d0*tocyc))**2)*
     &                 cos(w*tout(itout)+theta))
                  accim = accim + 
     &                 (y(jc)*y(ic+nst) - y(jc+nst)*y(ic))*
     &                 (accti(jc+maxst-nst,ic+maxst-nst) 
     &                 + tdacc(jc+maxst-nst,ic+maxst-nst)*v*
     &                 (sin(w*tout(itout)/(4.0d0*tocyc))**2)*
     &                 cos(w*tout(itout)+theta))
               end do
            end do
         else
            do ic = 1, nst
               do jc = 1, nst
                  accre = accre + 
     &                 (y(jc)*y(ic) + y(jc+nst)*y(ic+nst))*
     &                 (accti(jc+maxst-nst,ic+maxst-nst) 
     &                 + tdacc(jc+maxst-nst,ic+maxst-nst)*v*
     &                 cos(w*tout(itout)+theta))
                  accim = accim + 
     &                 (y(jc)*y(ic+nst) - y(jc+nst)*y(ic))*
     &                 (accti(jc+maxst-nst,ic+maxst-nst) 
     &                 + tdacc(jc+maxst-nst,ic+maxst-nst)*v*
     &                 cos(w*tout(itout)+theta))
               end do
            end do
         end if


c-----------------------------------------------------------------------
c     Tracing output: prob and accim
c-----------------------------------------------------------------------
         if (ltrace) then

            write(0,2000) itout, tout(itout), prob, 
     &           accim
 2000       format(' ','At step ',i5,', t = ',1pe10.3,', prob = ',
     &           1pe12.5,', accim = ',1pe12.5)

         end if

c-----------------------------------------------------------------------
c     Write time series to standard output and stop
c-----------------------------------------------------------------------
         write(*,*) tout(itout), accre
 
c-----------------------------------------------------------------------
c     END DO LOOP OVER OUTPUT TIMES
c-----------------------------------------------------------------------
      end do
      
      if (ptrace) then
         call gft_close_all()
      end if

      if (vartrace) then
         close(ufrom)
      end if

      stop

c-----------------------------------------------------------------------
c     Usage message
c-----------------------------------------------------------------------
 900  continue
      write(0,*) 'usage: accel <w> <V> <hmax> <tocyc> [<tol>]'
      write(0,*)
      write(0,*) 'w: frequency of driving field'
      write(0,*) 'V: strength of driving field'
      write(0,*) 'hmax: maximum harmonic order'
      write(0,*) 'tocyc: # of cycles in turn-on'
      write(0,*) 'tol: optional tolerance for lsoda'
      write(0,*)
      write(0,*) 'Output written to standard out'
      write(0,*) 'Probabilities written to ',fname,' if (ptrace)'
      write(0,*) 'Energies read from ',Efile

      stop

      end