c This is the quadruple precision version c c This version computes 100 modes per decade so |B|^2 can be plotted c c This program solves the mode equation for f_k for the initially and finally static spacetime in spatially flat sections c It does so for multiple values of k and only prints out |B2|^2 at the largest c value of u. c c This program uses the same variable stepsize as modes_sf_2.f. It allows for integrations for larger values of T and k. c This version also starts the modes out at a later time for large values of k c and gives them the exact starting values appropriate for a conformally invariant c field. IMPLICIT NONE Integer j,jt,ik, nu,ib,nb, JUNK, IT, nbad,nok,nk,itest Real*16 u, ui,uin,delu, T,Y, m, k,gamma2,ki,delk Real*16 eps,h1,hmin Real*16 a,h,hp,Omega,om,omp,om2p,W,V,fr,fpr,fi,fpi,beta2 Real*16 ain,aout,om0,fr0,fpr0,fi0,fpi0,dif1,dif2,dif3,dif4 Real*16 for,fopr,foi,fopi,B2 Real*16 pi,delumin,deluj,test,utest,atest Real*16 htest,mak DIMENSION Y(4) COMMON/RealP/k COMMON/IntP/m,T External DERIVS, BSSTEP OPEN(UNIT = 11, FILE = 'modes_sf_B2_qp_plot.dat') Open(unit=12,file = 'B2.out') READ(11,*) JUNK c READ(11,*) ki, delk, nk, uin, delu, nu READ(11,*) ki, nb, uin, delu, nu READ(11,*) JUNK Read(11,*) m,T READ(11,*) JUNK Read(11,*) eps,h1,hmin pi = 2.q0*qasin(1.q0) gamma2 = m**2 - 0.25q0 ain = qexp(-T) aout = qexp(T) c Do loop for the bands of k c For this version there are 101 values of k in each decade Do ib = 1,nb delk = ki*(10.q0**(ib-2)) k = ki*10.q0**(ib-1) Write(6,11) ib, k,delk Enddo Do ib = 1,nb delk = ki*(10.q0**(ib-2)) k = ki*10.q0**(ib-1)-delk Do ik = 1,91 k = k+delk ui = uin c This is the static approximation at the initial time c We'll adjust the phase so that it is the same as for the actual modes om0 = qsqrt(k**2/ain**2 + m**2) fr0 = 1.q0/qsqrt(2.q0*om0) fpr0 = 0.q0 fi0 = 0.q0 fpi0 = -om0/qsqrt(2.q0*om0) a = qexp(T*qtanh(ui/T)) h = 1.q0/qcosh(ui/T)**2 hp = - ((2.q0/qcosh(ui/T)**2)*qtanh(ui/T))/T Omega = qsqrt(k**2/a**2 + m**2 - hp/2.q0 - h**2/4.q0) om = qsqrt(k**2/a**2 + m**2) omp = -(k**2/qcosh(ui/T)**2)*qexp(-2.q0*T*qtanh(ui/T))/om om2p = (2.q0*k**2/qcosh(ui/T)**4)*qexp(-2.q0*T*qtanh(ui/T))/om + - (2.q0*k**2/qcosh(ui/T)**2)*qtanh(ui/T)*qexp(-2.q0*T*qtanh(ui/T))/ - (om*T) - ((k**4/qcosh(ui/T)**4)*qexp(-4.q0*T*qtanh(ui/T)))/om**3 W = Omega + 3.q0*omp**2/(8.q0*om**3) - om2p/(4.q0*om**2) V = -omp/om c Y(1) is the real part of f, Y(2) the real part of df/dt, Y(3) the imaginary part of f, and Y(4) the imaginary part of df/dt Y(1)= 1.q0/qsqrt(2.q0*W) Y(2)= V/(2.q0**1.5q0*qsqrt(W)) Y(3) = 0.q0 Y(4) = -qsqrt(W/2.q0) fr = Y(1) fpr = Y(2) fi = Y(3) fpi = Y(4) c Now check to see if the mode is effectively conformally invariant. If(m*a/k.lt.1.e-5) then itest = 0 Else itest = 1 Endif Print *,'k = ',k,' Initial value of itest is ',itest DO IT = 1,nu utest = ui atest = qexp(T*qtanh(utest/T)) htest = 1.q0/qcosh(utest/T)**2 c First check to see if we want to do a computation at this time step at all. If(itest.eq.0) then mak = m*atest/k c The 1.e-8 gives double precision accuracy to having a conformally invariant c field. However, we don't necessarily need that. In fact 7 or 8 digit c accuracy should be sufficient. I'll try that. c If(mak.lt.1.e-8) then c If(mak.lt.1.e-4) then If(mak.lt.1.q-5) then c If(mak.lt.1.e-6) then c If(mak.lt.1.e-7) then c Print*,'u = ',utest,' ma/k = ',mak ui = utest+delu Go to 100 Else itest = 1 om = k/atest c Note these are exact starting values for a conformally invariant scalar field Y(1)= 1.q0/qsqrt(2.q0*om) Y(2)= 0.5q0*htest/(qsqrt(2.q0*om)) Y(3) = 0.q0 Y(4) = -qsqrt(om/2.q0) Print*,'mode integration starts after u = ',ui Endif Endif c Now we want to decrease the time interval further if it is too large delumin = 2.q0*pi*atest/qsqrt(k**2+m**2*atest**2) test = delu/delumin If(test.gt.1.q3) then Print*, 'Extra values of u are being used for the numerical integration' jt = test deluj = delu/jt Print *,'u = ',u,' jt = ',jt,' delu = ',delu,' deluj = ',deluj Else jt = 1 deluj = delu Endif Do j = 1,jt u = ui + deluj CALL ODEINT(Y,4,ui,u,EPS,H1,HMIN,NOK,NBAD,DERIVS,BSSTEP) ui = u Enddo a = qexp(T*qtanh(u/T)) h = 1.q0/qcosh(u/T)**2 hp = - (2.q0/qcosh(u/T)**2)*qtanh(u/T)/T Omega = qsqrt(k**2/a**2 + m**2 - hp/2.q0 - h**2/4.q0) om = qsqrt(k**2/a**2 + m**2) omp = -(k**2/qcosh(u/T)**2)*qexp(-2.q0*T*qtanh(u/T))/om om2p = (2.q0*k**2/qcosh(u/T)**4)*qexp(-2.q0*T*qtanh(u/T))/om + - (2.q0*k**2/qcosh(u/T)**2)*qtanh(u/T)*qexp(-2.q0*T*qtanh(u/T))/ - (om*T) - ((k**4/qcosh(u/T)**4)*qexp(-4.q0*T*qtanh(u/T)))/om**3 W = Omega + 3.q0*omp**2/(8.q0*om**3) - om2p/(4.q0*om**2) V = -omp/om fr = Y(1) fpr = Y(2) fi = Y(3) fpi = Y(4) c I'll match with the static out state here as a function of time starting at u= 10 c This is different than in the previous version c Note that the value of the mode function in the out c state if I just choose the obvious phase is the c same as it is for the in state since we are talking c about modes in a static spacetime If(it.eq.nu) then om0 = qsqrt(k**2/aout**2 + m**2) for = 1.q0/qsqrt(2.q0*om0) fopr = 0.q0 foi = 0.q0 fopi = -om0/qsqrt(2.q0*om0) B2 = fi**2*fopi**2 + fi**2*fopr**2 - 2*fi*foi*fopi*fpi - - 2*fi*fopr*for*fpi + foi**2*fpi**2 + for**2*fpi**2 - - 2*fi*foi*fopr*fpr + 2*fi*fopi*for*fpr + foi**2*fpr**2 + - for**2*fpr**2 + 2*foi*fopr*fpi*fr - 2*fopi*for*fpi*fr - - 2*foi*fopi*fpr*fr - 2*fopr*for*fpr*fr + fopi**2*fr**2 + - fopr**2*fr**2 Write(12,35) u,k,B2 Endif 100 Continue Enddo Enddo Enddo 11 Format(1x,I5,1x,5(d12.5,1x)) 35 Format(1x,10(d14.7,1x)) STOP END ! -------------------------------------------------------------- ! -------------- SUBROUTINES ----------------------------------- ! -------------------------------------------------------------- SUBROUTINE DERIVS(u,Y,DYDT,NVAR) INTEGER nvar REAL*16 u, T, Y(4), DYDT(4), k, m, a, h, hp, Omega2 COMMON/RealP/k COMMON/IntP/m,T a = qexp(T*qtanh(u/T)) h = 1.q0/qcosh(u/T)**2 hp = - (2.q0/qcosh(u/T)**2)*qtanh(u/T)/T Omega2 = k**2/a**2 + m**2 - hp/2.q0 - h**2/4.q0 c Print*,u,a,h,hp,Omega2 DYDT(1) = Y(2) DYDT(2) = -Omega2*Y(1) DYDT(3) = Y(4) DYDT(4) = -Omega2*Y(3) RETURN END