c c This version uses double precision output from the programs c pmodesz.f and qmodesz.f. c This version is explicitly for m = 0 c This program has been modified to use continuous values of omega. c It is also in the form of pqewkb8.f instead of pqwkb8.f. c c It reads in the p and q modes and then normalizes the q modes using c the Wronskian. It then subtracts the WKB expression from various c products such as p*q. C NPTS MUST BE AT LEAST 2 FOR THE PROGRAM TO WORK C THE IMPORTANT OUTPUT IS PUT IN PQDIF.OUT IMPLICIT NONE Integer ifiles,nn,nl,nt,n,l,n1,j,ll,npts,ni,nf,li,lf,ntt,lt,jt,nmodes,nom,LI1,LF1,LI2,LF2 Integer Junk,nab,nnab Real*16 Q,QP,P,PP,PI,PPI,dels,m,qm,sp,sm,xkap,sin,eps,h1,hmin, 1 xim1,s,w,wcor,c,xn,pt,ppt,qt,qpt,mm,qq,sisi,st,omi,delom,omt,om Real*16 xi Real*16 xxi Real*16 r,f,fp,f2p,k,kp,Cerror,LogCe,sstart,sf Real*16 A22,A33,A34,A35,A44,B22,B33,B34,B35,B44 Real*16 AA22,AA33,AA34,AA35,AA44,BB22,BB33,BB34,BB35,BB44 Real*16 a2,a3,a4,a5,a6,a7,a8,b2,b3,b4,b5,b6,b7,b8 Real*16 aa2,aa3,bb2,bb3,aa4,bb4 PARAMETER(NN=8100,NL=15,NT=75) Dimension pt(nt),ppt(nt),qt(nt),qpt(nt) 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/EVQ/Q(NN,NL,NT),QP(NN,NL,NT) COMMON/EVP/P(NN,NL,NT),PP(NN,NL,NT) COMMON/DIFEQ/DELS,NPTS COMMON/PINT/PI(NN,NL),PPI(NN,NL) COMMON/PAR/M,xi,omi,delom,nom c This is where m is set to zero once and for all. m = 0.q0 OPEN(7,FILE='s.dat',STATUS='OLD') OPEN(17,FILE='series.dat',STATUS='OLD') open(8,file='pmodesz.a1',status='old') open(9,file='qmodesz.a1',status='old') OPEN(12,FILE='pqwz.dat',STATUS='OLD') OPEN(13,FILE='modes.dat',STATUS='OLD') open(16,file='roundoff.out') OPEN(19,FILE='pqz.out') open(20,file='pqtest.out') open(21,file='pqw8z.out') open(22,file='pqtestw8.out') open(23,file='pqtestw8m.out') open(24,file='pqtestw8pdw2.out') open(25,file='wkb8m0z.out',status='old') open(27,file='pqtestw8p2dw3.out') open(29,file='pqleq0.out') open(30,file='pqacc.out') Read(7,*) Junk Read(7,*) sin,dels,npts Read(13,*) Junk READ(13,*) omi,delom,nom,LI1,LF1 Read(12,*) Junk Read(12,*) ifiles,LI2,LF2 Read(17,*) Junk READ(17,*) xim1,nab If (qabs(xim1).gt.1.q10) then xi = 0.q0 Else xi = 1.q0/xim1 Endif Print*,'xi = ',xi If(ifiles.ge.2) then open(10,file='pmodesz.a2',status='old') open(11,file='qmodesz.a2',status='old') Endif Read(17,*) junk Read(17,*) a2,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(21,*) '0 xi a2 b2 a3 b3 a4 b4' Write(30,*) '0 xi a2 b2 a3 b3 a4 b4' Write(21,3) xi,a2,b2,a3,b3,a4,b4 Write(30,3) xi,a2,b2,a3,b3,a4,b4 Write(21,*) '0 NAB A22 B22 A33 B33 A44 B44' Write(30,*) '0 NAB A22 B22 A33 B33 A44 B44' Write(21,4) nab,A22,B22,A33,B33,A44,B44 Write(30,4) nab,A22,B22,A33,B33,A44,B44 c This is just to keep from having to change lsumz.f Write(21,*) '0' 3 Format(5x,9(F9.5,3x)) 4 Format(5x,I3,7x,9(F9.5,3x)) c Now read the values from the files and make sure they agree Read(25,*) Junk Read(25,*) xxi,aa2,bb2,aa3,bb3,aa4,bb4 c The problem is that the xi that is printed out is only to 3 digits so they don't agree. If(qabs(xi-xxi).gt.1.q-2.or.qabs(a2-aa2).gt.1.q-15.or. - qabs(a3-aa3).gt.1.q-15.or.qabs(a4-aa4).gt.1.q-15.or. - qabs(b2-bb2).gt.1.q-15.or.qabs(b3-bb3).gt.1.q-15.or.qabs(b4-bb4).gt.1.q-15) then print *,'Wrong values in wkb8m0z.out' print *, xi,a2,b2,a3,b3,a4,b4 print *, xxi,aa2,bb2,aa3,bb3,aa4,bb4 stop Endif Read(25,*) Junk Read(25,*) nnab,AA22,BB22,AA33,BB33,AA44,BB44 If((nab-nnab).ne.0.or.qabs(A22-AA22).gt.1.q-13.or. - qabs(B22-BB22).gt.1.q-13.or.qabs(A33-AA33).gt.1.q-13.or. - qabs(B33-BB33).gt.1.q-13.or.qabs(A44-AA44).gt.1.q-13.or. - qabs(B44-BB44).gt.1.q-13) then print *,'Wrong values in wkb8m0z.out' print *, nab,A22,B22,A33,B33,A44,B44 print *, xxi,AA22,BB22,AA33,BB33,AA44,BB44 stop endif Read(8,*) Junk Read(8,*) xxi,aa2,bb2,aa3,bb3,aa4,bb4 If(qabs(xi-xxi).gt.1.q-2.or.qabs(a2-aa2).gt.1.q-15.or. - qabs(a3-aa3).gt.1.q-15.or.qabs(a4-aa4).gt.1.q-15.or. - qabs(b2-bb2).gt.1.q-15.or.qabs(b3-bb3).gt.1.q-15.or.qabs(b4-bb4).gt.1.q-15) then print *,'Wrong values in pmodesz.a1' print *, xi,a2,b2,a3,b3,a4,b4 print *, xxi,aa2,bb2,aa3,bb3,aa4,bb4 stop Endif Read(8,*) Junk Read(8,*) nnab,AA22,BB22,AA33,BB33,AA44,BB44 If((nab-nnab).ne.0.or.qabs(A22-AA22).gt.1.q-13.or. - qabs(B22-BB22).gt.1.q-13.or.qabs(A33-AA33).gt.1.q-13.or. - qabs(B33-BB33).gt.1.q-13.or.qabs(A44-AA44).gt.1.q-13.or. - qabs(B44-BB44).gt.1.q-13) then print *,'Wrong values in pmodesz.a1' print *, nab,A22,B22,A33,B33,A44,B44 print *, xxi,AA22,BB22,AA33,BB33,AA44,BB44 stop endif Read(9,*) Junk Read(9,*) xxi,aa2,bb2,aa3,bb3,aa4,bb4 If(qabs(xi-xxi).gt.1.q-2.or.qabs(a2-aa2).gt.1.q-15.or. - qabs(a3-aa3).gt.1.q-15.or.qabs(a4-aa4).gt.1.q-15.or. - qabs(b2-bb2).gt.1.q-15.or.qabs(b3-bb3).gt.1.q-15.or.qabs(b4-bb4).gt.1.q-15) then print *,'Wrong values in qmodesz.a1' print *, xi,a2,b2,a3,b3,a4,b4 print *, xxi,aa2,bb2,aa3,bb3,aa4,bb4 stop Endif Read(9,*) Junk Read(9,*) nnab,AA22,BB22,AA33,BB33,AA44,BB44 If((nab-nnab).ne.0.or.qabs(A22-AA22).gt.1.q-13.or. - qabs(B22-BB22).gt.1.q-13.or.qabs(A33-AA33).gt.1.q-13.or. - qabs(B33-BB33).gt.1.q-13.or.qabs(A44-AA44).gt.1.q-13.or. - qabs(B44-BB44).gt.1.q-13) then print *,'Wrong values in qmodesz.a1' print *, nab,A22,B22,A33,B33 print *, xxi,AA22,BB22,AA33,BB33 stop endif If(ifiles.ge.2) then Print *,'THIS IS WRONG. I NEED TO CHANGE PQWKB8Z.F' Read(10,*) xxi,aa3,AA33,AA34,Bb2,Bb3,BB33,BB34 If(qabs(xi-xxi).gt.1.q-15.or.qabs(a3-aa3).gt.1.q-15.or. - qabs(A33-AA33).gt.1.q-15.or.qabs(A34-AA34).gt.1.q-15.or. - qabs(b2-bb2).gt.1.q-15.or.qabs(b3-bb3).gt.1.q-15.or. - qabs(B33-BB33).gt.1.q-15.or.qabs(B34-BB34).gt.1.q-15) then print *,'Wrong m,q and si values in pmodes.a2' print *, xi,a3,A33,A34,b2,b3,B33,B34 print *, xxi,aa3,AA33,AA34,Bb2,Bb3,BB33,BB34 stop endif Read(11,*) xxi,aa3,AA33,AA34,Bb2,Bb3,BB33,BB34 If(qabs(xi-xxi).gt.1.q-15.or.qabs(a3-aa3).gt.1.q-15.or. - qabs(A33-AA33).gt.1.q-15.or.qabs(A34-AA34).gt.1.q-15.or. - qabs(b2-bb2).gt.1.q-15.or.qabs(b3-bb3).gt.1.q-15.or. - qabs(B33-BB33).gt.1.q-15.or.qabs(B34-BB34).gt.1.q-15) then print *,'Wrong m,q and si values in qmodes.a2' print *, xi,a3,A33,A34,b2,b3,B33,B34 print *, xxi,aa3,AA33,AA34,Bb2,Bb3,BB33,BB34 stop endif endif c WRITE(19,1) xi,a3,A33,A34,b2,b3,B33,B34 c WRITE(20,1) xi,a3,A33,A34,b2,b3,B33,B34 c WRITE(21,1) xi,a3,A33,A34,b2,b3,B33,B34 c WRITE(22,1) xi,a3,A33,A34,b2,b3,B33,B34 c WRITE(23,1) xi,a3,A33,A34,b2,b3,B33,B34 c WRITE(24,1) xi,a3,A33,A34,b2,b3,B33,B34 c WRITE(27,1) xi,a3,A33,A34,b2,b3,B33,B34 1 FORMAT(1X,8(d12.5,2x)) DO 300 N=1,nom N1 = N + 1 om = omi + (n-1)*delom DO 100 L=LI1,LF1 LL = L + 1 - LI1 S = sin - DELS DO 50 J=1,NPTS S = S + DELS 12 Read(8,*) omt,lt,jt,st,pt(j),ppt(j) If(qabs(om-omt).gt.1.q-12.or.lt.ne.l.or.qabs(s-st).gt.1.q-12) Then c print *,'Wrong values read from pmodesa1.' c print *,om,omt,l,lt,j,jt,s,st Go to 12 endif c Allow the program to read through undesired modes. 14 Read(9,*) omt,lt,jt,st,qt(npts-j+1),qpt(npts-j+1) If(qabs(om-omt).gt.1.q-12.or.lt.ne.l.or.(npts-jt+1).ne.j) Then c print *,'Wrong values read from qmodesa1.' c print *,om,omt,l,lt,j,jt Go to 14 endif 50 Continue S = sin - DELS Do 100 J = 1,npts S = S + DELS c The p's and q's have the wrong normalization. The Wronskian c condition needs to be satisfied Call metric(s,f,fp,f2p,k,kp) r = s+1.q0 w = qt(j)*ppt(j) - pt(j)*qpt(j) wcor = 1.q0/(r**2*qsqrt(f*k)) c = WCOR/W c Print *,n,L,s,C c Now compute the roundoff error if there is one Cerror = qabs(w/(qt(j)*ppt(j))) LogCe = qlog10(Cerror) Write(16,88) s,om,L,Cerror,LogCe 88 Format(1x,2(e12.5,2x),I5,2x,2(e12.5,2x)) qt(j) = C*qt(j) qpt(j) = C*qpt(j) P(n1,ll,j) = pt(j) PP(n1,ll,j) = ppt(j) Q(N1,LL,J) = qt(j) QP(N1,LL,J) = qpt(j) 100 CONTINUE If(ifiles.lt.2) Go to 300 DO 200 L=LI2,LF2 LL = L + 1 - LI1 S = sin - DELS DO 150 J=1,NPTS S = S + DELS 22 Read(10,*) omt,lt,jt,st,pt(j),ppt(j) If(qabs(om-omt).gt.1.q-12.or.lt.ne.l.or.qabs(s-st).gt.1.q-12) Then c print *,'Wrong values read from pmodesa1.' c print *,om,omt,l,lt,j,jt,s,st Go to 22 endif c Allow the program to read through undesired modes. 24 Read(11,*) omt,lt,jt,st,qt(npts-j+1),qpt(npts-j+1) If(qabs(om-omt).gt.1.q-12.or.lt.ne.l.or.(npts-jt+1).ne.j) Then c print *,'Wrong values read from qmodesa1.' c print *,om,omt,l,lt,j,jt Go to 24 endif 150 Continue S = sin - DELS Do 200 J = 1,npts S = S + DELS c The q's have the wrong normalization. The Wronskian c condition needs to be satisfied. Call metric(s,f,fp,f2p,k,kp) r = s+1.q0 w = qt(j)*ppt(j) - pt(j)*qpt(j) wcor = 1.q0/(r**2*qsqrt(f*k)) c = WCOR/W c Now compute the roundoff error if there is one Cerror = qabs(w/(qt(j)*ppt(j))) LogCe = qlog10(Cerror) Write(16,88) s,om,L,Cerror,LogCe qt(j) = C*qt(j) qpt(j) = C*qpt(j) P(n1,ll,j) = pt(j) PP(n1,ll,j) = ppt(j) Q(N1,LL,J) = qt(j) QP(N1,LL,J) = qpt(j) 200 CONTINUE 300 CONTINUE li = li1 If(ifiles.eq.1) then lf = lf1 else lf = lf2 endif CALL PHISQ(sin,LI,LF) STOP END SUBROUTINE PHISQ(SINI,LI,LF) IMPLICIT NONE Integer nn,nl,nt,it,npts,ni,nf,li,lf,n,n1,lwrit,l,ll,nom,Lzero real*16 p,pp,q,qp,dels,eps,h1,hmin,sini,s,r,xn,om,omi,delom, 1 xl,pq,pqdif1,pqdif2,norm,ll1,tlp1,Tthth,ss, 2 phi2,term1,term2,term3,term4,term5, 3 term1d,term2d,term3d,term4d,term5d real*16 M,xi,pdqdr,qdpdr,dpdq real*16 w8,w8m,w8mll1,w8pdw2,w8p2dw3,pqdif3,pqdif4 real*16 m2,F,FP,F2P,F3P,F4P,F5P,F6P,H,HP,HPP,H3P, 1 H4P,H5P,h6p,V2,V2P,V2PP,V23P,V24P,V1,V1P,V1PP,V13P,V14P,v15p,X,XP, 2 xpp,X3P,X4P,Y,YP,YPP,Y3P,Y4P,v14p1,v14p2,v15p1,v15p2,v15p3 real*16 k,kp,diferror,pdqqdperr,pqacc,pdqqdpacc,dpdqacc real*16 RicciS,Rthth,dummy,pqdifh1,pqdifh2,pqdifh3,pqdifh4 real*16 phi2e,tththe,term1e,term2e,term3e,term4e,term5e,phi2d,tththd real*16 pqw8,pqw8m,pqw8mn,pqw8mll1,pqw8pdw2,pqw8p2dw3 real*16 w81,w82,w83,w84,w85,dif1,dif2,dif3,dif4,dif5,pqw810,pqw811,pqw812, 1 pqw813,pqw814,pqw815,pqtest,diftest,w8m1,w8m2,w8m3,w8m4,w8m5,w8pd1,w8pd2,w8pd3,w8pd4 Real*16 um10,u10,u30,u31,u50,u51,u52,u70,u71,u72,u91,u92,u93, 1 u112,u113,u114,u133,u134, 2 u154,u155,u175,u176,u196,u217,u238 Real*16 w10,w30,w50,w51,w70,w71,w72,w90,w91,w92,w111,w112,w113, 1 w132,w133,w134,w153,w154, 2 w174,w175,w195,w196, 3 w216,w237,w258 Real*16 x30,x31,x50,x51,x70,x71,x72,x91,x92,x93, 1 x112,x113,x133,x134, 2 x154,x155,x175,x196,x217 Real*16 y50,y51,y52,y70,y71,y72,y90,y91,y92,y93, 1 y110,y111,y112,y113,y114,y130,y131,y132,y133,y134,y150,y151,y152, 2 y153,y154,y155,y170,y171,y172,y173,y174,y175,y176,y196,y217,y238 Real*16 pqdif,pdqqdpdif,dpdqdif,pqwkb,pdqqdp,dpdqwkb,pdqqdpwkb Real*16 w8p1,w8p2,w8p3,w8p4,pqerror,logpqe PARAMETER(NN=8100,NL=15,NT=75) COMMON/EVP/P(NN,NL,NT),PP(NN,NL,NT) COMMON/EVQ/Q(NN,NL,NT),QP(NN,NL,NT) COMMON /DIFEQ/DELS,NPTS COMMON/WWK/w8,w8m,w8mll1,w8pdw2,w8p2dw3 COMMON/PAR/M,xi,omi,delom,nom COMMON/wkbl/R,F common/test/w81,w82,w83,w84,w85,w8m1,w8m2,w8m3,w8m4,w8m5,w8p1,w8p2,w8p3,w8p4, - w8pd1,w8pd2,w8pd3,w8pd4 Common/wkbterms/um10,u10,u30,u31,u50,u51,u52,u70,u71,u72,u91,u92,u93, 1 u112,u113,u114,u133,u134, 2 u154,u155,u175,u176,u196,u217,u238, 3 w10,w30,w50,w51,w70,w71,w72,w90,w91,w92,w111,w112,w113, 1 w132,w133,w134,w153,w154, 2 w174,w175,w195,w196, 3 w216,w237,w258, 3 x30,x31,x50,x51,x70,x71,x72,x91,x92,x93, 1 x112,x113,x133,x134, 2 x154,x155,x175,x196,x217, 3 y50,y51,y52,y70,y71,y72,y91,y92,y93, 1 y112,y113,y114,y133,y134, 2 y154,y155,y175,y176,y196,y217,y238 c This version is for massless fields 88 Format(1x,2(e12.5,2x),I5,2x,2(e12.5,2x)) Write(16,*) 'Now look at the error for pdqqpd' Write(30,*) '0 Number digits accuracy Number digits lost' Write(30,*) '0 s omega L pq pdqqdp dpdq subtraction error pdqqdp' m2 = 0.q0 S = SINI - DELS DO 500 IT=1,NPTS S = S + DELS R = S + 1.q0 Call metric(s,f,fp,f2p,k,kp) h = 1.q0/k WRITE(19,33) S 33 FORMAT(1X,'S = ',F10.7) c Read the wkb results for this s. 105 Read(25,*) ss,um10,u10,u30,u31,u50,u51,u52,u70,u71,u72,u91, 1 u92,u93,u112,u113,u114,u133,u134,u154,u155,u175,u176,u196, 2 u217,u238 Read(25,*) ss,w10,w30,w50,w51,w70,w71,w72,w90,w91,w92,w111, 1 w112,w113,w132,w133,w134,w153,w154,w174,w175,w195,w196, 2 w216,w237,w258 Read(25,*) ss, x30,x31,x50,x51,x70,x71,x72,x91,x92,x93,x112, 1 x113,x133,x134,x154,x155,x175,x196,x217 Read(25,*) ss,y50,y51,y52,y70,y71,y72,y91,y92,y93,y112,y113, 1 y114,y133,y134,y154,y155,y175,y176,y196,y217,y238 If(qabs(s-ss).gt.1.q-12) Then c Print *,'From wkb8, s = ',s,' ss = ',ss Go to 105 Endif DO 350 N=1,nom N1 = N + 1 om = omi + (n-1)*delom phi2e = 0.q0 term1e = 0.q0 term2e = 0.q0 term3e = 0.q0 term4e = 0.q0 term5e = 0.q0 phi2 = 0.q0 term1 = 0.q0 term2 = 0.q0 term3 = 0.q0 term4 = 0.q0 term5 = 0.q0 LWRIT = 0 Lzero = 0 DO 250 L=LI,LF LWRIT = LWRIT+1 LL = L+1-LI XL = L tlp1 = 2.q0*xl+1.q0 ll1 = xl*(xl+1.q0) pq = p(n1,ll,it)*q(n1,ll,it) pdqqdp = p(n1,ll,it)*qp(n1,ll,it) + pp(n1,ll,it)*q(n1,ll,it) c Now estimate any subtraction error that might occur. We'll use this towards the end of the program. diferror = qabs(pdqqdp/(p(n1,LL,it)*qp(n1,LL,it)-pp(n1,LL,it)*q(n1,LL,it))) dpdq = pp(n1,LL,it)*qp(n1,LL,it) c Now write the L=0 ones in a file If(L.eq.0) Write(29,*) s,om,pq,pdqqdp,dpdq CALL WKB(OM,XL) c Test whether I can predict the zeroth order wkb value of pq. c pqtest = (xl+0.5q0)*pq c diftest = (pqtest-1.q0/(2.q0*r*qsqrt(f)))/(1.q0/(2.q0*r*qsqrt(f))) c write(20,37) r,n,l,pqtest,diftest c 37 Format(1x,f6.3,2x,i4,2x,i4,2x,2(d12.5,2x)) c Test the convergence of the modes. This is done in sequence so that c you see the convergence betwee successive wkb approximations. c these are for w8. pqw810 = 1.q0/(2.q0*r**2*pq) pqw811 = pqw810-w81 pqw812 = pqw811 - w82 pqw813 = pqw812 - w83 pqw814 = pqw813 - w84 pqw815 = pqw814 - w85 dif1 = pqw811/w81 dif2 = pqw812/w82 dif3 = pqw813/w83 dif4 = pqw814/w84 dif5 = pqw815/w85 WRITE(22,49) om,L,pqw810,pqw811,pqw812,pqw813,pqw814,pqw815 WRITE(22,49) om,L,w81,w82,w83,w84,w85 WRITE(22,49) om,L,dif1,dif2,dif3,dif4,dif5 c These are for w8m. pqw810 = 2.q0*r**2*pq pqw811 = pqw810 - w8m1 pqw812 = pqw811 - w8m2 pqw813 = pqw812 - w8m3 pqw814 = pqw813 - w8m4 pqw815 = pqw814 - w8m5 dif1 = pqw811/w8m1 dif2 = pqw812/w8m2 dif3 = pqw813/w8m3 dif4 = pqw814/w8m4 dif5 = pqw815/w8m5 WRITE(23,49) om,L,pqw810,pqw811,pqw812,pqw813,pqw814,pqw815 WRITE(23,49) om,L,w8m1,w8m2,w8m3,w8m4,w8m5 WRITE(23,49) om,L,dif1,dif2,dif3,dif4,dif5 c These are for w8pdw2 pqw810 = -2.q0*r**2*pdqqdp-4.q0*r*pq pqw811 = pqw810 - w8p1 pqw812 = pqw811 - w8p2 pqw813 = pqw812 - w8p3 pqw814 = pqw813 - w8p4 dif1 = pqw811/w8p1 dif2 = pqw812/w8p2 dif3 = pqw813/w8p3 dif4 = pqw814/w8p4 WRITE(24,49) om,L,pqw810,pqw811,pqw812,pqw813,pqw814 WRITE(24,49) om,L,w8p1,w8p2,w8p3,w8p4 WRITE(24,49) om,L,dif1,dif2,dif3,dif4 c These are for w8p2dw3 pqw810 = 8.q0*r**2*dpdq+2.q0*h/(r**2*pq*f)+8.q0*r*pdqqdp 1 +8.q0*pq pqw811 = pqw810 - w8pd1 pqw812 = pqw811 - w8pd2 pqw813 = pqw812 - w8pd3 pqw814 = pqw813 - w8pd4 dif1 = pqw811/w8pd1 dif2 = pqw812/w8pd2 dif3 = pqw813/w8pd3 dif4 = pqw814/w8pd4 WRITE(27,49) om,L,pqw810,pqw811,pqw812,pqw813,pqw814 WRITE(27,49) om,L,w8pd1,w8pd2,w8pd3,w8pd4 WRITE(27,49) om,L,dif1,dif2,dif3,dif4 49 FORMAT(1X,D12.5,1X,I4,1X,7(D12.5,1X)) c w8pq = 1.q0/(2.q0*r**2*pq) c w8mpq = 2.q0*r**2*pq c w8pdw2pq = -2.q0*r**2*(pdqdr+qdpdr)-4.q0*r*pq c w8p2dw3pq = 8.q0*r**2*dpdq+2.q0*h/(r**2*pq*f)+8.q0*r*(qdpdr+pdqdr) c 1 +8.q0*pq pqwkb = w8m/(2.q0*r**2) pdqqdpwkb = -(w8pdw2/(2.q0*r**2)+w8m/r**3) dpdqwkb = -(w8*h/(2.q0*r**2*f)-w8p2dw3/(8.q0*r**2)-w8pdw2/(2.q0*r**3) 1 -w8m/(2.q0*r**4)) dif1 = (pq-w8m/(2.q0*r**2))/pq dif2 = (pdqqdp+w8pdw2/(2.q0*r**2)+w8m/r**3)/pdqqdp c This is because for xi = om = L = 0 I have dpdq = 0. If(qabs(xi).le.1.q-12.and.qabs(om).le.1.q-12.and.L.eq.0) then dif3 = 1.q0 Else dif3 = (dpdq+w8*h/(2.q0*r**2*f)-w8p2dw3/(8.q0*r**2)-w8pdw2/(2.q0*r**3) 1 -w8m/(2.q0*r**4))/dpdq Endif c Now we want the pq-wkb terms. c Note, there is a cancelation in the dominant part of the c pdqqdp term that makes it naturally less accurate than c the pq or dpdq terms. This is because of the form of p and c q near the horizon. c c What we are doing here is setting these to zero if the wkb c is canceling the modes to the accuracy they are computed. c Otherwise there are problems because the wkb is computed c to quadruple precision while the modes are not. Thus c pq-wkb has a random part if they cancel too well, which is c due to the modes and a nonrandom part due to wkb. The c nonrandom part will integrate to something other than c zero as the random part integrates to zero and this can c and does cause problems. c c We will also estimate the error here and print it into a file below. c If(Lzero.eq.1) then pqdif = 0.q0 pdqqdpdif = 0.q0 dpdqdif = 0.q0 pqacc = 0.q0 pdqqdpacc = 0.q0 dpdqacc = 0.q0 Else If(qabs(dif1).gt.3.q-1) then pqacc = 14.q0 Else pqacc = 14.q0 + qlog10(qabs(dif1)) Endif If(qabs(dif2).gt.3.q-1.and.diferror.gt.3.q-1) then pdqqdpacc = 14.q0 Else If(qabs(dif2).gt.3.q-1) then pdqqdpacc = 14.q0 + qlog10(diferror) Else If(diferror.gt.3.q-1) then pdqqdpacc = 14.q0 + qlog10(qabs(dif2)) Else pdqqdpacc = 14.q0 + qlog10(qabs(dif2)) + qlog10(diferror) Endif If(dif3.gt.3.q-1) then dpdqacc = 14.q0 Else dpdqacc = 14.q0 + qlog10(qabs(dif3)) Endif c This is so that higher values of L all get set to zero. If(pqacc.lt.1.q0.or.pdqqdpacc.lt.1.q0.or.dpdqacc.lt.1.q0) then pqdif = 0.q0 pdqqdpdif = 0.q0 dpdqdif = 0.q0 Lzero = 1 else pqdif = dif1*pq pdqqdpdif = dif2*pdqqdp c This is because dpdq = 0 in this case while the wkb approx. does not If(qabs(xi).le.1.q-12.and.qabs(om).le.1.q-12.and.L.eq.0) then dpdqdif = -dpdqwkb Else dpdqdif = dif3*dpdq Endif Endif Endif WRITE(21,*) S,om,L,pqdif,pdqqdpdif,dpdqdif c This is just so that I can see how much subtraction error there is for pdqqdp pdqqdperr = -qlog10(diferror) If(pdqqdperr.lt.0.3q0) pdqqdperr = 0.q0 Write(30,35) s,om,L,pqacc,pdqqdpacc,dpdqacc,pdqqdperr If(it.eq.1.or.it.eq.2.or.it.eq.npts) Then WRITE(19,29) s,om,L,DIF1,DIF2,dif3 29 FORMAT(1X,2(D12.5,1X),I4,1X,8(D12.5,1X)) 35 Format(1x,2(d12.5,1x),I4,2x,4(F5.1,2x)) c This is for the test. WRITE(20,29) s,om,L,pq,pdqqdp,dpdq WRITE(20,29) s,om,L,pqwkb,pdqqdpwkb,dpdqwkb WRITE(20,29) s,om,L,DIF1,DIF2,dif3 ENDIF 250 CONTINUE 350 continue 500 CONTINUE RETURN END