c=======================================================================
c     dhoqp.f: This routine generates a plot of the quasienergy 
c     density of states for the driven harmonic oscillator using
c     the results of analytical calculation.
c=======================================================================
      program dhoqp

      implicit none

c-----------------------------------------------------------------------
c     w: driving frequency
c     w0: natural frequency of the oscillator
c     E: quasienergy value
c     V: strength of the driving field
c     imax: number of output points
c     jmax: number of periodic orbits used in calculating the density
c           of states
c-----------------------------------------------------------------------
      real w, w0, a, b, E, dE, V, pi, m, hbar

      integer i, imax, j, jmax
      parameter(imax = 5000, jmax = 60)

      pi = 4.0d0*atan(1.0d0)
      w = 4.0d0*pi
      w0 = 4.0d0
      V = 1.0d-1
      m = 1.0d0
      hbar = 1.0d0

c-----------------------------------------------------------------------
c     calculate the density of states for each E and write to 
c     standard out
c-----------------------------------------------------------------------
      do i = 1, imax
         E = 0.0d0 + dble(i-1)*w*hbar/dble(imax-1)
         b = 0.0d0
         do j = 1, jmax
            a = 2.0d0*pi*j*(E-V**2/(4.0d0*m*(w**2-w0**2)))/
     &           (w*hbar)
            b = b - 2.0d0*cos(a)/(w*hbar*sqrt(2.0d0*(1.0d0
     &           - cos(2.0d0*pi*j*w0/w))))
         end do
         write(*,*) E, b
      end do
      
      stop
      
      end