c This program does the integral over omega for the modes and computes the c final value for and the various components of c This version is for a general zero temperature black hole. Implicit None Integer i,is,j,k,n,nom,npts,ndouble,nab,nnab,junk Real*16 S1t,S2t,S4t,dif1,dif2,oml,RI,RR,a,itest,itexact,itexint Real*16 s,sin,dels,omi,om,omt,dom,delom,r, 1 f,fp,fpp,f3p,f4p,h,hp,hpp,h3p,mu,Lam,LLog, 1 C,C1,difS1,difS2,difS3,difS4,difS5,S1w,S2w,S3w,S4w,S5w, 3 ss,S1,S2,S3,S4,S5,m,Zeta, 4 Rtt,Rrr,Rthth,RicciS,S1om,S1m,S2om,S2m,S3om,S3m,S4om, 5 S4m,S5om,S5m,Ttt0,Trr0,Tthth0 Real*16 S1band,S2band,S3band,S4band,S5band,S1btest,S2btest,S3btest,S4btest,S5btest Real*16 a2,a3,b2,b3,xxi,aa2,aa3,bb2,bb3,AA22,AA33,BB22,BB33 Real*16 a4,a5,a6,a7,a8,b4,b5,b6,b7,b8 Real*16 xim1,xi,A22,A33,A34,A35,A44,B22,B33,B34,B35,B44 Real*16 Phi2,Ttt,Trr,Tthth,Phi2a,Ttta,Trra,Tththa,Tdifa, 1 Phi2w,Tttw,Trrw,Tththw,Tdifw, 1 Phi2aw,Tttaw,Trraw,Tththaw,Tdifaw,Phi2m,Tttm,Trrm,Tththm,Tdifm, 1 Tdif,Tdifs,Tdiftest,Tracea,Tracew,Tracem,Traceanomaly,tracedif, 1 traceadif,tracewdif Real*16 error,dom2,S1m2,S2m2,S3m2,S4m2,S5m2,dif1m,dif2m,dif3m,dif4m,dif5m Real*16 accS1m,accS2m,accS3m,accS4m,accS5m,accS1,accS2,accS3,accS4,accS5 Real*16 accS1mtot,accS2mtot,accS3mtot,accS4mtot,accS5mtot Real*16 S1Lacc,S2Lacc,S3Lacc,S4Lacc,S5Lacc,S1accom,S2accom,S3accom,S4accom,S5accom Real*16 Ttt1,Ttt2,Ttt3,Ttt4,Ttt5,Tttmtest Real*16 Trr1,Trr2,Trr3,Trr4,Trr5,Trrmtest Real*16 Tthth1,Tthth2,Tthth3,Tthth4,Tthth5,Tththmtest Real*16 Tdif1,Tdif2,Tdif3,Tdif4,Tdif5,Tdifmtest Real*16 Tttmabs,Trrmabs,Tththmabs,Tdifmabs Real*16 Tttmro,Trrmro,Tththmro,Tdifmro Real*16 Tttro,Trrro,Tththro,Tdifro Real*16 dif1r,dif2r,dif3r,dif4r,dif5r Real*16 accTtt1,accTtt2,accTtt3,accTtt4,accTtt5,accTtt Real*16 accTrr1,accTrr2,accTrr3,accTrr4,accTrr5,accTrr Real*16 accTthth1,accTthth2,accTthth3,accTthth4,accTthth5,accTthth Real*16 accTdif1,accTdif2,accTdif3,accTdif4,accTdif5,accTdif Real*16 accS1cut,accS2cut,accS3cut,accS4cut,accS5cut Real*16 Tttneed,Trrneed,Tththneed,Tdifneed Real*16 RicciSp,S4a,S4aw,S4tot Real*16 sm10,s10,s11,s30,s31,s50,s51,s52,s70,s71,s72,s73, 1 s90,s91,s92,s93,s111,s112,s113,s114, 2 s132,s133,s134,s135,s153,s154,s155, 3 s174,s175,s176, 4 s195,s196,s197,s216,s217,s237,s238,s258,s259 Real*16 ns10,ns30,ns50,ns51,ns70,ns71,ns72, 1 ns90,ns91,ns92,ns111,ns112,ns113, 2 ns132,ns133,ns134,ns153,ns154, 3 ns174,ns175, 4 ns195,ns196,ns216,ns237,ns258 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,y91,y92,y93, 1 y112,y113,y114,y133,y134, 2 y154,y155,y175,y176,y196,y217,y238 Real*16 w8,w8m,w8mn,w8ml12,w8pdw2,w8p2dw3 Dimension omi(30),delom(30),nom(30),om(30,8005),dom(30,8005),dom2(30,8005),Zeta(7) Dimension S1band(30),S2band(30),S3band(30),S4band(30),S5band(30) Dimension S1btest(30),S2btest(30),S3btest(30),S4btest(30),S5btest(30) Dimension S1Lacc(30),S2Lacc(30),S3Lacc(30),S4Lacc(30),S5Lacc(30) 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/metr/F,FP,FPP,F3P,F4P,H,HP,HPP,H3P Common/lln/xi,mu,Lam OPEN(11,FILE='s.dat',STATUS='OLD') OPEN(12,FILE='series.dat',STATUS='OLD') Open(14,file='omintz.dat',status='old') Open(15,file='Slsumz.out',status='old') c Open(10,file='Slsum_acc.out',status='old') Open(16,file='wkb8m0z.out',status='old') Open(17,file='Sn.out') Open(18,file='Tmn.out') Open(20,file='Sn.tst') Open(21,file='Tmn.tst') Open(22,file='Tmna.out') Open(23,file='Tmnw.out') Open(24,file='snm.tst') Open(33,file='Tmn_digits.out') Open(34,file='Tmn_roundoff.out') Open(35,file='Tmn_acc.out') Open(36,file='integral_acc.out') Open(37,file='integral_error.out') Open(38,file='cutoff_acc.out') Write(33,*) '0 Number of extra digits of accuracy that we need for the stress tensor' Write(33,*) '0 s Ttt Trr Tthth Tdif' Write(34,*) '0 Total roundoff error in computing the components of the stress tensor' Write(34,*) '0 s Ttt Trr Tthth Tdif' Write(35,*) '0 Total accuracy of the components of the stress tensor' Write(35,*) '0 s Ttt Trr Tthth Tdif' Write(36,*) '0 Accuracy of the contributions of the bands to the mode integrals' Write(36,*) '0 s band no. S1 S2 S3 S4 S5' c We will set mu = 1 and lam = 1. lam is the lower limit cutoff on the c omega integrals for the 1/om counterterms. C is Euler's constant. mu = 1.q0 lam = 1.q0 m = 0.q0 C = 0.57721566490153286q0 Zeta(3) = 1.20205690315959429q0 Zeta(5) = 1.03692775514336993q0 Zeta(7) = 1.00834927738192283q0 c sin is the value of s you want to start computing the stress tensor at Read(11,*) Junk Read(11,*) sin,dels,npts Read(12,*) Junk READ(12,*) xim1,nab If (qabs(xim1).gt.1.q10) then xi = 0.q0 Else xi = 1.q0/xim1 Endif Print*,'xi = ',xi Read(12,*) junk Read(12,*) a2,b2 If(nab.ge.2) then Read(12,*) junk Read(12,*) a3,b3 Endif If(nab.ge.3) then Read(12,*) junk Read(12,*) 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(18,3) xi,a2,b2,a3,b3,a4,b4 Write(18,*) '0 NAB A22 B22 A33 B33 A44 B44' Write(18,4) nab,A22,B22,A33,B33,A44,B44 Write(18,*) '0 s S4total Ttt Trr Tthth Tuu' 3 Format(5x,9(F9.5,3x)) 4 Format(5x,I3,7x,9(F9.5,3x)) Read(14,*) j Read(14,*) omi(1),delom(1),nom(1) If(j.gt.1) then Read(14,*) omi(2),delom(2),nom(2) Endif If(j.gt.2) then Read(14,*) omi(3),delom(3),nom(3) Endif If(j.gt.3) then Read(14,*) omi(4),delom(4),nom(4) Endif If(j.gt.4) then Read(14,*) omi(5),delom(5),nom(5) Endif If(j.gt.5) then Read(14,*) omi(6),delom(6),nom(6) Endif If(j.gt.6) then Read(14,*) omi(7),delom(7),nom(7) Endif If(j.gt.7) then Read(14,*) omi(8),delom(8),nom(8) Endif If(j.gt.8) then Read(14,*) omi(9),delom(9),nom(9) Endif If(j.gt.9) then Read(14,*) omi(10),delom(10),nom(10) Endif If(j.gt.10) then Read(14,*) omi(11),delom(11),nom(11) Endif If(j.gt.11) then Read(14,*) omi(12),delom(12),nom(12) Endif If(j.gt.12) then Read(14,*) omi(13),delom(13),nom(13) Endif If(j.gt.13) then Read(14,*) omi(14),delom(14),nom(14) Endif If(j.gt.14) then Read(14,*) omi(15),delom(15),nom(15) Endif If(j.gt.15) then Read(14,*) omi(16),delom(16),nom(16) Endif If(j.gt.16) then Read(14,*) omi(17),delom(17),nom(17) Endif If(j.gt.17) then Read(14,*) omi(18),delom(18),nom(18) Endif If(j.gt.18) then Read(14,*) omi(19),delom(19),nom(19) Endif If(j.gt.19) then Read(14,*) omi(20),delom(20),nom(20) Endif If(j.gt.20) then Read(14,*) omi(21),delom(21),nom(21) Endif If(j.gt.21) then Read(14,*) omi(22),delom(22),nom(22) Endif If(j.gt.22) then Read(14,*) omi(23),delom(23),nom(23) Endif If(j.gt.23) then Read(14,*) omi(24),delom(24),nom(24) Endif If(j.gt.24) then Read(14,*) omi(25),delom(25),nom(25) Endif If(j.gt.25) then Read(14,*) omi(26),delom(26),nom(26) Endif If(j.gt.26) then Read(14,*) omi(27),delom(27),nom(27) Endif If(j.gt.27) then Read(14,*) omi(28),delom(28),nom(28) Endif If(j.gt.28) then Read(14,*) omi(29),delom(29),nom(29) Endif If(j.gt.29) then Read(14,*) omi(30),delom(30),nom(30) Endif c Now read the values from wkb8m0z.out and make sure they agree Read(16,*) Junk Read(16,*) xxi,aa2,bb2,aa3,bb3 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(b2-bb2).gt.1.q-15.or.qabs(b3-bb3).gt.1.q-15) then print *,'Wrong values in wkb8m0z.out' print *, xi,a2,b2,a3,b3 print *, xxi,aa2,bb2,aa3,bb3 stop Endif Read(16,*) Junk Read(16,*) nnab,AA22,BB22,AA33,BB33 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) then print *,'Wrong values in wkb8m0z.out' print *, nab,A22,B22,A33,B33 print *, xxi,AA22,BB22,AA33,BB33 stop endif c First set up the values of om and dom to use. Do i=1,j om(i,1) = omi(i) Do n=2,nom(i) om(i,n) = om(i,n-1)+delom(i) Enddo Print *,'i = ',i,'om(i,1) = ',om(i,1),'om(i,nom) = ',om(i,nom(i)) DO n=1,nom(i) If(n.eq.1.or.n.eq.nom(i)) then dom(i,n) = delom(i)/3.q0 k = 0 else if(k.eq.0) then dom(i,n) = 4.q0*delom(i)/3.q0 k = 1 else dom(i,n) = 2.q0*delom(i)/3.q0 k = 0 endif Enddo c Next set up the values of om and dom to use when there are half the number of modes. DO n=1,nom(i),2 If(n.eq.1.or.n.eq.nom(i)) then dom2(i,n) = 2.q0*delom(i)/3.q0 k = 0 else if(k.eq.0) then dom2(i,n) = 4.q0*2.q0*delom(i)/3.q0 k = 1 else dom2(i,n) = 2.q0*2.q0*delom(i)/3.q0 k = 0 endif Enddo Enddo s = sin-dels Do is=1,npts s = s+dels Print *,s c Now compute the analytical approximation r = s+1.q0 Call Tmna(r,Ttta,Trra,Tththa) Tdifa = (Trra-Ttta)/f RicciS = fp**2/(2*f**2*h) - fpp/(f*h) + fp*hp/(2*f*h**2) + 2/r**2 - - 2/(h*r**2) - 2*fp/(f*h*r) + 2*hp/(h**2*r) RicciSp = -(f3p/(f*h)) - fp**3/(f**3*h) + (2*fp*fpp)/(f**2*h) - - (fp**2*hp)/(f**2*h**2) + (3*fpp*hp)/(2.*f*h**2) - - (fp*hp**2)/(f*h**3) + (fp*hpp)/(2.*f*h**2) - 4.q0/r**3 + 4.q0/(h*r**3) + - (2*fp)/(f*h*r**2) + (2*fp**2)/(f**2*h*r) - (2*fpp)/(f*h*r) + - (2*fp*hp)/(f*h**2*r) - (4*hp**2)/(h**3*r) + (2*hpp)/(h**2*r) Rthth = 1.q0/r**2 - 1/(r**2*h) - fp/(2*r*f*h) + hp/(2*r*h**2) Rrr = fp**2/(4*f**2*h) - fpp/(2*f*h) + fp*hp/(4*f*h**2) + hp/(h**2*r) Rtt = fp**2/(4*f**2*h) - fpp/(2*f*h) + fp*hp/(4*f*h**2) - fp/(f*h*r) c Normalize Phi2 by multiplying by 16*pi^2 LLog = qlog(mu**2*f/(4.q0*Lam**2)) Phi2a = - fp**2/(6.q0*f**2*h) - fp*hp/(12.q0*f*h**2) 1 + fpp/(6.q0*f*h) + fp/(3.q0*r*f*h) - (xi - 1.q0/6.q0)*RicciS*LLog S4a = f3p/(6.*f*h) + fp**3/(3.*f**3*h) - (fp*fpp)/(2.*f**2*h) + - (fp**2*hp)/(4.*f**2*h**2) - (fpp*hp)/(4.*f*h**2) + - (fp*hp**2)/(6.*f*h**3) - (fp*hpp)/(12.*f*h**2) - fp/(3.*f*h*r**2) - - fp**2/(3.*f**2*h*r) + fpp/(3.*f*h*r) - (fp*hp)/(3.*f*h**2*r) - - (xi-1.q0/6.q0)*(RicciSp*LLog + RicciS*fp/f) Print*,'s = ',s,' RicciS = ',Riccis,' Phi2a = ',Phi2a Write(21,*) 'Analytic approximation for Phi2 and Tab' Write(21,17)s,Phi2a,Ttta,Trra,Tththa,Tdifa Write(22,17)s,Ttta,Trra,Tththa,Tdifa 17 Format(1x,10(e14.7,1x)) c Now compute the wkb terms. c Read the wkb results for this s. 105 Read(16,*) 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(16,*) 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(16,*) ss, x30,x31,x50,x51,x70,x71,x72,x91,x92,x93,x112, 1 x113,x133,x134,x154,x155,x175,x196,x217 Read(16,*) 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 sm10 = - -0.040121333422050007q0*f/r**2 + 17*f*qlog(2.q0)/(240*r**2) - - f*qlog(4.q0)/(64*r**2) + 7*f*qlog(f)/(1920*r**2) - - 7*f*qlog(lam)/(960*r**2) - 7*f*qlog(r)/(960*r**2) s11 = - -0.072951000177433347q0 + 17*qlog(2.q0)/120 - qlog(4.q0)/32 + 7*qlog(f)/960 - - 7*qlog(lam)/480 - 7*qlog(r)/480 s10 = - -0.0656279164861975073q0 + 5*qlog(2.q0)/6 - qlog(4.q0)/4 - qlog(f)/24 + - qlog(lam)/12 + qlog(r)/12 ns10 = - 0.032829666755383341q0*f/r**2 - 17*f*qlog(2.q0)/(240*r**2) + - f*qlog(4.q0)/(64*r**2) - 7*f*qlog(f)/(1920*r**2) + - 7*f*qlog(lam)/(960*r**2) + 7*f*qlog(r)/(960*r**2) s31 = r**2/(12*f) s30 = - 2*C*r**2/f + 6*r**2*qlog(2.q0)/f - r**2*qlog(f)/f + - 2*r**2*qlog(lam)/f + 2*r**2*qlog(r)/f ns30 = - -0.1489612498195308406q0 + 5*qlog(2.q0)/6 - qlog(4.q0)/4 - qlog(f)/24 + - qlog(lam)/12 + qlog(r)/12 s52 = r**4/(18*f**2) s51 = - 2*r**4/(3*f**2) + 4*C*r**4/(3*f**2) + 4*r**4*qlog(2.q0)/f**2 - - 2*r**4*qlog(f)/(3*f**2) + 4*r**4*qlog(lam)/(3*f**2) + - 4*r**4*qlog(r)/(3*f**2) s50 = 28*r**4*Zeta(3)/(3*f**2) ns51 = r**2/(36*f) ns50 = - -2*r**2/(3*f) + 2*C*r**2/(3*f) + 2*r**2*qlog(2.q0)/(3*f) + - 2*r**2*qlog(4.q0)/(3*f) - r**2*qlog(f)/(3*f) + 2*r**2*qlog(lam)/(3*f) + - 2*r**2*qlog(r)/(3*f) s73 = 2*r**6/(45*f**3) s72 = - 4*r**6/(5*f**3) + 16*C*r**6/(15*f**3) + - 16*r**6*qlog(2.q0)/(5*f**3) - 8*r**6*qlog(f)/(15*f**3) + - 16*r**6*qlog(lam)/(15*f**3) + 16*r**6*qlog(r)/(15*f**3) s71 = 112*r**6*Zeta(3)/(15*f**3) s70 = 496*r**6*Zeta(5)/(15*f**3) ns72 = r**4/(90*f**2) ns71 = - -2*r**4/(15*f**2) + 4*C*r**4/(15*f**2) + - 4*r**4*qlog(2.q0)/(15*f**2) + 4*r**4*qlog(4.q0)/(15*f**2) - - 2*r**4*qlog(f)/(15*f**2) + 4*r**4*qlog(lam)/(15*f**2) + - 4*r**4*qlog(r)/(15*f**2) ns70 = 28*r**4*Zeta(3)/(15*f**2) s93 = - 88*r**8/(105*f**4) + 32*C*r**8/(35*f**4) + - 96*r**8*qlog(2.q0)/(35*f**4) - 16*r**8*qlog(f)/(35*f**4) + - 32*r**8*qlog(lam)/(35*f**4) + 32*r**8*qlog(r)/(35*f**4) s92 = 32*r**8*Zeta(3)/(5*f**4) s91 = 992*r**8*Zeta(5)/(35*f**4) s90 = 4064*r**8*Zeta(7)/(35*f**4) ns92 = - -4*r**6/(105*f**3) + 16*C*r**6/(105*f**3) + - 16*r**6*qlog(2.q0)/(105*f**3) + 16*r**6*qlog(4.q0)/(105*f**3) - - 8*r**6*qlog(f)/(105*f**3) + 16*r**6*qlog(lam)/(105*f**3) + - 16*r**6*qlog(r)/(105*f**3) ns91 = 16*r**6*Zeta(3)/(15*f**3) ns90 = 496*r**6*Zeta(5)/(105*f**3) s114 = - 160*r**10/(189*f**5) + 256*C*r**10/(315*f**5) + - 256*r**10*qlog(2.q0)/(105*f**5) - 128*r**10*qlog(f)/(315*f**5) + - 256*r**10*qlog(lam)/(315*f**5) + 256*r**10*qlog(r)/(315*f**5) s113 = 256*r**10*Zeta(3)/(45*f**5) s112 = 7936*r**10*Zeta(5)/(315*f**5) s111 = 32512*r**10*Zeta(7)/(315*f**5) ns113 = - -8*r**8/(945*f**4) + 32*C*r**8/(315*f**4) + - 32*r**8*qlog(2.q0)/(315*f**4) + 32*r**8*qlog(4.q0)/(315*f**4) - - 16*r**8*qlog(f)/(315*f**4) + 32*r**8*qlog(lam)/(315*f**4) + - 32*r**8*qlog(r)/(315*f**4) ns112 = 32*r**8*Zeta(3)/(45*f**4) ns111 = 992*r**8*Zeta(5)/(315*f**4) s135 = - 8768*r**12/(10395*f**6) + 512*C*r**12/(693*f**6) + - 512*r**12*qlog(2.q0)/(231*f**6) - 256*r**12*qlog(f)/(693*f**6) + - 512*r**12*qlog(lam)/(693*f**6) + 512*r**12*qlog(r)/(693*f**6) s134 = 512*r**12*Zeta(3)/(99*f**6) s133 = 15872*r**12*Zeta(5)/(693*f**6) s132 = 65024*r**12*Zeta(7)/(693*f**6) ns134 = - 32*r**10/(10395*f**5) + 256*C*r**10/(3465*f**5) + - 256*r**10*qlog(2.q0)/(3465*f**5) + 256*r**10*qlog(4.q0)/(3465*f**5) - - 128*r**10*qlog(f)/(3465*f**5) + 256*r**10*qlog(lam)/(3465*f**5) + - 256*r**10*qlog(r)/(3465*f**5) ns133 = 256*r**10*Zeta(3)/(495*f**5) ns132 = 7936*r**10*Zeta(5)/(3465*f**5) s155 = 2048*r**14*Zeta(3)/(429*f**7) s154 = 63488*r**14*Zeta(5)/(3003*f**7) s153 = 260096*r**14*Zeta(7)/(3003*f**7) ns154 = 512*r**12*Zeta(3)/(1287*f**6) ns153 = 15872*r**12*Zeta(5)/(9009*f**6) s176 = 28672*r**16*Zeta(3)/(6435*f**8) s175 = 126976*r**16*Zeta(5)/(6435*f**8) s174 = 520192*r**16*Zeta(7)/(6435*f**8) ns175 = 2048*r**14*Zeta(3)/(6435*f**7) ns174 = 63488*r**14*Zeta(5)/(45045*f**7) s197 = 458752*r**18*Zeta(3)/(109395*f**9) s196 = 2031616*r**18*Zeta(5)/(109395*f**9) s195 = 8323072*r**18*Zeta(7)/(109395*f**9) ns196 = 28672*r**16*Zeta(3)/(109395*f**8) ns195 = 126976*r**16*Zeta(5)/(109395*f**8) s217 = 4063232*r**20*Zeta(5)/(230945*f**10) s216 = 16646144*r**20*Zeta(7)/(230945*f**10) ns216 = 2031616*r**18*Zeta(5)/(2078505*f**9) s238 = 16252928*r**22*Zeta(5)/(969969*f**11) s237 = 66584576*r**22*Zeta(7)/(969969*f**11) ns237 = 4063232*r**20*Zeta(5)/(4849845*f**10) s259 = 32505856*r**24*Zeta(5)/(2028117*f**12) s258 = 133169152*r**24*Zeta(7)/(2028117*f**12) ns258 = 16252928*r**22*Zeta(5)/(22309287*f**11) Write(24,*) 'The snm sums. ' Write(24,*) s,sm10,s11,s10 Write(24,*) s31,s30,s52,s51,s50 Write(24,*) s73,s72,s71,s70,s93 Write(24,*) s92,s91,s90,s114,s113 Write(24,*) s112,s111,s135,s134,s133 Write(24,*) s132,s155,s154,s153,s176 Write(24,*) s175,s174,s197,s196,s195 Write(24,*) s217,s216,s238,s237,s259 Write(24,*) s258 Write(24,*)'Now for the nsnm sums. ' Write(24,*) ns10,ns30,ns51,ns50,ns72 Write(24,*) ns71,ns70,ns92,ns91,ns90 Write(24,*) ns113,ns112,ns111,ns134,ns133 Write(24,*) ns132,ns154,ns153,ns175,ns174 Write(24,*) ns196,ns195,ns216,ns237,ns258 w8 = sm10*um10+s10*u10+s30*u30+s31*u31+s50*u50+s51*u51+s52*u52+ - s70*u70+s71*u71+s72*u72+s91*u91+s92*u92+s93*u93+ - s112*u112+s113*u113+s114*u114+ - s133*u133+s134*u134+ - s154*u154+s155*u155+ - s175*u175+s176*u176+s196*u196+s217*u217+s238*u238 w8m = s10*w10+s30*w30+s50*w50+s51*w51+s70*w70+s71*w71+s72*w72+ - s90*w90+s91*w91+s92*w92+s111*w111+ 1 s112*w112+s113*w113+s132*w132+s133*w133+s134*w134+ 2 s153*w153+s154*w154+s174*w174+ - s175*w175+ 3 s195*w195+s196*w196+s216*w216+ 4 s237*w237+s258*w258 w8mn = ns10*w10+ns30*w30+ns50*w50+ns51*w51+ns70*w70+ns71*w71+ns72*w72+ - ns90*w90+ns91*w91+ns92*w92+ns111*w111+ 1 ns112*w112+ns113*w113+ns132*w132+ns133*w133+ns134*w134+ 2 ns153*w153+ns154*w154+ns174*w174+ - ns175*w175+ 3 ns195*w195+ns196*w196+ns216*w216+ 4 ns237*w237+ns258*w258 w8ml12 = s11*w10+s31*w30+s51*w50+s52*w51+s71*w70+s72*w71+s73*w72+ - s91*w90+s92*w91+s93*w92+s112*w111+ 1 s113*w112+s114*w113+s133*w132+s134*w133+s135*w134+ 2 s154*w153+s155*w154+s175*w174+ - s176*w175+ 3 s196*w195+s197*w196+s217*w216+s238*w237+s259*w258 w8pdw2 = s30*x30+s31*x31+s50*x50+s51*x51+s70*x70+s71*x71+s72*x72+ - s91*x91+s92*x92+s93*x93+s112*x112+s113*x113+ 1 s133*x133+s134*x134+ - s154*x154+s155*x155+ 2 s175*x175+s196*x196+s217*x217 w8p2dw3 = s50*y50+s51*y51+s52*y52+s70*y70+s71*y71+s72*y72+ - s91*y91+s92*y92+s93*y93+ - s112*y112+s113*y113+s114*y114+ - s133*y133+s134*y134+ - s154*y154+s155*y155+ - s175*y175+s176*y176+s196*y196+s217*y217+s238*y238 S1w = w8mn/(2.q0*r**2) S2w = -w8*h/(2.q0*r**2*f)+w8p2dw3/(8.q0*r**2)+w8pdw2/(2.q0*r**3)+w8m/(2.q0*r**4) S3w = w8ml12/(2.q0*r**2) S4w = -w8pdw2/(2.q0*r**2)-w8m/r**3 S5w = w8m/(2.q0*r**2) Write(20,17) s,S1w,S2w,S3w,S4w,S5w c Now compute the various components normalized in the traditional way I have been doing c with 16*pi^2 Phi2 and 90*8^4*pi^2*Tmn c I didn't put this in the common statement so I'll just assign it now c since it is only for one variable Phi2w = 4.q0*S5w Tttw = 92160.q0*((2.q0*xi+0.5q0)*S1w/f + (2.q0*xi-0.5q0)*S2w/h - + (2.q0*xi-0.5q0)*S3w/r**2 - xi*fp*S4w/(2.q0*f*h) - + ((2.q0*xi-0.5q0)*(-.25q0/r**2 + m**2 + xi*RicciS) - + xi*Rtt)*S5w) Trrw = 92160.q0*(-S1w/(2.q0*f) + S2w/(2.q0*h) - S3w/(2.q0*r**2) - +xi*(2.q0/r + fp/(2.q0*f))*S4w/h - +(1.q0/(8.q0*r**2) - m**2/2.q0 - xi*RicciS/2.q0 + xi*Rrr)*S5w) Tththw = 92160.q0*((2.q0*xi-0.5q0)*S1w/f + (2.q0*xi -0.5q0)*S2w/h - + 2.q0*xi*S3w/r**2 - xi*S4w/(r*h) - + (-xi/(2.q0*r**2) + (2.q0*xi-0.5q0)*(m**2+xi*RicciS) - + xi*Rthth)*S5w) c This may be more accurate than doing the subtractions. Tdifw = 92160.q0*(-(2.q0*xi+1.q0)*S1w/f + (1.q0-2.q0*xi)*S2w/h - 2.q0*xi*S3w/r**2 - + xi*(2.q0/r + fp/f)*S4w/h + (2.q0*xi/(4.q0*r**2) - 2.q0*xi*m**2 - 2.q0*xi**2*RicciS - + xi*(Rrr-Rtt))*S5w)/f c Let's test it Tdiftest = (Trrw-Tttw)/f error = (Tdifw-Tdiftest)/Tdifw Print *,'% difference in the 2 ways of computing Tdifw = ',error c This is the combination of the analytic approximation and the WKBfin approximation Phi2aw = Phi2a+Phi2w S4aw = S4a+S4w Tttaw = Ttta+Tttw Trraw = Trra+Trrw Tththaw = Tththa+Tththw Tdifaw = (Trra-Ttta)/f + Tdifw Write(21,*) 'WKBfin values of Phi2 and Tab' Write(21,17)s,Phi2w,Tttw,Trrw,Tththw,Tdifw Write(21,*) 'Analytic plus WKBfin values of Phi2 and Tab' Write(21,17)s,Phi2aw,Tttaw,Trraw,Tththaw,Tdifaw Write(23,17)s,Tttaw,Trraw,Tththaw,Tdifaw c Write(23,17)s,Phi2w,Tttw,Trrw,Tththw,Tdifw c This is to test that the trace is correct in the case that xi = 1/6 If(qabs(xi-1.q0/6.q0).lt.1.q-10) then tracew = Tttw+Trrw+2.q0*Tththw tracea = Ttta+Trra+2.q0*Tththa Call Traceanom(r,Traceanomaly) c The trace of the wkb terms should be zero so compare with one of the components tracewdif = tracew/Tttw traceadif = (Traceanomaly-Tracea)/Traceanomaly Print *,'for the analytic approximation tracea = ',tracea,' traceanomaly = ',traceanomaly,' %diff = ',traceadif Print *,'for wkb zero trace = ',tracew,' trace/Tttw = ',tracewdif Endif c Now do the integral over omega for the modes. S1m = 0.q0 S2m = 0.q0 S3m = 0.q0 S4m = 0.q0 S5m = 0.q0 S1m2 = 0.q0 S2m2 = 0.q0 S3m2 = 0.q0 S4m2 = 0.q0 S5m2 = 0.q0 accS1mtot = 32. accS2mtot = 32. accS3mtot = 32. accS4mtot = 32. accS5mtot = 32. itest = 0.q0 a = 0.005q0 Write(37,*) 'Accuracy from the wkb subtractions and L sums' Do i=1,j c What I want to do is to check the accuracy of each band and check the contribution c of each band to the entire sum S1band(i) = 0.q0 S2band(i) = 0.q0 S3band(i) = 0.q0 S4band(i) = 0.q0 S5band(i) = 0.q0 S1btest(i) = 0.q0 S2btest(i) = 0.q0 S3btest(i) = 0.q0 S4btest(i) = 0.q0 S5btest(i) = 0.q0 ndouble = 1 S2Lacc(i) = 32. S3Lacc(i) = 32. S4Lacc(i) = 32. S5Lacc(i) = 32. Do n=1,nom(i) 77 Read(15,*,end=99) ss,omt,S1om,S2om,S3om,S4om,S5om If(qabs(s-ss).gt.1.q-12.or.qabs(om(i,n)-omt).gt.1.q-12) Go to 77 c Now read in the errors from the sums over L and WKB subtractions c 78 Read(10,*,end=999) ss,omt,S2accom,S3accom,S4accom,S5accom c If(qabs(s-ss).gt.1.q-12.or.qabs(om(i,n)-omt).gt.1.q-12) Go to 78 c S2Lacc(i) = Qmin1(S2Lacc(i),S2accom) c S3Lacc(i) = Qmin1(S3Lacc(i),S3accom) c S4Lacc(i) = Qmin1(S4Lacc(i),S4accom) c S5Lacc(i) = Qmin1(S5Lacc(i),S5accom) c Now compute the integral over omega for the band S1band(i) = S1band(i)+S1om*dom(i,n) S2band(i) = S2band(i)+S2om*dom(i,n) S3band(i) = S3band(i)+S3om*dom(i,n) S4band(i) = S4band(i)+S4om*dom(i,n) S5band(i) = S5band(i)+S5om*dom(i,n) itest = itest + dom(i,n)/(om(i,n)**2+a**2)**3.5q0 c Now cut the number of modes in half and do it If(ndouble.eq.1) then S1btest(i) = S1btest(i)+S1om*dom2(i,n) S2btest(i) = S2btest(i)+S2om*dom2(i,n) S3btest(i) = S3btest(i)+S3om*dom2(i,n) S4btest(i) = S4btest(i)+S4om*dom2(i,n) S5btest(i) = S5btest(i)+S5om*dom2(i,n) ndouble = 0 Else ndouble = 1 Endif Enddo c This is because the sum over L is the same for S1 as S5 S1Lacc(i) = S5Lacc(i) Write(37,67) s,S1Lacc(i),S2Lacc(i),S3Lacc(i),S4Lacc(i),S5Lacc(i) S1m = S1m+S1band(i) S2m = S2m+S2band(i) S3m = S3m+S3band(i) S4m = S4m+S4band(i) S5m = S5m+S5band(i) Print *,'i = ',i,' itest = ',itest S1m2 = S1m2+S1btest(i) S2m2 = S2m2+S2btest(i) S3m2 = S3m2+S3btest(i) S4m2 = S4m2+S4btest(i) S5m2 = S5m2+S5btest(i) Enddo c Now compute the relative difference between the sums Write(37,*) 'Errors in individual bands.' Write(37,*) ' 1st line integral for half the modes. 2nd line integral for all the modes' Write(37,*) ' 3rd line relative difference ' Do i=1,j c First do the error from cutting the modes in half dif1m = qabs((S1band(i)-S1btest(i))/S1band(i)) dif2m = qabs((S2band(i)-S2btest(i))/S2band(i)) dif3m = qabs((S3band(i)-S3btest(i))/S3band(i)) dif4m = qabs((S4band(i)-S4btest(i))/S4band(i)) dif5m = qabs((S5band(i)-S5btest(i))/S5band(i)) Write(37,57) s,i,S1band(i),S2band(i),S3band(i),S4band(i),S5band(i) Write(37,57) s,i,S1btest(i),S2btest(i),S3btest(i),S4btest(i),S5btest(i) Write(37,57) s,i,dif1m,dif2m,dif3m,dif4m,dif5m c Now compute the accuracy of the integral for the band accS1m = -qlog10(dif1m) accS2m = -qlog10(dif2m) accS3m = -qlog10(dif3m) accS4m = -qlog10(dif4m) accS5m = -qlog10(dif5m) c Now compare the accuracy of the mode integral with the accuracy of the sums over L and wkb subtractions c The miniumum is again a good conservative estimate of the real accuracy c Otherwise you could have a very accurate mode integral but not enough values of L and fool yourself c In fact this is not the case - I really should be integrating each value of L separately c accS1m = qmin1(accS1m,S1Lacc(i)) c accS2m = qmin1(accS2m,S2Lacc(i)) c accS3m = qmin1(accS3m,S3Lacc(i)) c accS4m = qmin1(accS4m,S4Lacc(i)) c accS5m = qmin1(accS5m,S5Lacc(i)) c Write(37,*) 'Estimated accuracy including wkb subtraction' c Write(37,*) ' and sum over L and integral over omega for band number ',i Write(37,*) 'Estimated accuracy from the integral over omega for band number ',i Write(37,67) s,accS1m,accS2m,accS3m,accS4m,accS5m c Now do the relative contribution to the entire mode integral dif1r = qabs(S1band(i)/S1m) dif2r = qabs(S2band(i)/S2m) dif3r = qabs(S3band(i)/S3m) dif4r = qabs(S4band(i)/S4m) dif5r = qabs(S5band(i)/S5m) Write(37,*) 'Relative contribution of this band to the mode integral' Write(37,57) s,i,dif1r,dif2r,dif3r,dif4r,dif5r Write(38,57) s,i,dif1r,dif2r,dif3r,dif4r,dif5r c Estimate the accuracy of the cutoff If(i.eq.j) then accS1cut = -qlog10(dif1r) - 1.q0 accS2cut = -qlog10(dif2r) - 1.q0 accS3cut = -qlog10(dif3r) - 1.q0 accS4cut = -qlog10(dif4r) - 1.q0 accS5cut = -qlog10(dif5r) - 1.q0 Write(38,*) 'Estimate of the accuracy of the cutoff' Write(38,67) s,accS1cut,accS2cut,accS3cut,accS4cut,accS5cut Endif c Now determine the accuracy of each contribution to the mode integral c Note that we include loss of accuracy by roundoff error by subtracting the log(difnr) if difnr is greater than 1 accS1m = accS1m - qlog10(dif1r) accS2m = accS2m - qlog10(dif2r) accS3m = accS3m - qlog10(dif3r) accS4m = accS4m - qlog10(dif4r) accS5m = accS5m - qlog10(dif5r) Write(37,*) 'Estimated total accuracy of mode integration for band number ',i Write(37,67) s,accS1m,accS2m,accS3m,accS4m,accS5m Write(36,66) s,i,accS1m,accS2m,accS3m,accS4m,accS5m c Finally take the total accuracy of the mode integral to be the minium from all of the bands c This is new in this version but necessary to take into account the accuracy of the sums over L and wkb subtractions accS1mtot = qmin1(accS1mtot,accS1m) accS2mtot = qmin1(accS2mtot,accS2m) accS3mtot = qmin1(accS3mtot,accS3m) accS4mtot = qmin1(accS4mtot,accS4m) accS5mtot = qmin1(accS5mtot,accS5m) Enddo c Now do the total difference due to half the modes Write(37,*) 'Total difference due to half the modes. Same first three lines as above' Write(37,*) 'Fourth line is accuracy estimate' dif1m = (S1m-S1m2)/S1m dif2m = (S2m-S2m2)/S2m dif3m = (S3m-S3m2)/S3m dif4m = (S4m-S4m2)/S4m dif5m = (S5m-S5m2)/S5m accS1m = -qlog10(qabs(dif1m)) accS2m = -qlog10(qabs(dif2m)) accS3m = -qlog10(qabs(dif3m)) accS4m = -qlog10(qabs(dif4m)) accS5m = -qlog10(qabs(dif5m)) Write(37,59) s,S1m2,S2m2,S3m2,S4m2,S5m2 Write(37,59) s,S1m,S2m,S3m,S4m,S5m Write(37,59) s,dif1m,dif2m,dif3m,dif4m,dif5m Write(37,67) s,accS1m,accS2m,accS3m,accS4m,accS5m Write(37,*) 'Total estimated accuracy of the mode sums is from the minimum of each individual one' Write(37,*) 'It is given without an estimate of the cutoff in omega' Write(37,67) s,accS1mtot,accS2mtot,accS3mtot,accS4mtot,accS5mtot c Now for the total accuracy of the mode sums take the smallest of the accuracy due to the density c of modes and the accuracy of the cutoff. c c I will skip this now because the accuracy of the cutoff is not a good indicator of c what is going on. If I increase the range of omega in the last band it doesn't change c this estimate much and it should. c accS1mtot = qmin1(accS1mtot,accS1cut) c accS2mtot = qmin1(accS2mtot,accS2cut) c accS3mtot = qmin1(accS3mtot,accS3cut) c accS4mtot = qmin1(accS4mtot,accS4cut) c accS5mtot = qmin1(accS5mtot,accS5cut) 57 Format(1x,F12.5,2x,I3,2x,6(E14.7,2x)) 59 Format(1x,F12.5,7x,6(E14.7,2x)) 67 Format(1x,F12.5,6x,7(F12.5,4x)) 66 Format(1x,F12.5,2x,I3,2x,7(F12.5,2x)) 64 Format(1x,F12.5,2x,5(F12.1,2x)) oml = om(j,nom(j)) RR = qsqrt(oml**2+a**2) RI = qsqrt(omi(1)**2+a**2) itexact = 8.q0/(15.q0*a**6) 1 -(omi(1)/RI-2.q0*omi(1)**3/(3.q0*RI**3)+omi(1)**5/(5.q0*RI**5))/a**6 itexint = (oml/RR-2.q0*oml**3/(3.q0*RR**3)+oml**5/(5.q0*RR**5))/a**6 1 -(omi(1)/RI-2.q0*omi(1)**3/(3.q0*RI**3)+omi(1)**5/(5.q0*RI**5))/a**6 Print *,'itexint = ',itexint dif1 = (itexact-itexint)/itexact dif2 = (itexint-itest)/itexint Print *,'s = ',s,' cutoff error = ',dif1,' density error = ',dif2 Print *,'omi = ',omi(1),' oml1 = ',om(1,nom(1)) c Print *,'domi = ',dom(1,1),' domf = ',dom(1,nom(1)) c Print *,'domi1 = ',dom(1,2),' domf1 = ',dom(1,nom(1)-1) c Print *,'domi2 = ',dom(1,3),' domf2 = ',dom(1,nom(1)-2) If(j.gt.1) then Print *,'omi2 = ',omi(2),' oml2 = ',om(2,nom(2)) c Print *,'domi = ',dom(2,1),' domf = ',dom(2,nom(2)) c Print *,'domi1 = ',dom(2,2),' domf1 = ',dom(2,nom(2)-1) c Print *,'domi2 = ',dom(2,3),' domf2 = ',dom(2,nom(2)-2) Endif S1 = S1w+S1m S2 = S2w+S2m S3 = S3w+S3m S4 = S4w+S4m S5 = S5w+S5m c These are computed this way because they are canceling for small s. difS1 = qabs(S1/S1w) difS2 = qabs(S2/S2w) difS3 = qabs(S3/S3w) difS4 = qabs(S4/S4w) difS5 = qabs(S5/S5w) accS1 = accS1mtot + qmin1(0.q0,qlog10(difS1)) accS2 = accS2mtot + qmin1(0.q0,qlog10(difS2)) accS3 = accS3mtot + qmin1(0.q0,qlog10(difS3)) accS4 = accS4mtot + qmin1(0.q0,qlog10(difS4)) accS5 = accS5mtot + qmin1(0.q0,qlog10(difS5)) Write(20,17) s,S1m,S2m,S3m,S4m,S5m Write(20,17) s,S1,S2,S3,S4,S5 Write(20,17) s,difS1,difS2,difS3,difS4,difS5 Write(20,17) s,accS1,accS2,accS3,accS4,accS5 Write(37,*) 'Total accuracy of the Sums' Write(37,67) s,accS1,accS2,accS3,accS4,accS5 c These are the combinations they occur in. S1t = S1/f S2t = S2/h S4t = s*S4 Write(17,17) s,S1t,S2t,S3,S4t,S5 Phi2m = 4.q0*S5m Tttm = 92160.q0*((2.q0*xi+0.5q0)*S1m/f + (2.q0*xi-0.5q0)*S2m/h - + (2.q0*xi-0.5q0)*S3m/r**2 - xi*fp*S4m/(2.q0*f*h) - + ((2.q0*xi-0.5q0)*(-.25q0/r**2 + m**2 + xi*RicciS) - + xi*Rtt)*S5m) Trrm = 92160.q0*(-S1m/(2.q0*f) + S2m/(2.q0*h) - S3m/(2.q0*r**2) - +xi*(2.q0/r + fp/(2.q0*f))*S4m/h - +(1.q0/(8.q0*r**2) - m**2/2.q0 - xi*RicciS/2.q0 + xi*Rrr)*S5m) Tththm = 92160.q0*((2.q0*xi-0.5q0)*S1m/f + (2.q0*xi -0.5q0)*S2m/h - + 2.q0*xi*S3m/r**2 - xi*S4m/(r*h) - + (-xi/(2.q0*r**2) + (2.q0*xi-0.5q0)*(m**2+xi*RicciS) - + xi*Rthth)*S5m) c This may be more accurate than doing the subtractions. Tdifm = 92160.q0*(-(2.q0*xi+1.q0)*S1m/f + (1.q0-2.q0*xi)*S2m/h - 2.q0*xi*S3m/r**2 - + xi*(2.q0/r + fp/f)*S4m/h + (2.q0*xi/(4.q0*r**2) - 2.q0*xi*m**2 - 2.q0*xi**2*RicciS - + xi*(Rrr-Rtt))*S5m)/f c Let's test it Tdiftest = (Trrm-Tttm)/f error = (Tdifm-Tdiftest)/Tdifm Print *,'% difference in the 2 ways of computing Tdif = ',error c Now check the value of the trace for the modes Tracem = Tttm+Trrm+2.q0*Tththm c This is because the trace should be zero so compare it to one component Tracedif = Tracem/Tttm Print *,'For the modes the trace should be zero. It is = ',Tracem,' %diff = ',Tracedif c Now compute the terms in the sums and use them to get estimates of roundoff error Ttt1 = 92160.q0*(2.q0*xi+0.5q0)*S1m/f Ttt2 = 92160.q0*(2.q0*xi-0.5q0)*S2m/h Ttt3 = 92160.q0*(2.q0*xi-0.5q0)*S3m/r**2 Ttt4 = 92160.q0*(-xi*fp*S4m/(2.q0*f*h)) Ttt5 = 92160.q0*((2.q0*xi-0.5q0)*(-.25q0/r**2 + m**2 + xi*RicciS) - + xi*Rtt)*S5m Trr1 = 92160.q0*(-S1m/(2.q0*f)) Trr2 = 92160.q0*(S2m/(2.q0*h)) Trr3 = 92160.q0*(- S3m/(2.q0*r**2)) Trr4 = 92160.q0*xi*(2.q0/r + fp/(2.q0*f))*S4m/h Trr5 = 92160.q0*(1.q0/(8.q0*r**2) - m**2/2.q0 - xi*RicciS/2.q0 + xi*Rrr)*S5m Tthth1 = 92160.q0*(2.q0*xi-0.5q0)*S1m/f Tthth2 = 92160.q0*(2.q0*xi -0.5q0)*S2m/h Tthth3 = 92160.q0*2.q0*xi*S3m/r**2 Tthth4 = 92160.q0*(-xi*S4m/(r*h)) Tthth5 = 92160.q0*(-xi/(2.q0*r**2) + (2.q0*xi-0.5q0)*(m**2+xi*RicciS) - + xi*Rthth)*S5m Tdif1 = 92160.q0*(-(2.q0*xi+1.q0)*S1m/f) Tdif2 = 92160.q0*(1.q0-2.q0*xi)*S2m/h Tdif3 = 92160.q0*(-2.q0*xi*S3m/r**2) Tdif4 = 92160.q0*xi*(2.q0/r + fp/f)*S4m/h Tdif5 = 92160.q0*(2.q0*xi/(4.q0*r**2) - 2.q0*xi*m**2 - 2.q0*xi**2*RicciS - + xi*(Rrr-Rtt))*S5m c Check that I have the correct terms here c Tttmtest = (Tttm-Ttt1-Ttt2-Ttt3-Ttt4-Ttt5)/Tttm c Trrmtest = (Trrm-Trr1-Trr2-Trr3-Trr4-Trr5)/Trrm c Tththmtest = (Tththm-Tthth1-Tthth2-Tthth3-Tthth4-Tthth5)/Tththm c Tdifmtest = (Tdifm-(Tdif1+Tdif2+Tdif3+Tdif4+Tdif5)/f)/Tdifm c Print*,' Test of the Tab terms' c Print*, Tttmtest,Trrmtest,Tththmtest,Tdifmtest c Now compute the roundoff errors Tttmabs = qabs(Ttt1)+qabs(Ttt2)+qabs(Ttt3)+qabs(Ttt4)+qabs(Ttt5) Trrmabs = qabs(Trr1)+qabs(Trr2)+qabs(Trr3)+qabs(Trr4)+qabs(Trr5) Tththmabs = qabs(Tthth1)+qabs(Tthth2)+qabs(Tthth3)+qabs(Tthth4)+qabs(Tthth5) Tdifmabs = (qabs(Tdif1)+qabs(Tdif2)+qabs(Tdif3)+qabs(Tdif4)+qabs(Tdif5))/f Tttmro = -qlog10(qabs(Tttm)/Tttmabs) Trrmro = -qlog10(qabs(Trrm)/Trrmabs) Tththmro = -qlog10(qabs(Tththm)/Tththmabs) Tdifmro = -qlog10(qabs(Tdifm)/Tdifmabs) c Now compute the final values Phi2 = Phi2aw+Phi2m S4tot = S4aw+S4m Ttt = Tttaw + Tttm Trr = Trraw + Trrm Tthth = Tththaw + Tththm Tdif = Tdifaw + Tdifm Tdifs = s*Tdif Write(21,*) 'Mode contributions to Phi2,S4, and Tab' Write(21,17)s,Phi2m,S4m,Tttm,Trrm,Tththm,Tdifm Write(21,*) 'Final values of Phi2,S4, and Tab' Write(21,17)s,Phi2,S4tot,Ttt,Trr,Tthth,Tdif Write(18,17)s,Phi2,S4tot,Ttt,Trr,Tthth,Tdif,Tdifs Write(21,*) 'Roundoff error in Tabm' Write(21,67)s,Tttmro,Trrmro,Tththmro,Tdifmro c Now compute the roundoff error again. Tttro = -qlog10(qabs(Ttt)/(qabs(Tttaw) + Tttmabs)) Trrro = -qlog10(qabs(Trr)/(qabs(Trraw) + Trrmabs)) Tththro = -qlog10(qabs(Tthth)/(qabs(Tththaw) + Tththmabs)) Tdifro = -qlog10(qabs(Tdif)/(qabs(Tdifaw) + Tdifmabs)) Write(21,*) 'Roundoff error in Tab' Write(21,67)s,Tttro,Trrro,Tththro,Tdifro c Now add the two together to get the total roundoff error and print in a separate file Tttro = Tttro + Tttmro Trrro = Trrro + Trrmro Tththro = Tththro + Tththmro Tdifro = Tdifro + Tdifmro Write(34,64) s,Tttro,Trrro,Tththro,Tdifro c Now compute the accuracy of each of the contributions Ttt1, ..., Ttt5, etc. accTtt1 = accS1mtot - qmin1(0.q0,qlog10(qabs(Ttt1/Tttm))) accTtt2 = accS2mtot - qmin1(0.q0,qlog10(qabs(Ttt2/Tttm))) accTtt3 = accS3mtot - qmin1(0.q0,qlog10(qabs(Ttt3/Tttm))) accTtt4 = accS4mtot - qmin1(0.q0,qlog10(qabs(Ttt4/Tttm))) accTtt5 = accS5mtot - qmin1(0.q0,qlog10(qabs(Ttt5/Tttm))) Write(21,*) 'Accuracy of Ttt1 to Ttt5' Write(21,67) s,accTtt1,accTtt2,accTtt3,accTtt4,accTtt5 c Now compute the accuracy of Ttt, etc. by first taking the minimum accuracy and then subtracting c of the roundoff errors for Tttm and Ttt accTtt = qmin1(accTtt1,accTtt2) accTtt = qmin1(accTtt,accTtt3) accTtt = qmin1(accTtt,accTtt4) accTtt = qmin1(accTtt,accTtt5) c Note Tttro is now the total roundoff error for computing Ttt accTtt = qmax1(0.q0,accTtt - Tttro) Write(21,*) 'Total Accuracy of Ttt = ',accTtt accTrr1 = accS1mtot - qmin1(0.q0,qlog10(qabs(Trr1/Trrm))) accTrr2 = accS2mtot - qmin1(0.q0,qlog10(qabs(Trr2/Trrm))) accTrr3 = accS3mtot - qmin1(0.q0,qlog10(qabs(Trr3/Trrm))) accTrr4 = accS4mtot - qmin1(0.q0,qlog10(qabs(Trr4/Trrm))) accTrr5 = accS5mtot - qmin1(0.q0,qlog10(qabs(Trr5/Trrm))) Write(21,*) 'Accuracy of Trr1 to Trr5' Write(21,67) s,accTrr1,accTrr2,accTrr3,accTrr4,accTrr5 c Now compute the accuracy of Trr, etc. by first taking the minimum accuracy and then subtracting c of the roundoff errors for Trrm and Trr accTrr = qmin1(accTrr1,accTrr2) accTrr = qmin1(accTrr,accTrr3) accTrr = qmin1(accTrr,accTrr4) accTrr = qmin1(accTrr,accTrr5) c We need the max here because we can get an accuracy estimate that is less than c zero due to the subtraction. accTrr = qmax1(0.q0,accTrr - Trrro) Write(21,*) 'Total Accuracy of Trr = ',accTrr accTthth1 = accS1mtot - qmin1(0.q0,qlog10(qabs(Tthth1/Tththm))) accTthth2 = accS2mtot - qmin1(0.q0,qlog10(qabs(Tthth2/Tththm))) accTthth3 = accS3mtot - qmin1(0.q0,qlog10(qabs(Tthth3/Tththm))) accTthth4 = accS4mtot - qmin1(0.q0,qlog10(qabs(Tthth4/Tththm))) accTthth5 = accS5mtot - qmin1(0.q0,qlog10(qabs(Tthth5/Tththm))) Write(21,*) 'Accuracy of Tthth1 to Tthth5' Write(21,67) s,accTthth1,accTthth2,accTthth3,accTthth4,accTthth5 c Now compute the accuracy of Tthth, etc. by first taking the minimum accuracy and then subtracting c of the roundoff errors for Tththm and Tthth accTthth = qmin1(accTthth1,accTthth2) accTthth = qmin1(accTthth,accTthth3) accTthth = qmin1(accTthth,accTthth4) accTthth = qmin1(accTthth,accTthth5) accTthth = qmax1(0.q0,accTthth - Tththro) Write(21,*) 'Total Accuracy of Tthth = ',accTthth accTdif1 = accS1mtot - qmin1(0.q0,qlog10(qabs(Tdif1/(f*Tdifm)))) accTdif2 = accS2mtot - qmin1(0.q0,qlog10(qabs(Tdif2/(f*Tdifm)))) accTdif3 = accS3mtot - qmin1(0.q0,qlog10(qabs(Tdif3/(f*Tdifm)))) accTdif4 = accS4mtot - qmin1(0.q0,qlog10(qabs(Tdif4/(f*Tdifm)))) accTdif5 = accS5mtot - qmin1(0.q0,qlog10(qabs(Tdif5/(f*Tdifm)))) Write(21,*) 'Accuracy of Tdif1 to Tdif5' Write(21,67) s,accTdif1,accTdif2,accTdif3,accTdif4,accTdif5 c Now compute the accuracy of Tdif, etc. by first taking the minimum accuracy and then subtracting c of the roundoff errors for Tdifm and Tdif accTdif = qmin1(accTdif1,accTdif2) accTdif = qmin1(accTdif,accTdif3) accTdif = qmin1(accTdif,accTdif4) accTdif = qmin1(accTdif,accTdif5) accTdif = qmax1(0.q0,accTdif - Tdifro) Write(21,*) 'Total Accuracy of Tdif = ',accTdif Write(35,64) s,accTtt,accTrr,accTthth,accTdif c Now let's just see how many more digits we need Tttneed = qmax1(0.q0,3.q0-accTtt) Trrneed = qmax1(0.q0,3.q0-accTrr) Tththneed = qmax1(0.q0,3.q0-accTthth) Tdifneed = qmax1(0.q0,3.q0-accTdif) Write(33,64) s,Tttneed,Trrneed,Tththneed,Tdifneed Enddo c Print *,'RicciS = ',RicciS,' Rtt = ',Rtt c Print *,'Rrr = ',Rrr,' Rthth = ',Rthth Stop 99 Print *,'Read through Slsum.out',' s = ',s,' ss = ',ss,' i = ',i,' n = ',n Print *, ' om = ',om(i,n),' omt = ',omt Stop 999 Print *,'Read through Slsum_acc.out',' s = ',s,' ss = ',ss,' i = ',i,' n = ',n Print *, ' om = ',om(i,n),' omt = ',omt Stop End