c   This is a modification which allows for continuous values of omega.
c
C     THIS PROGRAM COMPUTES THE Q'S FOR ZERO TEMPERATURE STATES
C     YOU COMPUTE THEM FOR A NUMBER OF VALUES OF S GIVEN BY NPTS
C     DELS DETERMINES THE SPACING BETWEEN THE POINTS
C     NPTS MUST BE AT LEAST 2 FOR THE PROGRAM TO WORK
	IMPLICIT NONE
	Integer nn,nl,nt,npts,nj,ni,nf,li,lf,nom,Junk,nab
	Real*8 dels,eps,h1,hmin,si,
     1         xi,xim1,sf,omi,delom
        Real*8 a2,a3,a4,a5,a6,a7,a8,b2,b3,b4,b5,b6,b7,b8
        Real*8 A22,A33,A34,A35,A44,B22,B33,B34,B35,B44
	Real*8 sstart
        Common/coeff/a2,a3,a4,a5,a6,a7,a8,b2,b3,b4,b5,b6,b7,b8
        Common/series/A22,A33,A34,A35,A44,B22,B33,B34,B35,B44
	COMMON/DIFEQ/DELS,EPS,H1,HMIN,NPTS
	COMMON/PAR/xi,omi,delom,nom
	Open(13,File='s.dat',status='old')
	OPEN(14,FILE='ode.dat',STATUS='OLD')
	OPEN(15,FILE='modes.dat',STATUS='OLD')
	OPEN(16,FILE='series.dat',STATUS='OLD')
	Open(17,file='qinit.out')
	OPEN(20,FILE='qmodesz.out')
c  si is the value of s you want to start computing the stress tensor at
        Read(13,*) Junk
        Read(13,*) si,dels,npts
	Read(14,*) Junk
	READ(14,*) EPS,H1,HMIN
	Read(15,*) Junk
	READ(15,*) omi,delom,nom,LI,LF
	Read(15,*) Junk
c  sstart is for pmodesz.  sf is the starting value for s for the q modes.
	Read(15,*) sstart,sf
	Print *,'sf = ',sf
	Read(16,*) Junk
	READ(16,*) xim1,nab
	If (dabs(xim1).gt.1.d10) then
	   xi = 0.d0
	Else
	   xi = 1.d0/xim1
	Endif
	Read(16,*) Junk
        Read(16,*) a2,b2
        If(nab.ge.2) then
          Read(16,*) junk
          Read(16,*) a3,b3
        Endif
        If(nab.ge.3) then
          Read(16,*) junk
          Read(16,*) a4,b4
        Endif
c  Now convert to have f and k in terms of An and Bn.
        Call Convert(nab)
c  This labels the file so as to identify it uniquely
        Write(20,*) '0     xi        a2         b2         a3        b3      a4    b4'
        Write(17,*) '0     xi        a2         b2         a3        b3      a4    b4'
        Write(20,3) xi,a2,b2,a3,b3,a4,b4
        Write(17,3) xi,a2,b2,a3,b3,a4,b4
        Write(20,*) '0     NAB       A22       B22         A33       B33,   A44    B44'
        Write(17,*) '0     NAB       A22       B22         A33       B33    A44    B44'
        Write(20,4) nab,A22,B22,A33,B33,A44,B44
        Write(17,4) nab,A22,B22,A33,B33,A44,B44
   3    Format(5x,9(F9.5,3x))
   4    Format(5x,I3,7x,9(F9.5,3x))
   1    FORMAT(1X,F7.3,2x,I3,2x,7(F7.3,2x))
	CALL EVOLVQ(si,SF,LI,LF)
	STOP
	END

  		SUBROUTINE EVOLVQ(si,SF,LI,LF)
	IMPLICIT None
	Integer nn,NL,nt,n1,nom,L,Li,Lf,LL,kmax,kount,I,J,npts,nbad,nok,N,iseries
	Real*8 om,XL,DXSAV, XP(200), YP(10,200)
	Real*8 QI,QPI,Q,QP,x2,xnpts,x1,Y,DYDX,difq
	Real*8 DELS,EPS,H1,HMIN,si,sf,DELX,f,fp,f2p,k,kp
	Real*8 xi,omi,delom,capom,capomp,capom2p,r,q2p
	Real*8 errest
        Real*8 a2,a3,a4,a5,a6,a7,a8,b2,b3,b4,b5,b6,b7,b8
        Common/coeff/a2,a3,a4,a5,a6,a7,a8,b2,b3,b4,b5,b6,b7,b8
        common/series/iseries
	DIMENSION Y(2),DYDX(2)
	COMMON /DER/om,XL/PATH/ KMAX, KOUNT, DXSAV, XP, YP
	COMMON /DIFEQ/DELS,EPS,H1,HMIN,NPTS
	COMMON/PAR/xi,omi,delom,nom
	EXTERNAL DERIVS
	EXTERNAL BSSTEP
C	WE WANT TO EVOLVE BACKWARDS
c       This is for the power series routine for omega = 0
	iseries = 0
	DELX = -DELS
	DO 1000 N=1,nom
	N1 = N + 1
	om = omi+(n-1)*delom
	DO 1000 L=LI,LF
	LL = L+1-LI
	XL = L
c  Use a power series solution for the initial value of q if om = 0.
	If(dabs(om).lt.1.d-12) then
	  r = sf+1.d0
	  Call Qinit(r,xl,q,qp,q2p)
	Endif
c  If iseries is equal to 1 then we don't want to use the series because the power alpha is
c   complex.
	If(iseries.eq.1) then
c  Let's check the accuracy of the power series
	  Y(1) = q
	  Y(2) = qp
	  qi = q
	  qpi = qp
	  Call Derivs(sf,Y,DYDX,2)
	  difq = dabs((DYDX(2)-q2p)/DYDX(2))
	  Print*,'q = ',q,' qp = ',qp,' q2p = ',q2p
	  Print*,'DYDX = ',DYDX(2),' difq = ',difq
	  iseries = 0
        Else
c   This gives the ratio of q'/q for the zeroth order WKB approximation
          Call Metric(sf,f,fp,f2p,k,kp)
c         Capom = dsqrt(b2*om**2+xl*(xl+1.d0)*f/(sf+1.d0)**2)
          Capom = dsqrt(om**2+(xl+0.5d0)**2*f/(sf+1.d0)**2)
	  Capomp = (xl+0.5d0)**2*(fp/(sf+1.d0)**2-2.d0*f/(sf+1.d0)**3)/(2.d0*Capom)
	  Capom2p = -(xl+0.5d0)**4*(fp/(sf+1.d0)**2-2.d0*f/(sf+1.d0)**3)**2/(4.d0*Capom**3)
     1        +(xl+0.5d0)**2*(f2p/(sf+1.d0)**2 -4.d0*fp/(sf+1.d0)**3 +6.d0*f/(sf+1.d0)**4)/
     2         (2.d0*Capom)
c	  Print *,'Capom = ',Capom,' Capomp = ',Capomp,' Capom2p = ',Capom2p
c	  Print*,'b2 = ',b2,' om = ',om
c	  Print*,'xl = ',xl,' f = ',f
c	  Print*,'sf = ',sf,' fp = ',fp
c	  Print*,'f2p = ',f2p
          Y(2) = -Capom*dsqrt(1.d0/(f*k)) - Capomp/(2.d0*Capom)
	  Y(1) = 1.d0
	  q = Y(1)
	  qp = Y(2)
c  This is a very rough estimate
	  q2p = Capom**2/(f*k) + Capom*(fp/(2.d0*f**1.5d0*k**0.5d0) + kp/(2.d0*f**0.5d0*k**1.5d0))
     1          - Capom2p/(2.d0*Capom) + 0.75d0*Capomp**2/Capom**2
	  Call Derivs(sf,Y,DYDX,2)
	  difq = dabs((DYDX(2)-q2p)/DYDX(2))
c	  Print*,'q = ',q,' qp = ',qp,' q2p = ',q2p
c	  Print*,'DYDX = ',DYDX(2),' difq = ',difq
        Endif
	KMAX = 0
	NOK = 0
	NBAD = 0
	X1 = SF
C        GO TO ONE STEP LARGER THAN WHERE YOU WANT TO START
	XNPTS = NPTS 
	X2 = si + XNPTS*DELS
	CALL ODEINT(Y,2,X1,X2,EPS,H1,HMIN,NOK,NBAD,DERIVS,BSSTEP)
	q = Y(1)
	qp = Y(2)
c  Check to see if qp is too large
	If(dabs(qp).gt.1.d299) then
          Print*,' qp is too large.  omega = ',om,' L = ',L,' qp = ',qp
          Stop
        Endif
c  Now check to see if the total accuracy is OK.  Note this is a double precision version
c    of the code so we'll require 14 digits of accuracy.
        If(difq.lt.1.d-15) then
          errest = 14.d0 + dlog10(q)
        Else
          errest = -dlog10(difq) + dlog10(q)
        Endif
	WRITE(6,12) om,L,sf,X2,q,qp,difq,errest
  12	FORMAT(1X,d12.5,2X,I4,2X,5(D15.7),2x,F5.0)
        If(errest.lt.14.d0) then
          Print *,' Too much error in the above mode. '
          Stop
        Endif
	Write(17,45) sf,om,L,q,qp,difq,errest
  45    Format(1x,2(d12.5,2x),I4,2x,3(d12.5,2x),2x,F5.0)
	X1 = X2
	DO 980 I=1,NPTS
	X2 = X1 + DELX
	CALL ODEINT(Y,2,X1,X2,EPS,H1,HMIN,NOK,NBAD,DERIVS,BSSTEP)
	Q = Y(1)
	QP = Y(2)
c  Check to see if qp is too large
	If(dabs(qp).gt.1.d299) then
          Print*,' qp is too large.  omega = ',om,' L = ',L,' qp = ',qp
          Stop
        Endif
	J = npts+1-I
	Write(20,*) om,L,J,x2
	Write(20,*) Q,QP
 980	X1 = X2
1000    CONTINUE
	WRITE(6,14) X1
  14	FORMAT(1X,'si FROM EVOLVQ =', F10.7)
	RETURN
 	END
 
 	SUBROUTINE DERIVS(X,Y,DYDX,NVAR)
	IMPLICIT None
	Integer NVAR,nom
	Real*8 om,XL,xi,omi,delom
	Real*8 X,Y,DYDX,s,r,f,fp,f2p,k,kp,q,qp,ScalarR
	DIMENSION Y(NVAR),DYDX(NVAR)
	COMMON/DER/om,XL
	COMMON/PAR/xi,omi,delom,nom
	s = X
	r = s+1.d0
	Call metric(s,f,fp,f2p,k,kp)
	q = Y(1)
	qp = Y(2)
	ScalarR = -k*f2p/f + k*fp**2/(2.*f**2) - fp*kp/(2.*f) 
     1         - 2.*k*fp/(r*f)
     1       - 2.*kp/r - 2.*k/r**2 + 2.d0/r**2       
	DYDX(1) = qp
	DYDX(2) = -(2.d0/r + fp/(2.d0*f) + kp/(2.d0*k))*qp 
     1               + (om**2/(k*f) 
     1             + xl*(xl+1.d0)/(k*r**2) + xi*ScalarR/k)*q
c	Print*,'From DERIVS:  r = ',r,' f = ',f,' fp = ',fp
c	Print*,' f2p = ',f2p,' k = ',k,' kp = ',kp
c	Print*,'ScalarR = ',ScalarR,' q = ',q,' qp = ',qp
c	Print*,'omega = ',om,' qpp = ',DYDX(2)
 	RETURN
	END