! This program takes the values of the Bunch-Davies modes at u = -1 and ! uses them as initial conditions for the mode equation. The mode equation is ! then solved in terms of proper time u to get the exact Bunch-Davies modes for all times. c I've decided for now the best way to do it is to have different versions. This version c is for m = 1 and k = 1000. There will be 6 combinations total. Implicit None Integer i,im, nu, it, nk, nbad,nok,junk real*8 ustart,uin,u,delu,m,uout,y(4) Real*8 pi,eps,h1,hmin,k Real*8 gamma Real*8 Ckr,Cki,vr,vpr,vi,vpi,B2k Real*8 W,V,om,Omega,fr,fi,fpr,fpi,beta2 common/intp/k common/mass/m EXTERNAL DERIVS EXTERNAL BSSTEP Open(unit = 9, file = "BDmodes_v.dat") open(unit = 10, file = "BDmodes_v.out") open(unit = 11, file = "C_m_1_k_1000.out") open(unit = 12, file = "B_m_1_k_1000.out") open(unit = 13, file = "beta2BD_m_1_k_1000.out") read(9,*) junk read(9,*) uin,delu,nu read(9,*) junk read(9,*),eps,h1,hmin pi = 2.d0*dasin(1.d0) Print *,'pi = ',pi m = 1.d0 gamma = dsqrt(m**2 - 0.25d0) k = 1.d3 c We start the integration at u = -1 ustart = -1.d0 c Y(1) is vr, Y(2) is vpr, Y(3) is vi, Y(4) is vpi y(1) = 0.003743463873367724d0 y(2) = -9.069303933746404d0 y(3) = -0.01399904993336666d0 y(4) = -2.436632256977701d0 CALL ODEINT(Y,4,ustart,uin,EPS,H1,HMIN,NOK,NBAD,DERIVS,BSSTEP) u = uin uout = uin + delu vr = Y(1) vpr = Y(2) vi = Y(3) vpi = Y(4) c Now compute the time dependent beta(u) using Emil's second order definition. om = dsqrt(k**2/dcosh(uin)**2 + m**2) Omega = dsqrt((k**2-0.25d0)/dcosh(uin)**2 + gamma**2) W = Omega + (dtanh(uin)**2/(8.d0*om))*(1.d0-6.d0*m**2/om**2 - + 5.d0*m**4/om**4) - + (1.d0/(4.d0*om*dcosh(uin)**2))*(1.d0-m**2/om**2) V = dtanh(uin)*(1.d0 - m**2/om**2) c Remember beta2 is in terms of f = a^(3/2) v so I need to convert. c This comes from stwkbin_emil which is in terms of f fr = dcosh(uin)**1.5d0*vr fi = dcosh(uin)**1.5d0*vi fpr = dcosh(uin)**1.5d0*vpr + 1.5d0*dsqrt(dcosh(uin))*dsinh(uin)*vr fpi = dcosh(uin)**1.5d0*vpi + 1.5d0*dsqrt(dcosh(uin))*dsinh(uin)*vi beta2 = (1.d0/(2.d0*W))*(fpr**2+fpi**2 - - 2.d0*W*fpr*fi + 2.d0*W*fpi*fr - - V*(fr*fpr+fi*fpi) - + (V**2/4.d0 + W**2)*(fr**2+fi**2)) write(10,*) k, u, y(1), y(2), y(3), y(4) write(13,29) k, u, beta2 c Now compute the coefficients of Ckr and Cki Ckr = -(k**2*vi**2)/(8.d0*Pi**2) - (gamma**2*k**2*vi**2)/(2.d0*Pi**2) - - (k**2*vpi**2)/(2.d0*Pi**2) + (k**2*vpr**2)/(2.d0*Pi**2) + - (k**2*vr**2)/(8.d0*Pi**2) + (gamma**2*k**2*vr**2)/(2.d0*Pi**2) - - (k**4*vi**2)/(2.d0*Pi**2*dcosh(u)**2) + (k**4*vr**2)/(2.d0*Pi**2*dcosh(u)**2) - - (k**2*vi*vpi*dtanh(u))/Pi**2 + (k**2*vpr*vr*dtanh(u))/Pi**2 - - (k**2*vi**2*dtanh(u)**2)/(2.d0*Pi**2) + (k**2*vr**2*dtanh(u)**2)/(2.d0*Pi**2) Cki = -((k**2*vpi*vpr)/Pi**2) - (k**2*vi*vr)/(4.d0*Pi**2) - - (gamma**2*k**2*vi*vr)/Pi**2 - (k**4*vi*vr)/(Pi**2*dcosh(u)**2) - - (k**2*vi*vpr*dtanh(u))/Pi**2 - (k**2*vpi*vr*dtanh(u))/Pi**2 - - (k**2*vi*vr*dtanh(u)**2)/Pi**2 Write(11,29) u,k,Ckr,Cki B2k = (k**2*vi**2)/(8.d0*Pi**2) + (gamma**2*k**2*vi**2)/(2.d0*Pi**2) + - (k**2*vr**2)/(8.d0*Pi**2) + (gamma**2*k**2*vr**2)/(2.d0*Pi**2) + - (k**4*vi**2/dcosh(u)**2)/(2.d0*Pi**2) + - (k**4*vr**2/dcosh(u)**2)/(2.d0*Pi**2) + - (k**2*vi**2*dtanh(u)**2)/(2.d0*Pi**2) + (k**2*vr**2*dtanh(u)**2)/(2.d0*Pi**2) - + k**2*(vpr**2+vpi**2 + 2.d0*dtanh(u)*(vr*vpr+vi*vpi))/(2.d0*Pi**2) Write(12,29) u,k,B2k do it = 1,nu CALL ODEINT(Y,4,u,uout,EPS,H1,HMIN,NOK,NBAD,DERIVS,BSSTEP) write(10,*) k, uout, y(1), y(2), y(3), y(4) u = uout uout = uout + delu vr = Y(1) vpr = Y(2) vi = Y(3) vpi = Y(4) c Now compute the time dependent beta(u) using Emil's second order definition. om = dsqrt(k**2/dcosh(u)**2 + m**2) Omega = dsqrt((k**2-0.25d0)/dcosh(u)**2 + gamma**2) W = Omega + (dtanh(u)**2/(8.d0*om))*(1.d0-6.d0*m**2/om**2 - + 5.d0*m**4/om**4) - + (1.d0/(4.d0*om*dcosh(u)**2))*(1.d0-m**2/om**2) V = dtanh(u)*(1.d0 - m**2/om**2) c Remember beta2 is in terms of f = a^(3/2) v so I need to convert. c This comes from stwkbin_emil which is in terms of f fr = dcosh(u)**1.5d0*vr fi = dcosh(u)**1.5d0*vi fpr = dcosh(u)**1.5d0*vpr + 1.5d0*dsqrt(dcosh(u))*dsinh(u)*vr fpi = dcosh(u)**1.5d0*vpi + 1.5d0*dsqrt(dcosh(u))*dsinh(u)*vi beta2 = (1.d0/(2.d0*W))*(fpr**2+fpi**2 - - 2.d0*W*fpr*fi + 2.d0*W*fpi*fr - - V*(fr*fpr+fi*fpi) - + (V**2/4.d0 + W**2)*(fr**2+fi**2)) write(10,*) k, u, y(1), y(2), y(3), y(4) write(13,29) k, u, beta2 c Now compute the coefficients of Ckr and Cki Ckr = -(k**2*vi**2)/(8.d0*Pi**2) - (gamma**2*k**2*vi**2)/(2.d0*Pi**2) - - (k**2*vpi**2)/(2.d0*Pi**2) + (k**2*vpr**2)/(2.d0*Pi**2) + - (k**2*vr**2)/(8.d0*Pi**2) + (gamma**2*k**2*vr**2)/(2.d0*Pi**2) - - (k**4*vi**2)/(2.d0*Pi**2*dcosh(u)**2) + (k**4*vr**2)/(2.d0*Pi**2*dcosh(u)**2) - - (k**2*vi*vpi*dtanh(u))/Pi**2 + (k**2*vpr*vr*dtanh(u))/Pi**2 - - (k**2*vi**2*dtanh(u)**2)/(2.d0*Pi**2) + (k**2*vr**2*dtanh(u)**2)/(2.d0*Pi**2) Cki = -((k**2*vpi*vpr)/Pi**2) - (k**2*vi*vr)/(4.d0*Pi**2) - - (gamma**2*k**2*vi*vr)/Pi**2 - (k**4*vi*vr)/(Pi**2*dcosh(u)**2) - - (k**2*vi*vpr*dtanh(u))/Pi**2 - (k**2*vpi*vr*dtanh(u))/Pi**2 - - (k**2*vi*vr*dtanh(u)**2)/Pi**2 Write(11,29) u,k,Ckr,Cki B2k = (k**2*vi**2)/(8.d0*Pi**2) + (gamma**2*k**2*vi**2)/(2.d0*Pi**2) + - (k**2*vr**2)/(8.d0*Pi**2) + (gamma**2*k**2*vr**2)/(2.d0*Pi**2) + - (k**4*vi**2/dcosh(u)**2)/(2.d0*Pi**2) + - (k**4*vr**2/dcosh(u)**2)/(2.d0*Pi**2) + - (k**2*vi**2*dtanh(u)**2)/(2.d0*Pi**2) + (k**2*vr**2*dtanh(u)**2)/(2.d0*Pi**2) - + k**2*(vpr**2+vpi**2 + 2.d0*dtanh(u)*(vr*vpr+vi*vpi))/(2.d0*Pi**2) Write(12,29) u,k,B2k end do 29 Format(1x,20(d12.5,2x)) stop end SUBROUTINE DERIVS(u,Y,DYDT,NVAR) INTEGER nvar REAL*8 u, Y(4), DYDT(4), k, m, a, h COMMON/IntP/k COMMON/mass/m a = dcosh(u) h = dtanh(u) DYDT(1) = Y(2) DYDT(2) = -3.d0*h*Y(2)-((k**2-1.d0)/dcosh(u)**2 + (m**2+2.d0))*Y(1) DYDT(3) = Y(4) DYDT(4) = -3.d0*h*Y(4)-((k**2-1.d0)/dcosh(u)**2 + (m**2+2.d0))*Y(3) RETURN END