c=======================================================================
c     This is a subroutine that is used with the 'accel' program
c     (accel.f).  It is used to evaluate the derivatives of the
c     coefficients in the expansion of a quantum state with 
c     Hamiltonian:
c
c     H = H0 - V*x*cos(w*t + theta)
c
c     where H0 is a time-independent, integrable Hamiltonian.
c
c     The wavefunction is expanded in the energy basis of H0.
c=======================================================================
      subroutine fcn(neq, t, y, yprime)

         implicit none

c-----------------------------------------------------------------------
c     Common block for communication with accel.f, getE.f, dipole.f
c-----------------------------------------------------------------------
         include 'fcn.inc'

c-----------------------------------------------------------------------
c     Arguments of subroutine and indices
c     Note: dimension of y is nst*(nst+3) (nst = dims(4))
c           dimension of yprime is neq = dims(1)
c-----------------------------------------------------------------------
         integer ic, jc
         integer minst, maxst, nst
         real*8 t
         integer neq
         real*8 y, yprime

         dimension neq(4)
         dimension y(1), yprime(1)

c-----------------------------------------------------------------------
c     Temporaries for dimension of system
c-----------------------------------------------------------------------
         minst = neq(2)
         maxst = neq(3)
         nst = neq(4)

c-----------------------------------------------------------------------
c     Define derivatives
c     Note: y(ic) = a(ic+maxst-nst)
c           y(ic+nst) = b(ic+maxst-nst)
c           y(ic+2*nst) = E(ic+maxst-nst)
c           y(nst*ic+jc+2*nst) = x(jc+maxst-nst,ic+maxst-nst)
c
c     Note: we will replace x(j,i) with x(i,j) in the calculations
c     below for better vectorization.  We can do this as long as
c     the x(j,i) are real.
c-----------------------------------------------------------------------
         if (t .lt. totime) then

            do ic = 1, nst
               yprime(ic) = y(ic+2*nst)*y(ic+nst)/hbar

               yprime(ic+nst) = - y(ic+2*nst)*y(ic)/hbar

               do jc = 1, nst
                  yprime(ic) = yprime(ic) -
     &                 V*(sin(w*t/(4.0d0*tocyc))**2)*
     &                 cos(w*t+theta)*
     &                 y(nst*ic+jc+2*nst)*y(jc+nst)/hbar

                  yprime(ic+nst) = yprime(ic+nst)
     &                 +V*(sin(w*t/(4.0d0*tocyc))**2)*
     &                 cos(w*t+theta)*
     &                 y(nst*ic+jc+2*nst)*y(jc)/hbar
               end do
            end do

         else

            do ic = 1, nst
               yprime(ic) = y(ic+2*nst)*y(ic+nst)/hbar

               yprime(ic+nst) = - y(ic+2*nst)*y(ic)/hbar

               do jc = 1, nst
                  yprime(ic) = yprime(ic) -
     &                 V*cos(w*t+theta)*y(nst*ic+jc+2*nst)
     &                 *y(jc+nst)/hbar

                  yprime(ic+nst) = yprime(ic+nst) +
     &                 V*cos(w*t+theta)*y(nst*ic+jc+2*nst)
     &                 *y(jc)/hbar
               end do

            end do
         end if

         return

      end

c=======================================================================
c     Optional Jacobian for lsoda (see accel.f)
c=======================================================================
      subroutine jac(neq,t,y,ml,mu,pd,nrowpd)

         implicit none

         include 'fcn.inc'

         integer ic, jc
         integer minst, maxst, nst
         real*8 t
         integer neq(4)
         real*8 y, pd
         integer ml, mu, nrowpd

         dimension y(1), pd(nrowpd, 1)

c----------------------------------------------------------------------
c     Set values of minst, maxst, nst from array neq
c-----------------------------------------------------------------------
         minst = neq(2)
         maxst = neq(3)
         nst = neq(4)

c-----------------------------------------------------------------------
c     Calculate Jacobian.  Note:  this routine only works because
c     the x(j,i) are real and thus x(j,i) = x(i,j).  If the x(j,i) are
c     complex then the entier "accel" program must be altered.
c-----------------------------------------------------------------------
         if (t .lt. totime) then

c-----------------------------------------------------------------------
c     If we are still in turn-on... calculate bottom left of jac
c-----------------------------------------------------------------------
            do ic = 1, nst
               do jc = nst+1, 2*nst
                  if (jc .eq. ic+nst) then
                     pd(jc,ic) = -y(ic+2*nst)/hbar + V*(sin(w*t/
     &                    (4.0d0*tocyc))**2)*cos(w*t+theta)*
     &                    y((nst+1)*ic+2*nst)/hbar
                  else
                     pd(jc,ic) = V*(sin(w*t/(4.0d0*tocyc))**2)*
     &                    cos(w*t+theta)*y(nst*ic+jc+nst)/hbar
                  end if
               end do
            end do

c-----------------------------------------------------------------------
c     and calculate top right of jac
c-----------------------------------------------------------------------
            do ic = nst+1, 2*nst
               do jc = 1, nst
                  pd(jc,ic) = -pd(ic,jc)
               end do
            end do

         else

c-----------------------------------------------------------------------
c     If we are out of turn-on... calculate bottom left of jac
c-----------------------------------------------------------------------
            do ic = 1, nst
               do jc = nst+1, 2*nst
                  if (jc .eq. ic+nst) then
                     pd(jc,ic) = -y(ic+2*nst)/hbar + V*
     &                    cos(w*t+theta)*y((nst+1)*ic+2*nst)/hbar
                  else
                     pd(jc,ic) = V*cos(w*t+theta)*
     &                    y(nst*ic+jc+nst)/hbar
                  end if
               end do
            end do

c-----------------------------------------------------------------------
c     and calculate top right of jac
c-----------------------------------------------------------------------
            do ic = nst+1, 2*nst
               do jc = 1, nst
                  pd(jc,ic) = -pd(ic,jc)
               end do
            end do
         end if

         return

      end