c This version uses the profile for the electric field with parameters t2 and b. c Note in the paper the notation t1 is used instead of t2. c IMPLICIT NONE Integer j,ik, nt, JUNK, IT, nbad,nok,nk,itest Real*8 t, ti,tin,delt, Y, m, gamma2,k Real*8 eps,h1,hmin Real*8 Omega,om,omp,om2p,W,V,fr,fpr,fi,fpi,beta2 Real*8 om0,fr0,fpr0,fi0,fpi0,dif1,dif2,dif3,dif4 Real*8 for,fopr,foi,fopi,B2 Real*8 pi,delumin,deluj,test,ttest Real*8 htest,mak Real*8 t2,b,lam Real*8 TwoWm1,f02 Real*8 sech,Are,Aim,Bre,Bim Real*16 kk,ki,delk DIMENSION Y(4) COMMON/RealP/k COMMON/IntP/lam,t2,b External DERIVS, BSSTEP OPEN(UNIT = 11, FILE = 'modes_efield.dat') Open(unit=12,file = 'B2.out') Open(unit=13,file='modes.tst') Open(unit=14,file='beta2.out') READ(11,*) JUNK READ(11,*) ki, delk, nk, tin, delt, nt READ(11,*) JUNK Read(11,*) lam,t2,b READ(11,*) JUNK Read(11,*) eps,h1,hmin Write(13,*) b,t2,ti pi = 2.d0*dasin(1.d0) c I need the quadruple precision to get enough accuracy but then double c precision in k to do the calculations kk = ki-delk Do ik = 1,nk kk = kk+delk k = kk ti = tin t = ti 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-t2/2.d0)**2 + lam) fr0 = 1.d0/dsqrt(2.d0*om0) fpr0 = 0.d0 fi0 = 0.d0 fpi0 = -dsqrt(om0/2.d0) om = dsqrt(lam + (k + dlog(dcosh(b*(t + t2))*sech(b*(t - t2)))/(4.d0*b))**2) omp = -((4*b*k + dlog(dcosh(b*(t + t2))*sech(b*(t - t2))))* - (dtanh(b*(t - t2)) - dtanh(b*(t + t2))))/ - (4.d0*b*dsqrt((16*b**2*(k**2 + lam) + - 8*b*k*dlog(dcosh(b*(t + t2))*sech(b*(t - t2))) + - dlog(dcosh(b*(t + t2))*sech(b*(t - t2)))**2)/b**2)) om2p = -(sech(b*(t - t2))**2*sech(b*(t + t2))**2*dsinh(2*b*t2)* - ((4*b*k + dlog(dcosh(b*(t + t2))*sech(b*(t - t2))))* - (16*b**2*(k**2 + lam) + - 8*b*k*dlog(dcosh(b*(t + t2))*sech(b*(t - t2))) + - dlog(dcosh(b*(t + t2))*sech(b*(t - t2)))**2)*dsinh(2*b*t) - - 16*b**2*lam*dsinh(2*b*t2)))/ - (4.d0*b**2*((16*b**2*(k**2 + lam) + - 8*b*k*dlog(dcosh(b*(t + t2))*sech(b*(t - t2))) + - dlog(dcosh(b*(t + t2))*sech(b*(t - t2)))**2)/b**2)**1.5d0) c Write(13,*) b,t2,ti c Write(13,*) om,omp,om2p W = om + 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 c Write(13,*) fr0,fpr0,fi0,fpi0 c Write(13,*) fr,fpr,fi,fpi c Write(13,*) k,t2,b,dif1,dif2,dif3,dif4 Do j = 1,nt t = ti + delt CALL ODEINT(Y,4,ti,t,EPS,H1,HMIN,NOK,NBAD,DERIVS,BSSTEP) ti = t c Enddo c The following is to compute |beta|^2 om = dsqrt(lam + (k + dlog(dcosh(b*(t + t2))*sech(b*(t - t2)))/(4.d0*b))**2) omp = -((4*b*k + dlog(dcosh(b*(t + t2))*sech(b*(t - t2))))* - (dtanh(b*(t - t2)) - dtanh(b*(t + t2))))/ - (4.d0*b*dsqrt((16*b**2*(k**2 + lam) + - 8*b*k*dlog(dcosh(b*(t + t2))*sech(b*(t - t2))) + - dlog(dcosh(b*(t + t2))*sech(b*(t - t2)))**2)/b**2)) om2p = -(sech(b*(t - t2))**2*sech(b*(t + t2))**2*dsinh(2*b*t2)* - ((4*b*k + dlog(dcosh(b*(t + t2))*sech(b*(t - t2))))* - (16*b**2*(k**2 + lam) + - 8*b*k*dlog(dcosh(b*(t + t2))*sech(b*(t - t2))) + - dlog(dcosh(b*(t + t2))*sech(b*(t - t2)))**2)*dsinh(2*b*t) - - 16*b**2*lam*dsinh(2*b*t2)))/ - (4.d0*b**2*((16*b**2*(k**2 + lam) + - 8*b*k*dlog(dcosh(b*(t + t2))*sech(b*(t - t2))) + - dlog(dcosh(b*(t + t2))*sech(b*(t - t2)))**2)/b**2)**1.5d0) W = om + 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 If(ik.eq.1) then Write(14,*) k,t,beta2 Endif Enddo 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 c If(it.eq.nt) then c Write(14,35) fr,fpr,fi,fpi c fr = Y(1) c fpr = Y(2) c fi = Y(3) c fpi = Y(4) c Write(14,35) fr,fpr,fi,fpi c Note that in the out region it is +t2 not -t2 as it is for the in region om0 = dsqrt((k+t2/2.d0)**2 + lam) for = 1.d0/dsqrt(2.d0*om0) fopr = 0.d0 foi = 0.d0 fopi = -om0/dsqrt(2.d0*om0) 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 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,*) t2,t,k,Are,Aim,Bre,Bim,B2 c Endif 100 Continue Enddo 35 Format(1x,10(d14.7,1x)) STOP END ! -------------------------------------------------------------- ! -------------- SUBROUTINES ----------------------------------- ! -------------------------------------------------------------- SUBROUTINE DERIVS(t,Y,DYDT,NVAR) INTEGER nvar REAL*8 t, t2, b,Y(4), DYDT(4), lam, om2,k COMMON/RealP/k COMMON/IntP/lam,t2,b om2 = (k + (1.d0/(4.d0*b))*dlog(dcosh(b*(t + t2))/dcosh(b*(t - t2))))**2 + lam c Print*,t,om2 DYDT(1) = Y(2) DYDT(2) = -om2*Y(1) DYDT(3) = Y(4) DYDT(4) = -om2*Y(3) RETURN END ! ! Double Precision Function sech(x) Implicit None Real*8 x sech = 1.d0/dcosh(x) Return End