c=======================================================================
c     diffp.f: This routine generates a graph of the difference
c     between the outputs of 'dhonqp.f' and 'dhoqp.f'.
c=======================================================================
      program dhonqp

      implicit none

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

      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
      w0 = 4.0d0
      V = 1.0d-1
      m = 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 density of states with the two methods and write the
c     difference to standard out
c-----------------------------------------------------------------------
      do i = 1, imax
         Q = 0.0d0 + dble(i-1)*w*hbar/dble(imax-1)
         b = 0.0d0
         c = 0.0d0
         do j = 1, jmax
            b = b - T*cos((dble(j)*Q*T+S(j))/hbar)/(2.0d0*pi*
     &           hbar*sqrt(R(j)))
            a = 2.0d0*pi*j*(Q-V**2/(4.0d0*m*(w**2-w0**2)))/
     &           (w*hbar)
            c = c - 2.0d0*cos(a)/(w*hbar*sqrt(2.0d0*(1.0d0
     &           - cos(2.0d0*pi*j*w0/w))))
         end do
         write(*,*) Q, b-c
      end do
      
      stop
      
      end