pmodesz saved on 2-19-2009

C    This program is a revision of pmodes.f in that it allows for a
c     continuous set of values of omega rather than a discrete set and
c     works for a more general zero temperature black hole metric.
c
C     THIS PROGRAM COMPUTES THE P'S  FOR CHARGED BLACK HOLES
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
C     THE IMPORTANT OUTPUT IS PUT IN PQDIF.OUT
	IMPLICIT NONE
	Integer nn,nl,npts,nj,ni,nf,li,lf,nom,Junk,nab
	Real*8 P,PP,PI,PPI,dels,eps,h1,hmin,xm,qm,sp,sm,xkap,si,xi,
     1         xim1,sf,delom,omi,sstart
        Real*8 A22,A33,A34,A35,A44,B22,B33,B34,B35,B44
        Real*8 a2,a3,a4,a5,a6,a7,a8,b2,b3,b4,b5,b6,b7,b8
	Real*8 ae2,be2
        Common/series/A22,A33,A34,A35,A44,B22,B33,B34,B35,B44
        Common/coeff/a2,a3,a4,a5,a6,a7,a8,b2,b3,b4,b5,b6,b7,b8
	COMMON/DIFEQ/DELS,EPS,H1,HMIN,NPTS
        COMMON/PAR/xi,omi,delom,nom
	common/erncheck/ae2,be2,nab
	OPEN(14,FILE='s.dat',STATUS='OLD')
	OPEN(15,FILE='ode.dat',STATUS='OLD')
	OPEN(16,FILE='modes.dat',STATUS='OLD')
	OPEN(17,FILE='series.dat',STATUS='OLD')
	OPEN(18,FILE='pinit.out')
	OPEN(19,FILE='pmodesz.out')
c  si is the value of s you want to start computing the stress tensor at
	Read(14,*) Junk
	Read(14,*) si,dels,npts
c  sstart is the value of s that you evaluate the power series at
	Read(15,*) Junk
	READ(15,*) EPS,H1,HMIN
	Read(16,*) Junk
	READ(16,*) omi,delom,nom,LI,LF
	Read(16,*) Junk
	Read(16,*) sstart
	Read(17,*) Junk
        READ(17,*) xim1,nab
        If (dabs(xim1).gt.1.d10) then
           xi = 0.d0
        Else
           xi = 1.d0/xim1
        Endif
	Print*,'xi = ',xi
	Read(17,*) junk
	Read(17,*) a2,b2
c  This is just because I need new names for a different common statement
	ae2 = a2
	be2 = b2
	If(nab.ge.2) then
	  Read(17,*) junk
	  Read(17,*) a3,b3
	Endif
	If(nab.ge.3) then
	  Read(17,*) junk
	  Read(17,*) 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(18,*) '0     xi        a2         b2         a3        b3     a4     b4'
	Write(19,*) '0     xi        a2         b2         a3        b3     a4     b4'
	Write(18,3) xi,a2,b2,a3,b3,a4,b4
	Write(19,3) xi,a2,b2,a3,b3,a4,b4
	Write(18,*) '0     NAB       A22       B22         A33       B33    A44    B44'
	Write(19,*) '0     NAB       A22       B22         A33       B33    A44    B44'
	Write(18,4) nab,A22,B22,A33,B33,A44,B44
	Write(19,4) nab,A22,B22,A33,B33,A44,B44
	Write(18,*) '0   Sstart         omega       L       p(si)        p1p(si)      difp2p     Acc'
   3    Format(5x,9(F9.5,3x))
   4    Format(5x,I3,7x,9(F9.5,3x))
	CALL EVOLVP(sstart,si,LI,LF)
	STOP
	END
 	SUBROUTINE EVOLVP(sstart,si,LI,LF)
	IMPLICIT NONE
        Integer i1,nn,NL,nt,n1,nom,L,Li,Lf,LL,kmax,kount,I,J,npts,nbad,nok,N,nab
        Real*8 om,XL,DXSAV, XP(200), YP(10,200)
        Real*8 p,p1p,p2p,x2,xnpts,x1,Y,DYDX,difp
        Real*8 DELS,EPS,H1,HMIN,si,sstart,DELX,f,fp,f2p,k,kp
        Real*8 xi,omi,delom,capom,r,errest
	Real*8 a2,b2
        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
	common/erncheck/a2,b2,nab
        EXTERNAL DERIVS
        EXTERNAL BSSTEP
	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 starting value
            If(dabs(om).lt.1.d-12) then
	       If(dabs(xi).le.1.d-12.and.L.eq.0) then
	         P = 1.d0
	         p1p = 0.d0
	         p2p = 0.d0
	         x2 = si
	         Go to 998
	       Endif
c  For ERN need to use analytic values if omega = 0.
	       If(nab.eq.1.and.dabs(a2-1.d0).lt.1.d-12.and.dabs(b2-1.d0).lt.1.d-12) then
	         P = si**L
	         p1p = L*si**(L-1)
	         p2p = L*(L-1)*si**(L-2)
	         x2 = si
	         Go to 998
	       Endif
               Call  PSeries0(sstart,L,p,p1p,p2p)
            else
               Call  PSeries(sstart,om,L,p,p1p,p2p)
	    Endif
c  Renormalize p so that it is equal to 1 at the starting value
c   This is so that I can easily see how much it has grown.
	    p1p = p1p/p
	    p2p = p2p/p
	    p = 1.d0
	    Y(1) = P
	    Y(2) = P1P
c  Let's check the accuracy of the power series
            Call Derivs(sstart,Y,DYDX,2)
            difp = dabs((DYDX(2)-p2p)/DYDX(2))
	    DXSAV = 0.d0
	    KMAX = 0
	    NOK = 0
	    NBAD = 0
	    X1 = sstart
	    x2 = si
            CALL ODEINT(Y,2,X1,X2,EPS,H1,HMIN,NOK,NBAD,DERIVS,BSSTEP)
	    P = Y(1)
	    P1P = Y(2)
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.
c  For \xi > 0 and omega = L = 0 the mode actually decreases because R is negative
c   So do not take that into account.
	    If(om.lt.1.d-12.and.L.eq.0) then
	      errest = -dlog10(difp)
	    Else if(difp.lt.1.d-15) then
	      errest = 14.d0 + dlog10(p)
	    Else
	      errest = -dlog10(difp) + dlog10(p)
	    Endif
            WRITE(6,12) om,L,sstart,si,p,p1p,difp,errest
  12    FORMAT(1X,d12.5,2X,I4,2X,5(D15.7),2x,F5.0)
	    If(errest.lt.13.d0) then
	      Print *,' Too much error in the above mode. '
	      Stop
	    Endif
            Write(18,45) sstart,om,L,p,p1p,difp,errest
 45     Format(1x,2(d12.5,2x),I4,2x,3(d12.5,2x),f5.0)
998	    i1 = 1
   	    write(19,*) om,l,i1,x2
	    write(19,*) p,p1p
	    P = Y(1)
	    P1P = Y(2)
            X1 = X2
	    DO 999 I=1,npts-1
	      X2 = X1 + DELX
	      If(dabs(xi).lt.1.d-12.and.dabs(om).lt.1.d-12.and.L.eq.0) then
	        p = 1.d0
	        p1p = 0.d0
                Go to 990
	      Endif
c  For ERN need to use analytic values if omega = 0.
	      If(dabs(om).lt.1.d-12.and.nab.eq.1.and.dabs(a2-1.d0).lt.1.d-12.and.dabs(b2-1.d0).lt.1.d-12) then
	        P = x2**L
	        p1p = L*x2**(L-1)
	        Go to 990
	      Endif
	      CALL ODEINT(Y,2,X1,X2,EPS,H1,HMIN,NOK,NBAD,DERIVS,BSSTEP)
	      P = Y(1)
	      P1P = Y(2)
c  Check to see if p1p is getting too large
	      If(p1p.gt.1.d299) then
	        Print*,' p1p is too large.  omega = ',om,' L = ',L,' p1p = ',p1p
	        Stop
	      Endif
 990	      i1 = i+1
	      write(19,*) om,l,i1,x2
	      write(19,*) p,p1p
 	      X1 = X2
 999        CONTINUE
1000    CONTINUE
	    WRITE(6,14) X1
  14	FORMAT(1X,'SF FROM EVOLVP =', 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,p,p1p,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)
        p = Y(1)
        p1p = 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) = p1p
c  I have scaled om so that om_new = r0*om_old/Sqrt(b2)
c  Of course r0 is scaled out of the equation.
c       DYDX(2) = -(2.d0/r + fp/(2.d0*f) + kp/(2.d0*k))*p1p
c    1               + (om**2/(k*f)
c    1             + xl*(xl+1.d0)/(k*r**2) + xi*ScalarR/k)*p
        DYDX(2) = -(2.d0/r + fp/(2.d0*f) + kp/(2.d0*k))*p1p
     1               + (om**2/(k*f)
     1             + xl*(xl+1.d0)/(k*r**2) + xi*ScalarR/k)*p
        RETURN
        END

