c=======================================================================
c     dhonqp.f: This routine generates a plot of the quasienergy
c     density of states for the driven harmonic oscillator, using
c     results from the program 'orbitnt.f'.  
c=======================================================================
      program dhonqp

      implicit none

c-----------------------------------------------------------------------
c     w: driving frequency
c     T: driving period
c     Q: quasienergy value
c     R(j): residue of the jth periodic orbit
c     S(j): action "
c     imax: total number of output points
c     jmax: number of orbits used in calculating density of states
c-----------------------------------------------------------------------
      real w, pi, T, Q, hbar, dQ, Rt, St, b

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

      real R(jmax), S(jmax)

      pi = 4.0d0*atan(1.0d0)
      w = 4.0d0*pi
      T = 2.0d0*pi/w
      hbar = 1.0d0

c-----------------------------------------------------------------------
c     read in output from 'orbitnt.f'
c-----------------------------------------------------------------------
      do j = 1, jmax
         read(*,*) N, Rt, St
         if (N .eq. j) then
            R(j) = Rt
            S(j) = St
         else
            write(0,*) 'Error reading input'
         end if
      end do

c-----------------------------------------------------------------------
c     calculate the density of states for each value of Q and output
c     to standard out
c-----------------------------------------------------------------------
      do i = 1, imax
         Q = 0.0d0 + dble(i-1)*w*hbar/dble(imax-1)
         b = 0.0d0
         do j = 1, jmax
            b = b - T*cos((dble(j)*Q*T+S(j))/hbar)/(2.0d0*pi*
     &           hbar*sqrt(R(j)))
         end do
         write(*,*) Q, b
      end do
      
      stop
      
      end