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 IMPLICIT NONE Integer j,ik, nu, JUNK, IT, nbad,nok,nT Real*8 u, ui,uin,delu, T,Y, m, k,gamma2 Real*8 eps,h1,hmin Real*8 a,h,hp,Omega,om,omp,om2p,W,V,fr,fpr,fi,fpi,beta2 Real*8 a0,om0,fr0,fpr0,fi0,fpi0,dif1,dif2,dif3,dif4 Real*8 Tin,delT Real*8 for,fopr,foi,fopi,B2 DIMENSION Y(4) COMMON/RealP/k COMMON/IntP/m,T External DERIVS, BSSTEP OPEN(UNIT = 11, FILE = 'modes_s.dat') c Note writes to 12 are commented out right now as arerights to 14 OPEN(UNIT = 12, FILE = 'mode_s.out') Open(unit=13,file = 'B2.out') Open(unit=14,file = 'modes_s.tst') Open(unit=15,file = 'beta2.tst') READ(11,*) JUNK READ(11,*) ik, uin, delu, nu READ(11,*) JUNK Read(11,*) m,Tin,delT,nT READ(11,*) JUNK Read(11,*) eps,h1,hmin c Write(14,*)' k T ui (fr-fr0)/fr fpr-fpr0 (fpi - fpi0)/fpi' gamma2 = m**2 - 0.25d0 T = Tin - delT Do j = 1,nT T = T + delT c This is for the static approximation a0 = dcosh(T) ui = uin k =ik 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/a0**2 + m**2) fr0 = 1.d0/dsqrt(2.d0*om0) fpr0 = 0.d0 fi0 = 0.d0 fpi0 = -om0/dsqrt(2.d0*om0) a = dcosh(T*dtanh(ui/T)) h = (1.d0/dcosh(ui/T)**2)*dtanh(T*dtanh(ui/T)) hp = (1.d0/dcosh(ui/T)**4)*(1.d0/dcosh(T*dtanh(ui/T))**2) - - ((2.d0/dcosh(ui/T)**2)*dtanh(ui/T)*dtanh(T*dtanh(ui/T)))/T 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(ui/T)**2)*(1.d0/dcosh(T*dtanh(ui/T))**2)*dtanh(T*dtanh(ui/T)))/om) om2p = -(((k**2/dcosh(ui/T)**4)*(1.d0/dcosh(T*Tanh(ui/T))**4))/om) + - ((2.d0*k**2/dcosh(ui/T)**2)*(1.d0/dcosh(T*dtanh(ui/T))**2)*dtanh(ui/T)*dtanh(T*dtanh(ui/T)))/ - (om*T) + ((2.d0*k**2/dcosh(ui/T)**4)*(1.d0/dcosh(T*dtanh(ui/T))**2)*dtanh(T*dtanh(ui/T))**2)/ - om - ((k**4/dcosh(ui/T)**4)*(1.d0/dcosh(T*dtanh(ui/T))**4)*dtanh(T*dtanh(ui/T))**2)/om**3 c This is the version from Aug. 2013 and there is a mistake in the last term. c W = Omega + 3.d0*omp**2/(8.d0*om**3) - om2p/(4.d0*om**3) W = Omega + 3.d0*omp**2/(8.d0*om**3) - om2p/(4.d0*om**2) c W = om 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) c WRITE(12,35) k, T, ui, Y(1), Y(2),Y(3),Y(4) 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. 3 agrees. dif1 = (fr-fr0)/fr dif2 = (fpr-fpr0) dif3 = (fr-fr0) dif4 = (fpi-fpi0)/fpi 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) k,T,ui,dif1,dif2,dif4 c Write(14,35) k,T,ui,fr,fpr,fi,fpi c Note I want to see what happens to beta2 as u increases Write(15,*) k,ui,beta2 DO IT = 1,nu u = ui+delu CALL ODEINT(Y,4,ui,u,EPS,H1,HMIN,NOK,NBAD,DERIVS,BSSTEP) c WRITE(12,*) k, u, Y(1), Y(2), Y(3), Y(4) a = dcosh(T*dtanh(u/T)) h = (1.d0/dcosh(u/T)**2)*dtanh(T*dtanh(u/T)) hp = (1.d0/dcosh(u/T)**4)*(1.d0/dcosh(T*dtanh(u/T))**2) - - ((2.d0/dcosh(u/T)**2)*dtanh(u/T)*dtanh(T*dtanh(u/T)))/T 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(u/T)**2)*(1.d0/dcosh(T*dtanh(u/T))**2)*dtanh(T*dtanh(u/T)))/om) om2p = -(((k**2/dcosh(u/T)**4)*(1.d0/dcosh(T*Tanh(u/T))**4))/om) + - ((2.d0*k**2/dcosh(u/T)**2)*(1.d0/dcosh(T*dtanh(u/T))**2)*dtanh(u/T)*dtanh(T*dtanh(u/T)))/ - (om*T) + ((2.d0*k**2/dcosh(u/T)**4)*(1.d0/dcosh(T*dtanh(u/T))**2)*dtanh(T*dtanh(u/T))**2)/ - om - ((k**4/dcosh(u/T)**4)*(1.d0/dcosh(T*dtanh(u/T))**4)*dtanh(T*dtanh(u/T))**2)/om**3 c This is the original incorrect one. c W = Omega + 3.d0*omp**2/(8.d0*om**3) - om2p/(4.d0*om**3) W = Omega + 3.d0*omp**2/(8.d0*om**3) - om2p/(4.d0*om**2) c W = om 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(15,35) k,u,om,Omega,W,V,beta2 Write(15,*) k,u,beta2 c Write(14,35) k,T,u,fr,fpr,fi,fpi ui = u 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 om0 = dsqrt(k**2/a0**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 Write(13,35) k,u,B2 Enddo Enddo 35 Format(1x,10(d12.5,1x)) STOP END ! -------------------------------------------------------------- ! -------------- SUBROUTINES ----------------------------------- ! -------------------------------------------------------------- SUBROUTINE DERIVS(u,Y,DYDT,NVAR) INTEGER nvar REAL*8 u, T, Y(4), DYDT(4), k, m, a, h, hp, Omega2 COMMON/RealP/k COMMON/IntP/m,T a = dcosh(T*dtanh(u/T)) h = (1.d0/dcosh(u/T)**2)*dtanh(T*dtanh(u/T)) hp = (1.d0/dcosh(u/T)**4)*(1.d0/dcosh(T*dtanh(u/T))**2) - - ((2.d0/dcosh(u/T)**2)*dtanh(u/T)*dtanh(T*dtanh(u/T)))/T 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