c This is a double precision version c This alternate version prints out beta2 for a given value of k. c This version uses the new profile with parameters u1 and b. c This version uses the 2nd order form of W that Emil prefers 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, JUNK, IT, nbad,nok,nk,itest Real*8 u, ui,uin,delu, Y, m, k,gamma2,ki,delk Real*8 eps,h1,hmin Real*8 a,h,hp,Omega,om,omp,om2p,W,V,fr,fpr,fi,fpi,beta2 Real*8 ain,aout,om0,fr0,fpr0,fi0,fpi0,dif1,dif2,dif3,dif4 Real*8 for,fopr,foi,fopi,B2 Real*8 pi,delumin,deluj,test,utest,atest Real*8 htest,mak Real*8 u1,b Real*8 TwoWm1,f02 Real*8 int,intrho,Are,Aim,Bre,Bim,B2new,B2dif,B2test,ABtest,om0test Real*8 statetest,hptest DIMENSION Y(4) COMMON/RealP/k COMMON/IntP/m,u1,b External DERIVS, BSSTEP OPEN(UNIT = 11, FILE = 'modes_sf_B2_new.dat') Open(unit=12,file = 'B2dp.out') Open(unit=13,file='modesdp.tst') Open(unit=14,file='B2dp.tst') Open(unit=15,file='ABdp.out') c Open(unit=14,file='beta2.out') READ(11,*) JUNK READ(11,*) ki, delk, nk, uin, delu, nu READ(11,*) JUNK Read(11,*) m,u1,b READ(11,*) JUNK Read(11,*) eps,h1,hmin pi = 2.d0*dasin(1.d0) gamma2 = m**2 - 0.25d0 ain = dexp(-u1) aout = dexp(u1) k = ki-delk Do ik = 1,nk 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 = dsqrt(k**2/ain**2 + m**2) fr0 = 1.d0/dsqrt(2.d0*om0) fpr0 = 0.d0 fi0 = 0.d0 fpi0 = -om0/dsqrt(2.d0*om0) a = dexp((-dlog(dcosh(b*(ui-u1)))+dlog(dcosh(b*(ui+u1))))/(2.d0*b)) h = -dtanh(b*(ui - u1))/2.d0 + dtanh(b*(ui + u1))/2.d0 hp = -(b/dcosh(b*(ui - u1))**2)/2.d0 + (b/dcosh(b*(ui + u1))**2)/2.d0 Omega = dsqrt(k**2/a**2 + m**2 - hp/2.d0 - h**2/4.d0) om = dsqrt(k**2/a**2 + m**2) omp = (k**2*dcosh(b*(ui - u1))**(-1.d0 + 1.d0/b)*dsinh(b*(ui - u1)))/ - (2.d0*dcosh(b*(ui + u1))**(1.d0/b)* - dsqrt(m**2 + (k**2*dcosh(b*(ui - u1))**(1.d0/b))/dcosh(b*(ui + u1))**(1.d0/b))) - - (k**2*dcosh(b*(ui - u1))**(1.d0/b)*dcosh(b*(ui + u1))**(-1.d0 - 1.d0/b)* - dsinh(b*(ui + u1)))/ - (2.d0*dsqrt(m**2 + (k**2*dcosh(b*(ui - u1))**(1.d0/b))/dcosh(b*(ui + u1))**(1.d0/b))) om2p = -(k**4*dcosh(b*(ui - u1))**(-2.d0 + 2.d0/b)*dsinh(b*(ui - u1))**2)/ - (4.d0*dcosh(b*(ui + u1))**(2.d0/b)* - (m**2 + (k**2*dcosh(b*(ui - u1))**(1.d0/b))/dcosh(b*(ui + u1))**(1.d0/b))**1.5d0) - + (k**2*dcosh(b*(ui - u1))**(-2.d0 + 1.d0/b)*dsinh(b*(ui - u1))**2)/ - (2.d0*dcosh(b*(ui + u1))**(1.d0/b)* - dsqrt(m**2 + (k**2*dcosh(b*(ui - u1))**(1.d0/b))/dcosh(b*(ui + u1))**(1.d0/b))) - - (b*k**2*dcosh(b*(ui - u1))**(-2.d0 + 1.d0/b)*dsinh(b*(ui - u1))**2)/ - (2.d0*dcosh(b*(ui + u1))**(1.d0/b)* - dsqrt(m**2 + (k**2*dcosh(b*(ui - u1))**(1.d0/b))/dcosh(b*(ui + u1))**(1.d0/b))) + - (k**4*dcosh(b*(ui - u1))**(-1.d0 + 2.d0/b)*dcosh(b*(ui + u1))**(-1.d0 - 2.d0/b)* - dsinh(b*(ui - u1))*dsinh(b*(ui + u1)))/ - (2.d0*(m**2 + (k**2*dcosh(b*(ui - u1))**(1.d0/b))/dcosh(b*(ui + u1))**(1.d0/b))**1.5d0) - - (k**2*dcosh(b*(ui - u1))**(-1.d0 + 1.d0/b)*dcosh(b*(ui + u1))**(-1.d0 - 1.d0/b)* - dsinh(b*(ui - u1))*dsinh(b*(ui + u1)))/ - dsqrt(m**2 + (k**2*dcosh(b*(ui - u1))**(1.d0/b))/dcosh(b*(ui + u1))**(1.d0/b)) - - (k**4*dcosh(b*(ui - u1))**(2.d0/b)*dcosh(b*(ui + u1))**(-2.d0 - 2.d0/b)* - dsinh(b*(ui + u1))**2)/ - (4.d0*(m**2 + (k**2*dcosh(b*(ui - u1))**(1.d0/b))/dcosh(b*(ui + u1))**(1.d0/b))**1.5d0) - + (k**2*dcosh(b*(ui - u1))**(1.d0/b)*dcosh(b*(ui + u1))**(-2.d0 - 1.d0/b)* - dsinh(b*(ui + u1))**2)/ - (2.d0*dsqrt(m**2 + (k**2*dcosh(b*(ui - u1))**(1.d0/b))/dcosh(b*(ui + u1))**(1.d0/b))) - + (b*k**2*dcosh(b*(ui - u1))**(1.d0/b)*dcosh(b*(ui + u1))**(-2.d0 - 1.d0/b)* - dsinh(b*(ui + u1))**2)/ - (2.d0*dsqrt(m**2 + (k**2*dcosh(b*(ui - u1))**(1.d0/b))/dcosh(b*(ui + u1))**(1.d0/b))) c Write(13,*) b,u1,ui c Write(13,*) a,h,hp,om,omp,om2p W = Omega + 3.d0*omp**2/(8.d0*om**3) - om2p/(4.d0*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.d0/dsqrt(2.d0*W) Y(2)= V/(2.d0**1.5d0*dsqrt(W)) Y(3) = 0.d0 Y(4) = -dsqrt(W/2.d0) fr = Y(1) fpr = Y(2) fi = Y(3) fpi = Y(4) c Compare with the static modes - % diff for 1 and 4 and absolute diff for 2. dif1 = (fr-fr0)/fr dif2 = (fpr-fpr0) dif3 = (fr-fr0) dif4 = (fpi-fpi0)/fpi Write(13,*) k,u1,b,dif1,dif2,dif3,dif4 c Now check to see if the mode is effectively conformally invariant. c Originally I used ma/k < 1.d-5. But what I really want to do is to c use the accuracy the modes are computed to as the accuracy for the c state. That is eps from the modes.dat file. statetest = dsqrt((m**2-hp/2.d0-h**2/4.d0)/eps)*a/k Print *,' k = ',k,' statetest = ',statetest If(statetest.le.1.d0) then itest = 0 Else itest = 1 Endif Print *,'k = ',k,' Initial value of itest is ',itest Do IT = 1,nu utest = ui atest = dexp((-dlog(dcosh(b*(utest-u1)))+dlog(dcosh(b*(utest+u1))))/(2.d0*b)) htest = -dtanh(b*(utest - u1))/2.d0 + dtanh(b*(utest + u1))/2.d0 hptest = -(b/dcosh(b*(utest- u1))**2)/2.d0 + (b/dcosh(b*(utest + u1))**2)/2.d0 statetest = dsqrt((m**2-hptest/2.d0-htest**2/4.d0)/eps)*atest/k c First check to see if we want to do a computation at this time step at all. If(itest.eq.0) then c mak = m*atest/k c I now use the statetest variable - see above c 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.d-8) then c If(mak.lt.1.d-4) then c If(mak.lt.1.d-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 If(statetest.lt.1.d0) then 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.d0/dsqrt(2.d0*om) Y(2)= 0.5d0*htest/(dsqrt(2.d0*om)) Y(3) = 0.d0 Y(4) = -dsqrt(om/2.d0) 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.d0*pi*atest/dsqrt(k**2+m**2*atest**2) test = delu/delumin If(test.gt.1.d3) then c Print*, 'Extra values of u are being used for the numerical integration' jt = test deluj = delu/jt c 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 c The following is to compute |beta|^2 a = dexp((-dlog(dcosh(b*(u-u1)))+dlog(dcosh(b*(u+u1))))/(2.d0*b)) h = -dtanh(b*(u - u1))/2.d0 + dtanh(b*(u + u1))/2.d0 hp = -(b/dcosh(b*(u - u1))**2)/2.d0 + (b/dcosh(b*(u + u1))**2)/2.d0 Omega = dsqrt(k**2/a**2 + m**2 - hp/2.d0 - h**2/4.d0) om = dsqrt(k**2/a**2 + m**2) omp = (k**2*dcosh(b*(u - u1))**(-1.d0 + 1.d0/b)*dsinh(b*(u - u1)))/ - (2.d0*dcosh(b*(u + u1))**(1.d0/b)* - dsqrt(m**2 + (k**2*dcosh(b*(u - u1))**(1.d0/b))/dcosh(b*(u + u1))**(1.d0/b))) - - (k**2*dcosh(b*(u - u1))**(1.d0/b)*dcosh(b*(u + u1))**(-1.d0 - 1.d0/b)* - dsinh(b*(u + u1)))/ - (2.d0*dsqrt(m**2 + (k**2*dcosh(b*(u - u1))**(1.d0/b))/dcosh(b*(u + u1))**(1.d0/b))) om2p = -(k**4*dcosh(b*(u - u1))**(-2.d0 + 2.d0/b)*dsinh(b*(u - u1))**2)/ - (4.d0*dcosh(b*(u + u1))**(2.d0/b)* - (m**2 + (k**2*dcosh(b*(u - u1))**(1.d0/b))/dcosh(b*(u + u1))**(1.d0/b))**1.5d0) - + (k**2*dcosh(b*(u - u1))**(-2.d0 + 1.d0/b)*dsinh(b*(u - u1))**2)/ - (2.d0*dcosh(b*(u + u1))**(1.d0/b)* - dsqrt(m**2 + (k**2*dcosh(b*(u - u1))**(1.d0/b))/dcosh(b*(u + u1))**(1.d0/b))) - - (b*k**2*dcosh(b*(u - u1))**(-2.d0 + 1.d0/b)*dsinh(b*(u - u1))**2)/ - (2.d0*dcosh(b*(u + u1))**(1.d0/b)* - dsqrt(m**2 + (k**2*dcosh(b*(u - u1))**(1.d0/b))/dcosh(b*(u + u1))**(1.d0/b))) + - (k**4*dcosh(b*(u - u1))**(-1.d0 + 2.d0/b)*dcosh(b*(u + u1))**(-1.d0 - 2.d0/b)* - dsinh(b*(u - u1))*dsinh(b*(u + u1)))/ - (2.d0*(m**2 + (k**2*dcosh(b*(u - u1))**(1.d0/b))/dcosh(b*(u + u1))**(1.d0/b))**1.5d0) - - (k**2*dcosh(b*(u - u1))**(-1.d0 + 1.d0/b)*dcosh(b*(u + u1))**(-1.d0 - 1.d0/b)* - dsinh(b*(u - u1))*dsinh(b*(u + u1)))/ - dsqrt(m**2 + (k**2*dcosh(b*(u - u1))**(1.d0/b))/dcosh(b*(u + u1))**(1.d0/b)) - - (k**4*dcosh(b*(u - u1))**(2.d0/b)*dcosh(b*(u + u1))**(-2.d0 - 2.d0/b)* - dsinh(b*(u + u1))**2)/ - (4.d0*(m**2 + (k**2*dcosh(b*(u - u1))**(1.d0/b))/dcosh(b*(u + u1))**(1.d0/b))**1.5d0) - + (k**2*dcosh(b*(u - u1))**(1.d0/b)*dcosh(b*(u + u1))**(-2.d0 - 1.d0/b)* - dsinh(b*(u + u1))**2)/ - (2.d0*dsqrt(m**2 + (k**2*dcosh(b*(u - u1))**(1.d0/b))/dcosh(b*(u + u1))**(1.d0/b))) - + (b*k**2*dcosh(b*(u - u1))**(1.d0/b)*dcosh(b*(u + u1))**(-2.d0 - 1.d0/b)* - dsinh(b*(u + u1))**2)/ - (2.d0*dsqrt(m**2 + (k**2*dcosh(b*(u - u1))**(1.d0/b))/dcosh(b*(u + u1))**(1.d0/b))) c Write(13,*) b,u1,u c Write(13,*) a,h,hp,om,omp,om2p W = Omega + 3.d0*omp**2/(8.d0*om**3) - om2p/(4.d0*om**2) V = -omp/om fr = Y(1) fpr = Y(2) fi = Y(3) fpi = Y(4) beta2 = -(fi*fpr) + fpi*fr + fpi**2/(2.d0*W) + fpr**2/(2.d0*W) - (fi*fpi*V)/(2.d0*W) - - (fpr*fr*V)/(2.d0*W) + (fi**2*V**2)/(8.d0*W) + (fr**2*V**2)/(8.d0*W) + - (fi**2*W)/2.d0 + (fr**2*W)/2.d0 c Write(14,35) a,k,u,om,Omega,W,V,beta2 c Write(6,35) a,k,u,om,omp,om2p,W c If(it.eq.nu) then c Write(14,35) k,u,beta2 c Write(6,35) a,u,k,beta2 c Endif 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 c Write(14,35) fr,fpr,fi,fpi fr = Y(1) fpr = Y(2) fi = Y(3) fpi = Y(4) c Write(14,35) fr,fpr,fi,fpi om0 = dsqrt(k**2/aout**2 + m**2) for = 1.d0/dsqrt(2.d0*om0) fopr = 0.d0 foi = 0.d0 fopi = -om0/dsqrt(2.d0*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 c Now do all of the coefficients Are = fopr*fi - for*fpi + foi*fpr - fopi*fr Aim = -(fopi*fi) + foi*fpi + for*fpr - fopr*fr Bre = -(fopr*fi) + for*fpi + foi*fpr - fopi*fr Bim = -(fopi*fi) + foi*fpi - for*fpr + fopr*fr c Now let's test them B2test = Bre**2+Bim**2 B2dif = (B2-B2test)/B2 ABtest = Are**2+Aim**2-B2test-1.d0 Write(14,35) k, B2dif,ABtest If(dabs(B2dif).gt.1.d-12) then Write(14,*) ' Problem with old and new B2, Relative difference = ',B2dif Write(14,*), ' B2old = ',B2,' B2new = ',B2test Endif c This in some sense checks the accuracy of the mode solutions c I've decided 10 digits is enough for now. If(dabs(ABtest).gt.1.d-8) then Write(14,*), ' Problem with ABtest. It is ',ABtest Endif int = k**2*dlog(1.d0+B2) intrho = k**2*om0*B2 c This is for plotting the integrands only. I'll write a driver c program to do the numerical integrations using the exact modes Write(12,35) aout,u,k,om0,B2,int,intrho Write(6,35) aout,u,k,B2,int Write(15,*) k,Are,Aim,Bre,Bim c Note I did this and found no difference to quadruple precision c Now let's try matching at the specific value of u by using the c correct scale factor for the omega for the matching coefficients c om0test = dsqrt(k**2/a**2 + m**2) c for = 1.d0/dsqrt(2.d0*om0test) c fopr = 0.d0 c foi = 0.d0 c fopi = -om0test/dsqrt(2.d0*om0test) c Now do all of the coefficients c Are = fopr*fi - for*fpi + foi*fpr - fopi*fr c Aim = -(fopi*fi) + foi*fpi + for*fpr - fopr*fr c Bre = -(fopr*fi) + for*fpi + foi*fpr - fopi*fr c Bim = -(fopi*fi) + foi*fpi - for*fpr + fopr*fr c Now let's test them c B2test = Bre**2+Bim**2 c B2dif = (B2-B2test)/B2 c ABtest = Are**2+Aim**2-B2test-1.d0 c Write(14,*) ' These are matched with the actual value of the scale factor used ' c Write(14,*) ' New ABtest = ',ABtest c Write(12,35) for,fopr,foi,fopi c Write(12,35) fr,fpr,fi,fpi c f02 = for**2+foi**2 c TwoWm1 = 1.d0/(2.d0*W) c Write(12,35) f02,TwoWm1 Endif 100 Continue Enddo Enddo 35 Format(1x,10(e16.7,1x)) STOP END ! -------------------------------------------------------------- ! -------------- SUBROUTINES ----------------------------------- ! -------------------------------------------------------------- SUBROUTINE DERIVS(u,Y,DYDT,NVAR) INTEGER nvar REAL*8 u, u1, b,Y(4), DYDT(4), k, m, a, h, hp, Omega2 COMMON/RealP/k COMMON/IntP/m,u1,b a = dexp((-dlog(dcosh(b*(u-u1)))+dlog(dcosh(b*(u+u1))))/(2.d0*b)) h = -dtanh(b*(u - u1))/2.d0 + dtanh(b*(u + u1))/2.d0 hp = -(b/dcosh(b*(u - u1))**2)/2.d0 + (b/dcosh(b*(u + u1))**2)/2.d0 Omega2 = k**2/a**2 + m**2 - hp/2.d0 - h**2/4.d0 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