c This is set up for a type III singularity c c This version computes I1 in the same way as in Paper III and prints it out in the file ktrack.out c It does so using the same subtractions as are in that paper. c c This is a corrected version of Bradley's program to compute the c stress tensor for a radiation dominated universe. I've corrected c the mistakes in the starting values for the WKB approximation. c c I'm also changing the names of his input and output files to distinguish them from his. program BigRip4 implicit none external f, jac character junk*5, outfile*16, testfile*16, 1 ratfile*16 integer i, j, nk, npart, nt, tnum, u, otype, ktrck, mxstep, ktest real*16 k, delk, dk, m, pi, z, aTable, uTable, mu real*16 tstart, tstep, w, tfinal, tf, tol, h real*16 yrnsq, trnsq, wkbsq, tk, lm, kcut, old real*16 kmax dimension k(4200001), delk(4200000), nk(101), dk(100), old(100) dimension tk(100) common/cutoff/nk, npart, u common/param/m, pi, z, mu common/ks/k, ktest common/scale/w, tf common/sq/yrnsq, trnsq, wkbsq, tol common/ktr/ktrck common/h0/h common/o/old common/km/kmax c This is to print out the analytic approximation for m = 0 Open(10,file='rho.out') open(16,FILE='radbr.in') c open(16,FILE='/home/hickbb0/Summer/revisedprogram/typeIII/ c -radbr.in') open(18,FILE='radbr.ad.cfg') c open(18,FILE='/home/hickbb0/Summer/revisedprogram/typeIII/ c -adiabatic-reg/radbr.ad.cfg') read(18,*) junk, outfile read(18,*) junk, testfile read(18,*) junk, ratfile read(18,*) junk, mxstep read(18,*) junk, wkbsq read(18,*) junk, yrnsq read(18,*) junk, trnsq read(18,*) junk, ktrck read(18,*) junk, tol read(18,*) junk, otype open(15,FILE=trim(adjustl(outfile))) open(17,FILE='testfile') open(19,FILE='ktrack_AD.out') open(21,FILE=trim(adjustl(ratfile))) open(22,file='modecons.out') open(25,file='digiterror.out') open(80,file='psy2.dat') open(89,file='wkb.dat') open(101,file='scalevalues_AD.dat') write(80,334) 'time','k','rho','trace' 334 format(5a20) write(*,*) "OUTFILE=", outfile write(*,*) "TESTFILE=", testfile write(22,801) 'Time','Modes','check','SintError','SrintError' 801 Format(5a17) write(21,801) 'Time','Quantum','Analytic' c *** Initialization block. pi = 2.q0*qasin(1.q0) z = 0 write(*,*) 'line 55' read(16,*) junk, tstart, tstep read(16,*) junk, tfinal read(16,*) junk, tnum read(16,*) junk, kcut read(16,*) junk, lm read(16,*) junk, m read(16,*) junk, mu read(16,*) junk, z read(16,*) junk, w read(16,*) junk, h read(16,*) npart write(*,*) 'line 67' do i = 1, npart read(16,*) nk(i+1), dk(i), tk(i) end do k(1) = 0 nk(1) = 1 do i = 1, 50 old(i) = 0 end do Print *,'The cutoff lm = ',lm do i = 1, npart do j = nk(i)+i, nk(i+1)+i-1 delk(j-1) = dk(i) k(j) = k(j-1) + delk(j-1) If(i.eq.npart.and.j.eq.nk(i+1)+i-1) then kmax = k(j) Endif c Print*,i,j,' k(j) = ',k(j),' delk(j-1) = ',delk(j-1) end do c k(nk(i+1)+i) = k(nk(i+1)+i-1) end do ktest=j write(*,*) 'end k DO' write(15,*) "Start: ", tstart write(15,*) "Stop: ", tfinal write(15,*) "Step Size: ", tstep write(15,*) "Num Steps: ", tnum write(15,*) "Coupling Constant: ", z write(15,*) "Mass: ", m write(15,*) "W: ", w write(15,*) "H0: ", h write(15,801) 't','Quantum Rho','Quantum Trace', 'Classical','Scalar 1 Cuvature' write(17,*) "Start: ", tstart write(17,*) "Step Size: ", tstep write(17,*) "Num Steps: ", tnum write(17,*) "Coupling Constant: ", z write(17,*) "Mass: ", m write(17,*) "W: ", w write(17,*) "H0: ", h write(17,*) "Init done, calling phi_ODE" tf = tfinal write(*,*) 'call ODE' call phi_ODE(delk, tnum, tfinal, tstart, tstep, tk, otype, h, lm, 1 kcut, mxstep) stop end subroutine phi_ODE (delk, tnum, tfinal, tstart, tstep, tk, otype, 1 h, lm, kcut, mxstep) implicit none external f, jac character junk*4 integer neq, i, itol, itask, istate, iopt, lrw, iwork, 1 liw, mf, out_ios, nk, npart, tnum, tfinal, j, u, 2 otype, uold, mxstep, q, tcount, ktest real*16 y, t, tout, rtol, atol, rwork, error, k, 1 bessel1, bessel1dot, bessel2, bessel2dot, 2 delk, m, pi, z, e1, e2, tstart, tstep, 3 yint, lm, tk, h, wi, kcut, mu real*16 yr, yrdot, yi, yidot, 1 cons, tderive, tmp, a0, a1, sint, srint, a, san, 2 sran, cons2, srinterror, sinterror, san2, 3 san0, sran0, t_start,t_stop real*16 w0 real*16 yrnsq, trnsq, wkbsq, old,tol real*16 kcheck,fre,fpre,fim,fpim dimension nk(101), tk(100) dimension k(4200001), delk(4200000), iwork(20), wi(4200001), 1 rwork(45680100), y(40000000) common/cutoff/nk, npart, u common/param/m, pi, z, mu common/ks/k, ktest common/wis/wi common/error/sinterror,srinterror common/sq/yrnsq, trnsq, wkbsq, tol write(*,*) 'ODE' neq = 4*(nk(npart+1)+npart-1) c Write(*,*) neq mf = 10 itol = 1 read(18,*) junk, rtol write(17,*) junk, rtol read(18,*) junk, atol write(17,*) junk, atol itask = 1 istate = 1 iopt = 1 rwork(5) = 0.q0 rwork(6) = 0.q0 rwork(7) = 0.q0 rwork(8) = 0.q0 rwork(9) = 0.q0 rwork(10) = 0.q0 iwork(5) = 0 iwork(6) = mxstep iwork(7) = 0 iwork(8) = 0 iwork(9) = 0 iwork(10) = 0 lrw = 20 + 16*neq liw = 20 t = tstart tout = t y(1) = 0 y(2) = 1/2 y(3) = 1 y(4) = 0 write(17,*) "Initial Data (yr, yrdot, yidot):" u = 0 uold = 0 c This is from the radiation dominated version c a0 = qsqrt(tstart) c a1 = 0.5q0/qsqrt(tstart) c Print*,'a0 = ',a0 c Below is the original version do i=2,nk(npart+1)+npart-1 call wkb(tstart, k(i), y(4*i-3), y(4*i-2), y(4*i-1), y(4*i), 1 kcut) wi(i) = 0 end do c Let's test it c Print*,'Lowest k values of f, etc.',y(5),y(6),y(7),y(8) c write(19,*) 'This is to see what the first 4 are ',y(1),y(2),y(3),y(4) c write(19,333) tstart, k(1001),y(4001),y(4002),y(4003),y(4004) c do q=1,2101 c yr = y(4*q-3) c yrdot = y(4*q-2) c yi = y(4*q-1) c yidot = y(4*q) c write(*,*) m, a(tstart,0) c san0 = k(q)**2*(yr**2+yi**2 - a(tout,0)/(2.q0*k(q)) + c 1 (m**2)*a(tout,0)**3/4*k(q)**3) c san2 = k(q)**2*(yr**2+yi**2 - a(tout,0)/(2.q0*k(q))) c sran0 = k(q)**2*(yr**2+yi**2) c san0= sint(k(q), yr, yrdot, yi, yidot, tout,lm, 0) c sran0= srint(k(q), yr, yrdot, yi, yidot, tout, lm, 0) c c write(80,333) tout, k(q), yr**2+yi**2, san0, sran0, sran(tstart,lm,0) c 333 format(6e20.10) c end do do j=1,npart if (tout.ge.tk(j)) then u = j end if end do c print out integrand at the initial time trnsq = 1 call phi_int(tout, y, delk, lm, tfinal, otype, h, tstep) do i=1,tnum tout = tout + tstep do j=1,npart if (tout.ge.(tk(j)-tstep/2.q0)) then u = j end if end do if (u.ne.uold) then istate = 1 uold = u end if c write(*,*) tout c call cpu_time(t_start) call dlsode(f, neq, y, t, tout, itol, rtol, atol, itask, 1 istate, iopt, rwork, lrw, iwork, liw, jac, mf) c call cpu_time(t_stop) c write(*,*) 'solver time=',t_stop-t_start c This appears to follow the contribution to the integrand of one mode c in time c Make sure I don't print out sint from this call trnsq = 0 c do q=1001,1001 do q=2,2 yr = y(4*q-3) yrdot = y(4*q-2) yi = y(4*q-1) yidot = y(4*q) san0= sint(2,k(q), yr, yrdot, yi, yidot, tout,lm, 0) sran0= srint(2,k(q), yr, yrdot, yi, yidot, tout, lm, 0) c san0 = k(q)**2*(yr**2+yi**2 - a(tout,0)/(2.q0*k(q)) + c 1 (m**2)*a(tout,0)**3/4*k(q)**3) c san2 = k(q)**2*(yr**2+yi**2 - a(tout,0)/(2.q0*k(q))) c sran0 = k(q)**2*(yr**2+yi**2) c c write(19,333) tout, k(q),yr,yrdot,yi,yidot write(80,333) tout, k(q),yr,yrdot,yi,yidot,san0 c write(80,333) tout, qsqrt(tout), sran0, san0 333 format(7(e20.10,1x)) end do c Now set the quantity that tells you when to print out the integrand c for . If(i.eq.1.or.i.eq.300.or.i.eq.5000) then trnsq = 1 Else trnsq = 0 Endif c write (17,*) j, tout, tk(j) if (u.lt.npart) then do j=nk(u+1)+1+u-1,nk(npart+1)+npart-1 end do end if call phi_int(tout, y, delk, lm, tfinal, otype, h, tstep) end do return end subroutine phi_int (t, y, delk, lm, tfinal, otype, h, tstep) implicit none integer i, j, nk, npart, l, kend, otype, ktrck,q integer cerror real*16 k, delk, y, t, endpsum, psum, lm, yan, 1 val1, val2, val3, yint, m, pi, z, 2 sval1, sval2, sval3, sint, san, tsum, 3 endtsum, nlm, w, rho, tfinal, r, kt, 4 yrk, yik, yprk, ypik, h, c, c2, hot, wi, 5 propt, a, mu, srval1, srval2, srval3, 6 trsum, endtrsum, sran, srint, old, tstep, 7 cons, a0, a1, cons2, tderiv, cons3, tmp, constime real*16 rhoint,trint, tmp1,sinterror,srinterror,t_start,t_stop real*16 intrho, intrhotest,a3endtsum,a3endtrsum real*16 lasterror,s1,s2,s3 real*16 rhoGtt,rhom0,rhoI1,rhotot real*16 kmax,rhoUV,lam dimension nk(101),rhoint(1000),trint(1000) common/cutoff/nk, npart dimension k(4200001), delk(4200000), y(40000000), 1 psum(100), tsum(100), wi(4200001), trsum(100), old(100) common/param/m, pi, z, mu common/ks/k common/scale/w common/ktr/ktrck common/wis/wi common/o/old common/new/rhoint,trint common/error/sinterror, srinterror common/rhoa/rhoGtt,rhom0,rhoUV common/km/kmax c Set cerror to 0 for each new time and lam = kmax c lam is the actual upper limit cutoff used for k at a given time cerror = 0 lam = kmax c Make lasterror large initially so it won't kill the first value of k lasterror = 1.q1 endpsum = 0 endtsum = 0 endtrsum = 0 a0 = a(t,0) a1 = a(t,1) c Checking scale factor values write(101,990) t, a(t,0), a(t,1), a(t,2), a(t,3), a(t,4) 990 format(6es20.10) c write(17,110) t 110 format(//'t = ',f16.6) 121 format(ES22.11, ES22.11, ES22.11) 122 format(I1, ES22.11, ES22.11) 123 format(ES22.11, ES22.11) 124 format(ES15.5, ES15.5, ES15.5, ES15.5, ES15.5) 125 format(ES17.7, ES17.7, ES17.7, ES17.7,ES17.7,ES17.7,ES17.7,ES17.7) do j = 1, npart psum(j) = 0 tsum(j) = 0 trsum(j) = 0 do i=nk(j)+2+j-1, nk(j+1)+j-1, 2 c This is the original way, but I need to be able to test that not too c much error has occurred. So I need to do it a different way. c sval1 = sint(,i-2,k(i-2),y(4*(i-2)-3),y(4*(i-2)-2),y(4*(i-2)-1), c 1 y(4*(i-2)),t,lm,0)/3.q0 c sval2 = 4.q0*sint(i-1,k(i-1),y(4*(i-1)-3),y(4*(i-1)-2), c 1 y(4*(i-1)-1),y(4*(i-1)),t,lm,0)/3.q0 c sval3 = sint(i,k(i),y(4*(i)-3),y(4*(i)-2),y(4*(i)-1),y(4*(i)), c 1 t,lm,0)/3.q0 If(cerror.eq.1) then s1 = 0.q0 s2 = 0.q0 s3 = 0.q0 Else c The first number in sint and srint tells you where in the integration c interval that you are. s1 = sint(1,k(i-2),y(4*(i-2)-3),y(4*(i-2)-2),y(4*(i-2)-1), 1 y(4*(i-2)),t,lm,0) c I was having some problem with the error getting larger too soon when c things should still be OK so I'll try making sure it is small before c triggering the cutoff If(lasterror.lt.sinterror.and.lasterror.lt.1.q-14) then cerror = 1 Print *,'cerror has been set to 1' Print *,t,k(i),lasterror,sinterror Endif lasterror = sinterror s2 = sint(2,k(i-1),y(4*(i-1)-3),y(4*(i-1)-2), 1 y(4*(i-1)-1),y(4*(i-1)),t,lm,0) If(lasterror.lt.sinterror.and.lasterror.lt.1.q-10) then cerror = 1 Print *,'cerror has been set to 1' Print *,t,k(i),lasterror,sinterror Endif lasterror = sinterror s3 = sint(3,k(i),y(4*(i)-3),y(4*(i)-2),y(4*(i)-1),y(4*(i)), 1 t,lm,0) If(lasterror.lt.sinterror.and.lasterror.lt.1.q-14) then cerror = 1 Print *,'cerror has been set to 1' Print *,t,k(i),lasterror,sinterror Endif lasterror = sinterror c If cerror has turned to 1 then we shouldn't count this set of c contributions and we want to set the UV cutoff for this time If(cerror.eq.1) then s1 = 0.q0 s2 = 0.q0 s3 = 0.q0 lam = k(i-2) Print *,'t = ',t,' lam = ',lam Endif Endif sval1 = s1/3.q0 sval2 = 4.q0*s2/3.q0 sval3 = s3/3.q0 tsum(j) = tsum(j) + (sval1 + sval2 + sval3)*delk(i-2) c This was for checking the integration routine c If(qabs(t-6.02q0).lt.1.q-10) then c Write(17,*) 'This is for the sum over k for j = ',j,' i = ',i c Write(17,*) 'The time is ',t c Write(17,124) k(i-2),k(i-1),k(i),delk(i-2) c Write(17,124) s1,s2,s3,tsum(j) c Endif c Test how the integrand for a given mode changes in time c If(j.eq.1.and.i.eq.5) then c intrho = sval1*3.q0 c intrhotest = a0**3*intrho c Write(19,*) t,a0,j,i c Write(19,*) k(i-2),k(i-1),k(i),delk(i-2) c Write(19,190)k(i-2),t,a0,intrho,intrhotest c Write(19,190) y(4*(i-2)-3),y(4*(i-2)-2),y(4*(i-2)-1), c 1 y(4*(i-2)) c Endif If(cerror.eq.1) then srval1 = 0.q0 srval2 = 0.q0 srval3 = 0.q0 Else srval1 = srint(1,k(i-2),y(4*(i-2)-3),y(4*(i-2)-2),y(4*(i-2)-1), 1 y(4*(i-2)),t,lm,0)/3.q0 srval2 = 4.q0*srint(2,k(i-1),y(4*(i-1)-3),y(4*(i-1)-2), 1 y(4*(i-1)-1),y(4*(i-1)),t,lm,0)/3.q0 srval3 = srint(3,k(i),y(4*(i)-3),y(4*(i)-2),y(4*(i)-1),y(4*(i)), 1 t,lm,0)/3.q0 Endif trsum(j) = trsum(j) + (srval1 + srval2 + srval3)*delk(i-2) end do endtsum = endtsum + tsum(j) endtrsum = endtrsum + trsum(j) c write(17,124) j, tsum(j), endtsum, trsum(j), endtrsum end do c This is I1 for this version c Note that rhoUV is the contribution from modes above the UV cutoff to c rhos. endtsum = endtsum + san(t,lm,lam,0) endtrsum = endtrsum + sran(t,lm,0) a3endtsum = a0**3*endtsum a3endtrsum = a0**3*endtrsum Write(19,124) t,endtsum,a3endtsum,endtrsum,a3endtrsum c endtsum is the contribution of I1 to rho rhoI1 = endtsum rhotot = rhoI1+rhoGtt+rhom0+rhoUV Write(10,224) t,rhoI1,rhoUV,rhoGtt,rhom0,rhotot,lam 224 Format(1x,9(e14.7,1x)) if (old(2).eq.0) then cons = 0 else tderiv = (endtsum - old(2))/(2.q0*tstep) cons = tderiv - old(5) + old(7) tmp = qABS(tderiv) + qABS(old(5)) + qABS(old(7)) cons= cons/tmp end if old(1) = tderiv old(2) = old(3) old(3) = endtsum old(5) = endtrsum*a1/a0 old(7) = 4.q0*endtsum*a1/a0 if (old(9).eq.0) then cons3 = 0 else tderiv = (san(t,lm,kmax,0) - old(9))/(2.q0*tstep) cons3 = tderiv - old(12) + old(14) tmp1 = qABS(tderiv) + qABS(old(12)) + qABS(old(14)) cons3=cons3/tmp end if old(8) = tderiv old(9) = old(10) old(10) = san(t,lm,kmax,0) old(12) = sran(t,lm,0)*a1/a0 old(14) = 4.q0*san(t,lm,kmax,0)*a1/a0 rho = 3.q0/(8.q0*pi)*(a(t,1)/a(t,0))**2 hot = a(t,1)**2/a(t,0)**2 propt = t constime = t-tstep r =6*(a(t,2)/a(t,0) + a(t,1)**2/a(t,0)**2) 120 format(f7.4, f8.3, f10.3, f10.3) c write(17,121) san(t, lm, 0), endtsum, endtsum - san(t,lm,1) c write(17,121) sran(t, lm, 0), endtrsum, endtrsum - sran(t,lm,1) write(15,125) propt, endtsum, endtrsum, rho, r write(21,190) constime, cons, cons3 190 format(8es17.7) c end do return end subroutine f (neq, t, y, ydot) implicit none integer neq, i, nk, npart, j real*16 t, y, ydot, k, m, pi, z, a, a0, a1, a2, tk, mu dimension y(40000000), ydot(40000000), k(4200001), nk(101) common/cutoff/nk, npart, j common/param/m, pi, z, mu common/ks/k a0 = a(t,0) a1 = a(t,1) a2 = a(t,2) c write (17,*) t, j do i=1,nk(j+1)+j-1 ydot(4*i-3) = y(4*i-2) ydot(4*i-2) = -(k(i)**2/a0**2 + (6.q0*z-3.q0/2.q0)*(a2/a0 + a1**2/a0**2) 1 + 3.q0*a1**2/(4.q0*a0**2) + m**2) * y(4*i-3) ydot(4*i-1) = y(4*i) ydot(4*i) = -(k(i)**2/a0**2 + (6.q0*z-3.q0/2.q0)*(a2/a0 + a1**2/a0**2) 1 + 3.q0*a1**2/(4.q0*a0**2) + m**2) * y(4*i-1) c write(17,*) i end do do i=nk(j+1)+1+j-1,nk(npart+1)+npart-1 ydot(4*i-3) = 0 ydot(4*i-2) = 0 ydot(4*i-1) = 0 ydot(4*i) = 0 end do return end subroutine jac () return end real*16 function a (t, d) implicit none real*16 w, c, t, tf, h, n, ts, as, pi, m, z, mu integer d common/scale/w, tf common/h0/h common/param/m, pi, z, mu c c = -2.q0/(3.q0*w+1.q0) c n = 1.5q0 ts = 75.q0*qsqrt(3.q0)*qlog(10.q0)**(3.q0)/(8.q0*qsqrt(2.q0*pi)) as = 10.q0 c Type III Scale Factor with A=10 and B=5/4 if (d.eq.0) then a= as*QEXP((-(2.q0)*((2.q0*pi)**(1.q0/6.q0))*(-t ~ +ts)**(1.q0/3.q0))/(qsqrt(3.q0)*5.q0**(2.q0/3.q0))) end if if (d.eq.1) then a= 2.q0*as*QEXP((-(2.q0)*((2.q0*pi)**(1.q0/6.q0))*(-t ~+ts)**(1.q0/3.q0))/(qsqrt(3.q0)*5.q0**(2.q0/3.q0)))*(2.q0*pi) ~**(1.q0/6.q0)/(3.q0*qsqrt(3.q0)*5.q0**(2.q0/3.q0)*(-t+ts)**(2.q0/3.q0)) end if if (d.eq.2) then a= 4.q0*as*QEXP((-(2.q0)*((2.q0*pi)**(1.q0/6.q0))*(-t ~+ts)**(1.q0/3.q0))/(qsqrt(3.q0)*5.q0**(2.q0/3.q0)))*(5.q0*qsqrt(3.q0)* ~(2.q0*pi)**(1.q0/6.q0)+(10.q0*pi*(-t+ts))**(1.q0/3.q0))/ ~(135.q0*(25.q0*(-t+ts)**5)**(1.q0/3.q0)) end if if (d.eq.3) then a= 4.q0*as*QEXP((-(2.q0)*((2.q0*pi)**(1.q0/6.q0))*(-t ~+ts)**(1.q0/3.q0))/(qsqrt(3.q0)*5.q0**(2.q0/3.q0))) ~*(75.q0*qsqrt(3.q0)*5.q0**(1.q0/3.q0)* ~(2.q0*pi)**(1.q0/6.q0)+18.q0*(50.q0*pi*(-t+ts)) ~**(1.q0/3.q0)+2.q0*qsqrt(6.q0*pi)*(-t+ts) ~**(2.q0/3.q0))/(6075.q0*(-t+ts)**(8.q0/3.q0)) end if if (d.eq.4) then a = -16.q0*as*QEXP((-(2.q0)*((2.q0*pi)**(1.q0/6.q0))*(-t ~+ts)**(1.q0/3.q0))/(qsqrt(3.q0)*5.q0**(2.q0/3.q0)))*(2.q0*pi)**(1.q0/6.q0)* ~(5.q0**(1.q0/3.q0)*qsqrt(2.q0*pi)*(t-ts)-15.q0*(50.q0*qsqrt(3.q0)*5.q0** ~(1.q0/3.q0)+13.q0*5.q0**(2.q0/3.q0)*(2.q0*pi)**(1.q0/6.q0)*(-t+ts) ~**(1.q0/3.q0)+2.q0*qsqrt(3.q0)*(2.q0*pi)**(1.q0/3.q0)*(-t+ts)**(2.q0/3.q0)))/ ~(91125.q0*(-t+ts)**(11.q0/3.q0)) end if return end subroutine wkb (t, k, yr, yrdot, yi, yidot, kcut) implicit none real*16 k, t, yr, yrdot, yi, yidot, w0, w2, w4, w2p, m, 1 pi, z, kcut, tol, 2 yrnsq, trnsq, wkbsq, w, wp, w0p, mu common/param/m, pi, z, mu common/sq/yrnsq, trnsq, wkbsq, tol c call wkb3(t, k, w0, w0p, w2p, w4) call wkb3(t, k, w0, w0p, w2, w2p, w4) if (k.lt.kcut) then w = w0 wp = w0p c w = w2 c wp = w2p else c w = w2 c wp = w0p w = w4 wp = w2p end if 121 format(F22.11, ES22.11, ES22.11) if (wkbsq.eq.1) then c write(17, 121) k, w, wp end if yr = qsqrt(1.q0/(2.q0*w)) yrdot = -wp*(qsqrt(1.q0/(2.q0*w)))**3 yi = 0 yidot = -qsqrt(w/2.q0) write(89,79) k, yr, yrdot, yi, yidot 79 format(5e17.7) end c subroutine wkb3 (t, k, w0, w0p, w2p, w4) subroutine wkb3 (t, k, w0, w0p, w2, w2p, w4) implicit none real*16 k, t, w0, w2, w4, w2p, m, 1 pi, z, a, lm, a0, a1, a2, a3, a4, tol, 2 yrnsq, trnsq, wkbsq, w, wp, w0p, w0pp, w0ppp, 3 w0pppp, w2pp, mu, Hd, Hdd, test common/param/m, pi, z, mu a0 = a(t,0) a1 = a(t,1) a2 = a(t,2) a3 = a(t,3) a4 = a(t,4) Hd = 2.q0*a1*a2/a0**2 - 2.q0*a1**3/a0**3 Hdd = 2.q0*(-2.q0*a2*a1**2/a0**3 + (a2**2+a1*a3)/a0**2) 1 - 2.q0*(-3.q0*a1**4/a0**4+3.q0*a1**2*a2/a0**3) w0 = qsqrt(k**2/a0**2 + m**2) w0p = -((a1*k**2)/(a0**3*w0)) w0pp = (-3.q0*a1*w0p)/a0 + (a2*w0p)/a1 - w0p**2/w0 w0ppp = (3.q0*a1**2*w0p)/a0**2 - (3.q0*a2*w0p)/a0 - (a2**2*w0p)/a1**2 + 1 (a3*w0p)/a1 + w0p**3/w0**2 - (3.q0*a1*w0pp)/a0 + (a2*w0pp)/a1 - 2 (2.q0*w0p*w0pp)/w0 w0pppp = (-6.q0*a1**3*w0p)/a0**3 + (9.q0*a1*a2*w0p)/a0**2 + 1(2.q0*a2**3*w0p)/a1**3 - (3.q0*a3*w0p)/a0 - (3.q0*a2*a3*w0p)/a1**2 + (a4*w0p)/a1 - 2 (2.q0*w0p**4)/w0**3 + (6.q0*a1**2*w0pp)/a0**2 - (6.q0*a2*w0pp)/a0 - 3 (2.q0*a2**2*w0pp)/a1**2 + (2.q0*a3*w0pp)/a1 + (5.q0*w0p**2*w0pp)/w0**2 - 4 (2.q0*w0pp**2)/w0 - (3.q0*a1*w0ppp)/a0 + (a2*w0ppp)/a1 - 5 (2.q0*w0p*w0ppp)/w0 w2 = qsqrt((3.q0*a1**2)/(4.q0*a0**2) + w0**2 + ((3.q0*w0p**2)/(2.q0*w0**2) - 1 w0pp/w0)/2.q0 + (a1**2/a0**2 + a2/a0)*(-1.5q0 + 6.q0*z)) w2p = ((3.q0*Hd)/4.q0 + 2.q0*w0*w0p + ((w0p*w0pp)/w0**2 + (3.q0*((-2.q0*w0p**3)/w0**3 + 2 (2.q0*w0p*w0pp)/w0**2))/2.q0 - w0ppp/w0)/2.q0 + (-((a1*a2)/a0**2) + a3/a0 + 3 Hd)*(-1.5q0 + 6.q0*z))/(2.q0*w2) w2pp = -(w2p**2/w2) + ((3.q0*Hdd)/4.q0 + 2.q0*(w0p**2 + w0*w0pp) + ((-2.q0*w0p**2*w0pp)/w0**3 + 1 (w0p*w0ppp)/w0**2 + (w0pp**2 + w0p*w0ppp)/w0**2 + (3.q0*((-4.q0*w0p**2*w0pp)/w0**3 - 2 2.q0*((-3.q0*w0p**4)/w0**4 + (3.q0*w0p**2*w0pp)/w0**3) + (2.q0*(w0pp**2 + w0p*w0ppp))/w0**2))/2.q0 - 3 w0pppp/w0)/2.q0 + ((2.q0*a1**2*a2)/a0**3 - (a1*a3)/a0**2 - 4(a2**2 + a1*a3)/a0**2 + a4/a0 + Hdd)*(-1.5q0 + 6.q0*z))/(2.q0*w2) w4 = qsqrt((3.q0*a1**2)/(4.q0*a0**2) + w0**2 + ((3.q0*w2p**2)/(2.q0*w2**2) - 1 w2pp/w2)/2.q0 + (a1**2/a0**2 + a2/a0)*(-1.5q0 + 6.q0*z)) return end real*16 function sint (i,k, yr, yrdot, yi, yidot, t, 1 lm, mtest) implicit none integer i,mtest real*16 k, yr, yrdot, yi, yidot, t, lm, m, pi, z, a, mass, 1 a0, a1, a2, a3, a4, sun, srn, h1, wk2, 2 y2, yp2, yyp, tol, yrnsq, trnsq, wkbsq, 3 t1, t2, mu, sinterror,etest,kold real*16 sun1,sun2,sun3,sun4,suntot,yyp1,yyp2 common/param/m, pi, z, mu common/sq/yrnsq, trnsq, wkbsq, tol common y2,yp2,yyp common/error/sinterror,etest 150 format(F16.6, ES22.11, ES22.11) a0 = a(t,0) a1 = a(t,1) a2 = a(t,2) a3 = a(t,3) c a4 = a(t,4) if (mtest.eq.0) then m = m else m = 0 end if h1 = - 36.q0*(a3*a0**3 + 4.q0*a2*a1*a0**2+a1**3*a0)*a1*a0/a0**6 1 + 72.q0*(a2*a0**2+a1**2*a0)*a1**2*a0**2/a0**7 2 + 18.q0*(a2*a0**2+a1**2*a0)**2/a0**6 y2 = yr**2 + yi**2 yp2 = yrdot**2 + yidot**2 c Note this is for (fdot fstar + fdotstar f) but written in terms of c real and imaginary parts. yyp = 2.q0*(yr*yrdot + yi*yidot) c yyp1 = 2.q0*yr*yrdot c yyp2 = 2.q0* yi*yidot c I am going to comment this out and try the quantities I1 and I2 from c my paper III. if (k.eq.0.) then sint = 0.q0 c In this case I don't want to trigger the cutoff in the other program c There sinterro is initially set to 10. So here I make it large but c not that large sinterror = 9.q0 else wk2 = k**2 + m**2*a0**2 sun = 1.q0/(4.q0 * pi**2 * a0**4) * k**2*(a0**2*(yp2/a0 1 - a1*yyp/(2.q0*a0**2) + a1**2*y2/(4.q0*a0**3)) 2 + wk2*y2/a0 + (6.q0*z-1.q0)*(a1*(yyp - a1*y2/a0) 3 - (a1**2)*y2/a0)) c sun1 = 1.q0/(4.q0 * pi**2 * a0**4) * k**2*a0**2*yp2/a0 c sun2 = 1.q0/(4.q0 * pi**2 * a0**4) * k**2*(a0**2*( c 1 - a1*yyp/(2.q0*a0**2) )) c sun3 = 1.q0/(4.q0 * pi**2 * a0**4) * k**2*(a0**2*( c 1 a1**2*y2/(4.q0*a0**3))) c sun4 = 1.q0/(4.q0 * pi**2 * a0**4) * k**2*( c 2 + wk2*y2/a0) c suntot = sun1+sun2+sun3+sun4 c This if statement does not work right - because it gives the wrong c subtraction term when k = lm. It needs to be more sophisticated c because it depends on which integration interval you are on - the c one below lm or the one just above it. Both contain k = lm c if (k.lt.lm) then If(qabs(k-lm).lt.1.q-10) then If(i.eq.2) then Print*,'lm is in the middle of an integration interval' Print*,'k = ',k,' lm = ',lm Stop Else if(i.eq.3) then srn = 1.q0/(4.q0 * pi**2 * a0**4) * k**2 * 1 (k + 1.q0/k*(m**2*a0**2/2.q0 - (3.q0*z-1.q0/2.q0)* 2 (a1**2))) Else if(i.eq.1) then srn = 1.q0/(4.q0 * pi**2 * a0**4) * k**2 * 1 (k + 1.q0/k*(m**2*a0**2/2.q0 - (3.q0*z-1.q0/2.q0)* 2 (a1**2)) + (1.q0/k**3)*(-m**4*a0**4/8.q0 3 - (z-1.q0/6.q0)*(3.q0*m**2*a0**2*a1**2/2.q0) + 4 (z-1.q0/6.q0)**2*h1*a0**4/4.q0)) c This last case should never happen Else Print *,'From sint, i = ',i Stop Endif Else if(k.lt.lm) then srn = 1.q0/(4.q0 * pi**2 * a0**4) * k**2 * 1 (k + 1.q0/k*(m**2*a0**2/2.q0 - (3.q0*z-1.q0/2.q0)* 2 (a1**2))) Else srn = 1.q0/(4.q0 * pi**2 * a0**4) * k**2 * 1 (k + 1.q0/k*(m**2*a0**2/2.q0 - (3.q0*z-1.q0/2.q0)* 2 (a1**2)) + (1.q0/k**3)*(-m**4*a0**4/8.q0 3 - (z-1.q0/6.q0)*(3.q0*m**2*a0**2*a1**2/2.q0) + 4 (z-1.q0/6.q0)**2*h1*a0**4/4.q0)) End if c This is the exact subtraction term for I1. srn = (1.q0/(4.q0*pi**2*a0**4))*k**2*qsqrt(wk2) sinterror=qabs(sun-srn)/(qabs(sun)+qabs(srn)) c The following was not working to get rid of terms with too much c roundoff error. I use cerror in the subroutine that sums over the c sint contributions to do that. c Then this subroutine never gets called again once the roundoff error c is too large c c if(sinterror.lt.tol) then c sint = 0 c This is just a flag to say that sint has been set to zero c sinterror = -1.q0 c else c Here we are just doing the same subtraction for all values of k. c This gives I1 if we set the analytic approximation to zero as we do c below sint = sun - srn c sinterror = -qlog10(qabs(sint)/(qabs(sun)+qabs(srn))) c sinterror = -qlog10(sinterror) c end if c sinterror=qabs(sun-srn)/(qabs(sun)+qabs(srn)) end if if (trnsq.eq.1) then c write(17,35) t, k, sint, qabs((sun-srn)/srn) c Write(17,35) t,a0,a1,a2,a3 c Write(17,35),t, k, yr, yrdot,yi,yidot c Write(17,35),t, k, y2,yp2,yyp,yyp1,yyp2 c write(17,35) t, k, sun1,sun2,sun3,sun4,suntot,sun write(17,35) t, k, sint, sinterror end if c 35 Format(1X,10(e12.5,1X)) 35 Format(1X,10(e17.10,1X)) return end real*16 function san(t, lm,lam, mtest) implicit none integer mtest real*16 t, lm, m, pi, z, a, a0, a1, a2, a3, a4, r, u, 1 h1, h3, w, tf, c, l, gm, mass, mu real*16 rhoGtt,rhom0,rhoUV real*16 lam,omlam,kmax c Note lm is the cutoff between where the subtraction is to order k^0 c and where it is order 1/k in the anderson/eaker scheme - which we c are not using in this version of the program - so it is irrelevant. c lam is the largest value of k that is used for a given time t. This c changes with t. common/param/m, pi, z, mu common/scale/w common/rhoa/rhoGtt,rhom0,rhoUV a0 = a(t,0) a1 = a(t,1) a2 = a(t,2) a3 = a(t,3) a4 = a(t,4) h1 = - 36.q0*(a3*a0**3 + 4.q0*a2*a1*a0**2+a1**3*a0)*a1*a0/a0**6 1 + 72.q0*(a2*a0**2+a1**2*a0)*(a1*a0)**2/a0**7 2 + 18.q0*(a2*a0**2+a1**2*a0)**2/a0**6 h3 = 3.q0*a1**4/a0**4 c = -2.q0/(3.q0*w + 1.q0) if (mtest.eq.1) then u = mu else u = m end if if (lm.gt.0) then l = qlog(u**2*a0**2/(4.q0*lm**2)) else l = 0 end if if (mtest.eq.0) then m = m gm = -3.q0*a1**2*m**2/a0**2 san = 1.q0/(2880.q0*pi**2)*((-1.q0/6.q0)*h1 + h3) 1 + gm/(288.q0*pi**2) - m**4*(1.q0/2.q0+l)/(64.q0*pi**2) 2 + (z-1.q0/6.q0)*(h1/(288.q0*pi**2) + gm*(3.q0+l)/(16.q0*pi**2)) 3 + (z-1.q0/6.q0)**2*(h1*(2.q0+l)/(32.q0*pi**2) 4 - 9.q0/(4.q0*pi**2)*(a1**4/a0**4 + a1**2*a2/a0**3)) else m = 0 gm = -3.q0*a1**2*m**2/a0**2 san = 1.q0/(2880.q0*pi**2)*((-1.q0/6.q0)*h1 + h3) 1 + gm/(288.q0*pi**2) - m**4*(1.q0/2.q0+l)/(64.q0*pi**2) 2 + (z-1.q0/6.q0)*(h1/(288.q0*pi**2) + gm*(3.q0+l)/(16.q0*pi**2)) 3 + (z-1.q0/6.q0)**2*(h1*(2.q0+l)/(32.q0*pi**2) 4 - 9.q0/(4.q0*pi**2)*(a1**4/a0**4 + a1**2*a2/a0**3)) end if c I want the part propotional to G_t^t rhoGtt = gm/(288.q0*pi**2) c I want the analytic approximation with m = 0 rhom0 = 1.q0/(2880.q0*pi**2)*((-1.q0/6.q0)*h1 + h3) 35 Format(1X,10(e17.10,1X)) c There is no analytic approximation for the adiabatic term we are c subtracting here since we are just concerned with I1. san = 0.q0 c Now I want to estimate the contribution from modes above the cutoff c I do this assuming it is a fourth order adiabatic state - i.e. that c kcut is smaller than the largest values of k used. omlam = qsqrt(lam**2+m**2*a0**2) rhoUV = - a1**4/(160.q0*a0**4*Pi**2) + (a1**2*a2)/(120.q0*a0**3*Pi**2) + - (a1**2*m**2)/(96.q0*a0**2*Pi**2) - - (35*a0**4*a1**4*lam*m**8)/(1536.q0*omlam**9*Pi**2) - - (a1**4*lam**7)/(60.q0*a0**4*omlam**7*Pi**2) - - (a1**2*a2*lam**7)/(120.q0*a0**3*omlam**7*Pi**2) - - (7*a1**4*lam**5*m**2)/(120.q0*a0**2*omlam**7*Pi**2) - - (7*a1**2*a2*lam**5*m**2)/(240.q0*a0*omlam**7*Pi**2) - - (7*a1**4*lam**3*m**4)/(96.q0*omlam**7*Pi**2) - - (7*a0*a1**2*a2*lam**3*m**4)/(192.q0*omlam**7*Pi**2) + - (5*a0**2*a1**4*lam*m**6)/(1536.q0*omlam**7*Pi**2) + - (a1**4*lam**5)/(240.q0*a0**4*omlam**5*Pi**2) + - (a1**2*a2*lam**5)/(96.q0*a0**3*omlam**5*Pi**2) - - (a2**2*lam**5)/(960.q0*a0**2*omlam**5*Pi**2) + - (a1*a3*lam**5)/(480.q0*a0**2*omlam**5*Pi**2) + - (a1**4*lam**3*m**2)/(96.q0*a0**2*omlam**5*Pi**2) + - (5*a1**2*a2*lam**3*m**2)/(192.q0*a0*omlam**5*Pi**2) - - (a2**2*lam**3*m**2)/(384.q0*omlam**5*Pi**2) + - (a1*a3*lam**3*m**2)/(192.q0*omlam**5*Pi**2) + - (a1**4*lam*m**4)/(256.q0*omlam**5*Pi**2) - - (a1**4*lam**4)/(240.q0*a0**4*omlam**4*Pi**2) - - (a1**2*a2*lam**4)/(96.q0*a0**3*omlam**4*Pi**2) + - (a2**2*lam**4)/(960.q0*a0**2*omlam**4*Pi**2) - - (a1*a3*lam**4)/(480.q0*a0**2*omlam**4*Pi**2) - - (a1**4*lam**2*m**2)/(120.q0*a0**2*omlam**4*Pi**2) - - (a1**2*a2*lam**2*m**2)/(48.q0*a0*omlam**4*Pi**2) + - (a2**2*lam**2*m**2)/(480.q0*omlam**4*Pi**2) - - (a1*a3*lam**2*m**2)/(240.q0*omlam**4*Pi**2) - - (a1**4*m**4)/(240.q0*omlam**4*Pi**2) - - (a0*a1**2*a2*m**4)/(96.q0*omlam**4*Pi**2) + - (a0**2*a2**2*m**4)/(960.q0*omlam**4*Pi**2) - - (a0**2*a1*a3*m**4)/(480.q0*omlam**4*Pi**2) + - (a1**4*lam*m**2)/(192.q0*a0**2*omlam**3*Pi**2) - - (a1**2*lam**3*m**2)/(96.q0*a0**2*omlam**3*Pi**2) + - (a1**4*lam)/(96.q0*a0**4*omlam*Pi**2) return end real*16 function srint (i,k, yr, yrdot, yi, yidot, t, 1 lm, mtest) implicit none integer mtest,i real*16 k, yr, yrdot, yi, yidot, t, lm, m, pi, z, a, mass, 1 a0, a1, a2, a3, a4, sun, srn, h1, wk2, 2 y2, yp2, yyp, tol, yrnsq, trnsq, wkbsq, 3 t1, t2, mu, R, srinterror,ertest,kold common/param/m, pi, z, mu common/sq/yrnsq, trnsq, wkbsq, tol common/error/srinterror,ertest 150 format(F16.6, ES22.11, ES22.11) a0 = a(t,0) a1 = a(t,1) a2 = a(t,2) a3 = a(t,3) a4 = a(t,4) if (mtest.eq.0) then m = m else m = 0 end if h1 = - 36.q0*(a4*a0**4+7.q0*a3*a1*a0**3+11.q0*a2*a1**2*a0**2 1 + 4.q0*a2**2*a0**3+a1**4*a0)/a0**5 + 144.q0*a1*a0*(a3*a0**3 2 + 4.q0*a2*a1*a0**2 + a1**3*a0)/a0**6 3 - 216.q0*a1**2*a0**2*(a2*a0**2+a1**2*a0)/a0**7 + 108.q0*(a2*a0**2 4 + a1**2*a0)**2/a0**6 R = 6.q0*(a2*a0**2 + a1**2*a0)/a0**3 y2 = yr**2 + yi**2 yp2 = yrdot**2 + yidot**2 yyp = 2.q0*(yr*yrdot + yi*yidot) if (k.eq.0) then srint = 0 else wk2 = k**2 + m**2 * a0**2 sun = 1.q0/(2.q0 * pi**2 * a0**4) * k**2 * 1 (m**2*a0*y2 + (6.q0*z-1.q0)*(a0**2*(yp2/a0 - 2 yyp*a1/(2.q0*a0**2) + y2*a1**2/(4.q0*a0**3)) 3 - a1*(yyp - y2*a1/a0) - (wk2 + a2*a0 + 4 (z-1.q0/6.q0)*a0**2*R)*y2/a0)) c if (k.lt.lm) then If(qabs(k-lm).lt.1.q-10) then If(i.eq.2) then Print*,'lm is in the middle of an integration interval' Print*,'k = ',k,' lm = ',lm Stop Else if(i.eq.3) then srn = 1.q0/(4.q0 * pi**2 * a0**4) * k**2 * 1 ((m**2*a0**2 - (6.q0*z-1.q0)* 2 (a2*a0))/k) Else if(i.eq.1) then srn = 1.q0/(4.q0 * pi**2 * a0**4) * k**2 * 1 (((m**2*a0**2 - (6.q0*z-1.q0)* 2 (a2*a0))/k) + ((-m**4*a0**4/2.q0 3 - (z-1.q0/6.q0)*(3.q0*m**2*a0**2*(a2*a0+a1**2) + 4 (z-1.q0/6.q0)**2*h1*a0**4/4.q0))/k**3)) c This last case should never happen Else Print *,'From sint, i = ',i Stop Endif Else if(k.lt.lm) then srn = 1.q0/(4.q0 * pi**2 * a0**4) * k**2 * 1 ((m**2*a0**2 - (6.q0*z-1.q0)* 2 (a2*a0))/k) Else srn = 1.q0/(4.q0 * pi**2 * a0**4) * k**2 * 1 (((m**2*a0**2 - (6.q0*z-1.q0)* 2 (a2*a0))/k) + ((-m**4*a0**4/2.q0 3 - (z-1.q0/6.q0)*(3.q0*m**2*a0**2*(a2*a0+a1**2) + 4 (z-1.q0/6.q0)**2*h1*a0**4/4.q0))/k**3)) end if c This is the exact subtraction for I2 in Paper III srn = (m**2/(4.q0*pi**2*a0**2))*k**2/(2.q0*qsqrt(wk2)) c This was not working. See sint for more of an explanation c if(qabs((sun-srn))/qabs(sun+srn).lt.tol) then c srint = 0 c else srint = sun-srn c srinterror = -qlog10(qabs(srint)/(qabs(sun)+qabs(srn))) c end if end if return end real*16 function sran(t, lm, mtest) implicit none integer mtest real*16 t, lm, m, pi, z, a, a0, a1, a2, a3, a4, r, u, 1 h1, h3, w, tf, c, l, gm, mass, mu common/param/m, pi, z, mu common/scale/w a0 = a(t,0) a1 = a(t,1) a2 = a(t,2) a3 = a(t,3) a4 = a(t,4) h1 = - 36.q0*(a4*a0**4+7.q0*a3*a1*a0**3+11.q0*a2*a1**2*a0**2 1 + 4.q0*a2**2*a0**3+a1**4*a0)/a0**5 + 144.q0*a1*a0*(a3*a0**3 2 + 4.q0*a2*a1*a0**2 + a1**3*a0)/a0**6 3 - 216.q0*a1**2*a0**2*(a2*a0**2+a1**2*a0)/a0**7 + 108.q0*(a2*a0**2 4 + a1**2*a0)**2/a0**6 h3 = 12.q0*a2*a1**2/a0**3 c = -2.q0/(3.q0*w + 1.q0) if (mtest.eq.1) then u = mu else u = m end if if (lm.gt.0) then l = qlog(u**2*a0**2/(4.q0*lm**2)) else l = 0 end if if (mtest.eq.0) then m = m gm = m**2*(-6.q0*a2/a0 - 6.q0*a1**2/a0**2) sran = 1.q0/(2880.q0*pi**2)*((-1.q0/6.q0)*h1 + h3) 1 + gm/(288.q0*pi**2) - m**4*(1.q0+l)/(16.q0*pi**2) 2 + (z-1.q0/6.q0)*(h1/(288.q0*pi**2) + gm*(3.q0+l)/(16.q0*pi**2) 3 - 3.q0/(8.q0*pi**2)*m**2*a1**2/a0**2) 4 + (z-1.q0/6.q0)**2*(h1*(2.q0+l)/(32.q0*pi**2) 5 - 9.q0*(4.q0*a1*(a3*a0**3+4.q0*a2*a1*a0**2 6 + a1**3*a0)/a0**5 - 10.q0*a1**2*(a2*a0**2+a1**2*a0)/a0**5 6 + 3.q0*(a2*a0**2+a1**2*a0)**2/a0**6)/(8.q0*pi**2)) else m = 0 gm = m**2*(-6.q0*a2/a0 - 6.q0*a1**2/a0**2) sran = 1.q0/(2880.q0*pi**2)*((-1.q0/6.q0)*h1 + h3) 1 + gm/(288.q0*pi**2) - m**4*(1.q0+l)/(16.q0*pi**2) 2 + (z-1.q0/6.q0)*(h1/(288.q0*pi**2) + gm*(3.q0+l)/(16.q0*pi**2) 3 - 3.q0/(8.q0*pi**2)*m**2*a1**2/a0**2) 4 + (z-1.q0/6.q0)**2*(h1*(2.q0+l)/(32.q0*pi**2) 5 - 9.q0/(8.q0*pi**2)*(4.q0*a1*(a3*a0**2+4.q0*a2*a1*a0 6 + a1**3)/a0**5 - 10.q0*a1**2*(a2*a0**2+a1**2*a0)/a0**5 6 + 3.q0*(a2*a0**2+a1**2*a0)**2/a0**6)) end if c There is no contribution from the analytic approximation for this c value of I2. sran = 0.q0 return end