c the energy density in a de Sitter background using the c WKB expansion to obtain starting values for the modes c c NOTE THAT EVERYTHING IN THIS PROGRAM IS FOR CONFORMAL COUPLING!!! c IMPLICIT None character junk Integer j,i,ip,it,ik,nk,n,nu,nt,ntot,nbad,nok,ntest,nn,kk,ck Integer n2 Real*8 y,dydu,k,u,uin,ui,delu,eps,h1,hmin,wronski,a,rho,adot,T Real*8 W,Omega,drho,rhos,rhos2,Wd Real*8 pi,m,ci,cons,cons2,drho2,rhod,rhos3,conan Real*8 rhoan,H01,H03,G0,Tre,conan2 Real*8 h,tk,ad,addr Real*8 Td,Hu1,Hu3,Gu,a1,a2,a3,a4,Traan Real*8 reldiffr,reldifft,reldiff1 Real*8 m2,G00,Eulergamma,endenan Real*8 dsech,Tran,lamda Real*8 G,H3,f,Aq,Bq,Cq,cfit Real*8 u1,zeta3,ap2,ad1,ad2,ad3,ad4 PARAMETER(NN=60000,nt=15000) DIMENSION k(nt),Y(nn),DYDU(nn),f(nn),cfit(15,3) Dimension wronski(NN),rho(NN,nt),T(NN,nt) Dimension cons(NN,nt),reldiffr(NN,nt),reldifft(NN,nt) Dimension drho(NN,nt),cons2(NN,nt),drho2(NN,nt),rhos(NN) Dimension rhod(NN,nt),rhoan(NN),Traan(NN) Dimension G0(NN),rhos2(NN),rhos3(NN),Td(NN,nt),a4(NN),conan(NN) Dimension Gu(NN),a1(NN),a2(NN),a3(NN),Tre(NN) Dimension conan2(NN),ck(NN) Common/param/m,k,reldiff1,pi,lamda,Eulergamma Common /cinfo/cfit EXTERNAL DERIVS EXTERNAL BSSTEP c modes2.dat is the input file that contain the information you need to c run the code. Yresult2.txt is the output file. initial*.dat is the c initial value of Y(i), Y(i+1), Y(i+2),Y(i+3). pi,eulergamma are c constant. OPEN(15,FILE='modes0m4i1t4.dat') open(16,FILE='SETresult0m4i1t4.dat') open(14,FILE = 'scalefresult0m4i1t4.dat ') open(13,FILE='rhotest0m4i1t4.dat') open(17,FILE='Ttest0m4i1t4.dat') open(12,FILE='ap2fit0m4i1t4.dat') open(18,FILE='ap1fit0m4i1t4.dat') open(19,FILE='afit0m4i1t4.dat') do j = 1,15 Read(19,*) cfit(j,1) enddo do j = 1,15 Read(18,*) cfit(j,2) enddo do j = 1,15 Read(12,*) cfit(j,3) enddo c Open(14,File='initialD4.dat') c Open(13,File='initialDY4.dat') pi = 2.d0*dasin(1.d0) c lamda = 0.5d0 c zeta = 1.202056903159594d0 Eulergamma = 0.5772156649015329d0 c junk allows for comments in the data file c nk is the number of values of k c uin is the initial time, delu is the time step, and nu is the number of times c gamma is the one free parameter in the mode equation when everything is scaled. c Note that gamma2 = gamma^2 = - nu^2 c eps is the desired accuracyof the differential equation solver, h1 is the initial time c step for it and hmin is the minimum size for a time step READ(15,*) junk READ(15,*) nk,uin,delu,nu,tk,reldiff1 READ(15,*) junk Read(15,*) m,lamda Read(15,*) junk READ(15,*) EPS,H1,HMIN !$$$ntest is some parameter in the differential equation solver. This solver !$$$comes from Numerical Recipes. ci = 1.d0/6.d0 ntest = 0 c ci is a constant. ntot is total number of Y. it=1 Y(1)=a(uin) c Y(1) =dsqrt(3.d0/lamda)*dcosh(dsqrt(lamda/3.d0)*uin) a1(it) = a(uin)*ad1(uin) a2(it) = a(uin)*ad1(uin)**2+a(uin)**2*ad2(uin) a3(it) = a(uin)*ad1(uin)**3+4.d0*a(uin)**2*ad1(uin)*ad2(uin) + +a(uin)**3*ad3(uin) a4(it) = a(uin)*ad1(uin)**4+11.d0*a(uin)**2*ad1(uin)**2*ad2(uin) + +4.d0*a(uin)**3*ad2(uin)**2+7.d0*a(uin)**3*ad1(uin)*ad3(uin) + +a(uin)**4*ad4(uin) H01 =-36.d0*a3(it)*a1(it)/a(uin)**6+72.d0*a2(it)*a1(it)**2 / /a(uin)**7+18.d0*a2(it)**2/a(uin)**6+36.d0*a1(it)**2/a(uin)**6 - -18.d0/a(uin)**4 Hu1=-36.d0*a4(it)/a(uin)**5+144.d0*a3(it)*a1(it)/a(uin)**6 - -216.d0*a2(it)*a1(it)**2/a(uin)**7+108.d0*a2(it)**2/a(uin)**6 + +72.d0*a2(it)/a(uin)**5-72.d0*a1(it)**2/a(uin)**6 c Note that Y(1) = Real f, Y(2) = Real fp, Y(3) = Im f, c Y(4) = Im fp. ntot = nk*4+2 n = 0 Print*,'ntot = ',ntot Do i = 3,ntot-3,4 n = n+1 k(n) = n Omega = dsqrt(k(n)**2 + m**2*Y(1)**2) c this is the initial condition for second order wkb W = Omega c Wd = -(k(n)**2-0.25d0)*adot(uin)/(W*a(uin)**3) c The only way to get beta = 0 at the initial time is to use 1/W rather c than Wm1 c This is so I can use them later f(i) = dsqrt(1.d0/(2.d0*W)) f(i+1) = 0.d0 c -Wd/(2.d0*W)**1.5d0 f(i+2) = 0.d0 f(i+3) = -dsqrt(W/2.D0) c Read(14,*) Y(i),Y(i+2) c Read(13,*) Y(i+1),Y(i+3) c this part is used to check wronski condition c wronski(i) = Y(i+1)*Y(i+2)-Y(i+3)*Y(i)-0.5d0 c write(16,*) k(n),wronski(i) ENDdo c DERIVS is a subroutine that steps the equations forward in time. c It is written by me and is used for the ODE solver. 21 format(3(E22.15,2x)) 23 format(I5,2x,4(E22.15,2x)) 22 format(2(I5,2x),E22.15,2x,I5,2x,E22.15) 20 format(I3,2x,I5,2x,3(E22.15E1,2x)) c this part is used to generate the initial energy density n = 0 it = 1 DO i = 3,ntot-3,4 n = n+1 c this is the method in the paper PhysRevD.61.024003 c and PhysRevD.36.2963 rho(it,n) =k(n)**2*((f(i+1)**2+f(i+3)**2) + + (f(i)**2+f(i+2)**2)*(k(n)**2+m**2*Y(1)**2)) / /(4.d0*pi**2*Y(1)**4) rhod(it,n) =(k(n)**3+k(n)*m**2*Y(1)**2/2-m**4*Y(1)**4/(8.d0*k(n))) c - -(adot(u)**2+1)*(3.d0*m**2*a(u)**2)/(2.d0*k(n))) / /(4*pi**2*Y(1)**4) c calculating the relative difference between c the unregularized rho and the subtraction term c to check if we are losing too many digits of accuracy. c it it is, we exit the cycle reldiffr(it,n) = abs((rho(it,n)-rhod(it,n))/rhod(it,n)) if(reldiffr(it,n)