c For now some output files are commentted out to save space
c energy.tst beta.out pressure.tst
c This version does an arbitrarily large number of values of k starting with k = 1
c It does second order WKB according to Emil's approx.
c
c NOTE THAT EVERYTHING IN THIS PROGRAM IS FOR CONFORMAL COUPLING!!!
c
IMPLICIT None
Integer i,ip,it,ik,nk,n,nu,nt,ntot,nbad,nok,ntest,junk,nn,nm,nkvac,nkwkb2
Real*8 y,dydu,k,u,uin,ui,delu,eps,h1,hmin
Real*8 W,Wp,Wm1,V,Omega,beta2,beta2T
Real*8 pi,gamma2
Real*8 T1,T2,T3,T4,dif
Real*8 Rek,Imk,Ref1,h,epN,epR,epI,prN,prR,prI
Real*8 om,WW,VV
Real*8 enden,enden1,enden2,press,press1,press2,ep,p,delep,delep1,delep2,delp,delp1,delp2
Real*8 ep1,ep2,p1,p2,epratio,pratio
Real*8 kk,omk,m2,denvac,densub,endenvac,diff,diffv,H3tt,G00,Eulergamma,endenan
Real*8 summand,dsech,summpr,pressvac,pressan,diffvpr,prvac,prsub,Trsub,Tran
Real*8 diff1,G,H3,p2a,p2b,press2a,press2b,delp2a,delp2b
Real*8 epnum,epnum1,epnum2,prnum,prnum1,rest,depnumdu,u1
PARAMETER(NN=10000,nt=700)
DIMENSION k(nn),Y(NN),DYDU(NN),beta2(nn,nt),Rek(nn,nt),Imk(nn,nt),WW(nn,nt),VV(nn,nt)
Dimension om(nn,nt),epnum1(nn),epnum2(nn),prnum1(nn)
Common/param/gamma2,k
EXTERNAL DERIVS
EXTERNAL BSSTEP
OPEN(15,FILE='modes.dat')
OPEN(16,FILE='k.out')
OPEN(17,FILE='k.tst')
Open(18,file='energy.out')
Open(19,file='energy.tst')
Open(20,file='pressure.tst')
Open(21,file='denvac.tst')
Open(22,file='beta2.out')
Open(23,file='prvac.tst')
Open(24,file='pressure.out')
Open(25,file='test.out')
Open(26,file='mode_conservation.out')
Open(27,file='k1.out')
Open(28,file='k10.out')
Open(29,file='k100.out')
Open(30,file='k1000.out')
pi = 2.d0*dasin(1.d0)
Eulergamma = 0.5772156649015329d0
c junk allows for comments in the data file
c nk is the number of values of k - we won't make use of it in this version!
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
READ(15,*) junk
Read(15,*) gamma2
Read(15,*) junk
READ(15,*) EPS,H1,HMIN
!$$$ntest is some parameter in the differential equation solver. This solver
!$$$comes from Numerical Recipes.
m2 = gamma2 + 0.25d0
ntest = 0
ntot = 4*nk
c Note that Y(1) = Real f, Y(2) = Real fp, Y(3) = Im f,
c Y(4) = Im fp.
ntot = nk*4
n = 0
Print*,'ntot = ',ntot
Do i = 1,ntot-3,4
n = n+1
k(n) = n
om(n,1) = dsqrt(k(n)**2/dcosh(uin)**2 + m2)
Omega = dsqrt((k(n)**2-0.25d0)/dcosh(uin)**2 + gamma2)
W = Omega + (dtanh(uin)**2/(8.d0*om(n,1)))*(1.d0-6.d0*m2/om(n,1)**2
- + 5.d0*m2**2/om(n,1)**4)
- + (1.d0/(4.d0*om(n,1)*dcosh(uin)**2))*(1.d0-m2/om(n,1)**2)
V = dtanh(uin)*(1.d0 - m2/om(n,1)**2)
If(W.le.0.d0) then
Print *,'n = ',n,' W = ',W
Stop
Endif
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
WW(n,1) = W
VV(n,1) = V
Y(i) = dsqrt(1.d0/(2.d0*W))
Y(i+1) = dsqrt(1.d0/(2.d0*W))*(V/2.d0)
Y(i+2) = 0.d0
Y(i+3) = -dsqrt(W/2.D0)
T1 = (Y(i+1)**2+Y(i+3)**2)*(1.d0/(2.d0*W))
T2 = - Y(i+1)*Y(i+2) + Y(i+3)*Y(i)
T3 = - V*(Y(i)*Y(i+1)+Y(i+2)*Y(i+3))*(1.d0/(2.d0*W))
T4 = + (V**2/(8.d0*W) + W/2.d0)*(Y(i)**2+Y(i+2)**2)
beta2T = T1+T2+T3+T4
beta2(n,1) = (1.d0/(2.d0*W))*(Y(i+1)**2+Y(i+3)**2
- - 2.d0*W*Y(i+1)*Y(i+2) + 2.d0*W*Y(i+3)*Y(i)
- - V*(Y(i)*Y(i+1)+Y(i+2)*Y(i+3))
- + (V**2/4.d0 + W**2)*(Y(i)**2+Y(i+2)**2))
If(n.eq.10) then
dif = ((Y(i)**2+Y(i+2)**2)-1.d0/(2.d0*W))/(1.d0/(2.d0*W))
Write(16,27)uin,k(n),Omega,Y(i),Y(i+1),Y(i+2),Y(i+3),beta2(n,1)
Write(17,27)uin,k(n),W,V,dif,T1,T2,T3,T4,beta2T
Endif
c Now print out the initial values of beta2 to be plotted for certain modes
If(n.eq.1) then
Write(27,27) k(1),uin,beta2(n,1)
else if(n.eq.10) then
Write(28,27) k(10),uin,beta2(n,1)
else if(n.eq.100) then
Write(29,27) k(100),uin,beta2(n,1)
else if(n.eq.1000) then
Write(30,27) k(1000),uin,beta2(n,1)
Endif
Ref1 = (Y(i)*Y(i+1)+Y(i+2)*Y(i+3))
Rek(n,1) = -(Y(i+1)**2+Y(i+3)**2)*(1.d0/(2.d0*W))+ (V/(2.d0*W))*Ref1
- +(W/2.d0- V**2/(8.d0*W))*(Y(i)**2+Y(i+2)**2)
Imk(n,1) = Ref1 - 0.5d0*V*(Y(i)**2+Y(i+2)**2)
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.
c Now start the do loop that integrates the equations forward in time.
u = uin
DO it = 1,nu
ui = u
u = u+delu
CALL ODEINT(Y,ntot,ui,u,EPS,H1,HMIN,NOK,NBAD,DERIVS,BSSTEP)
n = 0
Do i = 1,ntot-3,4
n = n+1
om(n,it+1) = dsqrt(k(n)**2/dcosh(u)**2+m2)
Omega = dsqrt((k(n)**2-0.25d0)/dcosh(u)**2 + gamma2)
W = Omega + (dtanh(u)**2/(8.d0*om(n,it+1)))*(1.d0-6.d0*m2/om(n,it+1)**2
- + 5.d0*m2**2/om(n,it+1)**4)
- + (1.d0/(4.d0*om(n,it+1)*dcosh(u)**2))*(1.d0-m2/om(n,it+1)**2)
V = dtanh(u)*(1.d0 - m2/om(n,it+1)**2)
c This is so I can use them later
WW(n,it+1) = W
VV(n,it+1) = V
beta2(n,it+1) = (1.d0/(2.d0*W))*(Y(i+1)**2+Y(i+3)**2
- - 2.d0*W*Y(i+1)*Y(i+2) + 2.d0*W*Y(i+3)*Y(i)
- - V*(Y(i)*Y(i+1)+Y(i+2)*Y(i+3))
- + (V**2/4.d0 + W**2)*(Y(i)**2+Y(i+2)**2))
If(n.eq.10) then
Write(16,27) u,Omega,Y(i),Y(i+1),Y(i+2),Y(i+3),beta2(n,it+1)
T1 = (Y(i+1)**2+Y(i+3)**2)*(1.d0/(2.d0*W))
T2 = - Y(i+1)*Y(i+2) + Y(i+3)*Y(i)
T3 = - V*(Y(i)*Y(i+1)+Y(i+2)*Y(i+3))*(1.d0/(2.d0*W))
T4 = + (V**2/(8.d0*W) + W/2.d0)*(Y(i)**2+Y(i+2)**2)
beta2T = T1+T2+T3+T4
dif = ((Y(i)**2+Y(i+2)**2)-1.d0/(2.d0*W))/(1.d0/(2.d0*W))
Write(17,27)u,k(n),W,V,dif,T1,T2,T3,T4,beta2T
Endif
c Now print out the initial values of beta2 to be plotted for certain modes
If(n.eq.1) then
Write(27,27) k(1),u,beta2(n,it+1)
else if(n.eq.10) then
Write(28,27) k(10),u,beta2(n,it+1)
else if(n.eq.100) then
Write(29,27) k(100),u,beta2(n,it+1)
else if(n.eq.1000) then
Write(30,27) k(1000),u,beta2(n,it+1)
Endif
Ref1 = (Y(i)*Y(i+1)+Y(i+2)*Y(i+3))
Rek(n,it+1) = -(Y(i+1)**2+Y(i+3)**2)*(1.d0/(2.d0*W))+ (V/(2.d0*W))*Ref1
- +(W/2.d0- V**2/(8.d0*W))*(Y(i)**2+Y(i+2)**2)
Imk(n,it+1) = Ref1 - 0.5d0*V*(Y(i)**2+Y(i+2)**2)
Enddo
Enddo
c Now compute the energy density and the pressure
u = uin-delu
Do it= 1,nu+1
u = u + delu
h = dtanh(u)
c Let's print out beta2 for all times.
Do n=1,nk
Write(22,27) u,k(n),beta2(n,it)
Enddo
H3tt = 3.d0
H3 = 12.d0
G00 = -3.d0
G = -12.d0
endenan = H3tt/(2880.d0*pi**2)+m2*G00/(288.d0*pi**2) -(m2**2/(64.d0*pi**2))*
- (0.5d0 + dlog(m2*dcosh(u)**2/4.d0)+2.d0*Eulergamma)
c Note that the sign of the trace is the opposite from the Anderson-Eaker paper.
Tran = -H3/(2880.d0*pi**2) - m2*G/(288.d0*pi**2) + (m2**2/(16.d0*pi**2))*
- (1.d0 + dlog(m2*dcosh(u)**2/4.d0)+2.d0*Eulergamma)
pressan = (endenan + Tran)/3.d0
endenvac = 0.d0
enden1 = 0.d0
enden2 = 0.d0
pressvac = 0.d0
press1 = 0.d0
c press2 = 0.d0
press2a = 0.d0
press2b = 0.d0
Do n = 1,nk
epN = (1.d0/(2.d0*WW(n,it)))*(om(n,it)**2+WW(n,it)**2 + (VV(n,it)-h)**2/4.d0)
epR = (1.d0/(2.d0*WW(n,it)))*(om(n,it)**2-WW(n,it)**2 + (VV(n,it)-h)**2/4.d0)
epI = (VV(n,it)-h)/2.d0
prN = (1.d0/(6.d0*WW(n,it)))*(om(n,it)**2-2.d0*gamma2-0.5d0+WW(n,it)**2
- + (VV(n,it)-h)**2/4.d0)
prR = (1.d0/(6.d0*WW(n,it)))*(om(n,it)**2-2.d0*gamma2-0.5d0-WW(n,it)**2
- + (VV(n,it)-h)**2/4.d0)
prI = (VV(n,it)-h)/6.d0
ep1 = (1.d0/(2.d0*pi**2*dcosh(u)**3))*k(n)**2*epN*beta2(n,it)
ep2 = (1.d0/(2.d0*pi**2*dcosh(u)**3))*k(n)**2*(epR*Rek(n,it) + epI*Imk(n,it))
p1 = (1.d0/(2.d0*pi**2*dcosh(u)**3))*k(n)**2*prN*beta2(n,it)
p2 = (1.d0/(2.d0*pi**2*dcosh(u)**3))*k(n)**2*(prR*Rek(n,it) + prI*Imk(n,it))
p2a = (1.d0/(2.d0*pi**2*dcosh(u)**3))*k(n)**2*prR*Rek(n,it)
p2b = (1.d0/(2.d0*pi**2*dcosh(u)**3))*k(n)**2*prI*Imk(n,it)
c Write(25,27) u,k(n),epR,prR,epI,prI
c Write(25,27) u,k(n),ep1,p1,ep2,p2
Write(25,27) u,k(n),prR,Rek(n,it)
denvac = (1.d0/(4.d0*pi**2*dcosh(u)**3))*k(n)**2*epN
prvac = (1.d0/(4.d0*pi**2*dcosh(u)**3))*k(n)**2*prN
c This is the WKB subtraction term from the Anderson-Eaker paper
densub = (1.d0/(4.d0*pi**2*dcosh(u)**4))*k(n)**2*(k(n) + m2*dcosh(u)**2/(2.d0*k(n))
- - m2**2*dcosh(u)**4/(8.d0*k(n)**3))
c Note that the sign of the trace is the opposite from the Anderson-Eaker paper.
Trsub = -m2*k(n)/(4.d0*pi**2*dcosh(u)**2) + m2**2/(8.d0*pi**2*k(n))
prsub = (densub+Trsub)/3.d0
c Now check conservation
If(it.eq.1) then
epnum = ep1+ep2+denvac-densub
epnum2(n) = epnum
Else if(it.eq.2) then
u1 = u
epnum = ep1+ep2+denvac-densub
epnum1(n) = epnum
prnum = p1+p2+prvac-prsub
prnum1(n) = prnum
Else
epnum = ep1+ep2+denvac-densub
prnum = p1+p2+prvac-prsub
depnumdu = (epnum - epnum2(n))/(2.d0*delu)
rest = -3.d0*dtanh(u1)*(epnum1(n)+prnum1(n))
diff = (depnumdu-rest)/depnumdu
c Write(26,27) u1,k(n),depnumdu,rest,diff
u1 = u
epnum2(n) = epnum1(n)
epnum1(n) = epnum
prnum1(n) = prnum
Endif
diff = dabs((dabs(denvac)-dabs(densub))/(dabs(denvac)+dabs(densub)))
diff1 = dabs((dabs(prvac)-dabs(prsub))/(dabs(prvac)+dabs(prsub)))
If(diff.lt.1.d-15.or.diff1.lt.1.d-15) then
Print*,'u = ',u,' k = ',k(n)
c Print*,'diff energy density = ',diff,' diff pressure = ',diff1
c Print*,'denvac = ',denvac,' densub = ',densub
c Print*,'prvac = ',prvac,' prsub = ',prsub
Go to 150
Endif
endenvac = endenvac + denvac - densub
pressvac = pressvac + prvac - prsub
If(it.eq.nu) Then
summand = denvac-densub
diffv = (denvac-densub)/endenvac
summpr = prvac-prsub
diffvpr = (prvac-prsub)/pressvac
c Write(21,27) u,k(n),diff,denvac,densub,summand,endenvac,diffv
c Write(23,27) u,k(n),diff1,prvac,prsub,summpr,pressvac,diffvpr
Endif
c There are no subtraction terms for these so it makes sense to keep all the modes
150 Continue
enden1 = enden1+ ep1
enden2 = enden2+ ep2
press1 = press1+ p1
c press2 = press2 + p2
press2a = press2a + p2a
press2b = press2b + p2b
c If(it.eq.10) Then
delep1 = ep1/enden1
delep2 = ep2/enden2
delp1 = p1/press1
c delp2 = p2/press2
delp2a = p2a/press2a
delp2b = p2b/press2b
c Write(19,27) u,k(n),ep1,enden1,delep1,ep2,enden2,delep2
c Write(19,27) u,k(n),enden1,delep1,enden2,delep2
c Write(20,27) u,k(n),press1,delp1,press2,delp2
c Write(20,27) u,k(n),p1,press1,delp1,p2a,press2a,delp2a,p2b,press2b,delp2b
c Endif
Enddo
enden = enden1+enden2+endenvac+endenan
epratio = dabs(enden2/(enden1+endenvac+endenan))
press2 = press2a+press2b
press = press1+press2+pressvac+pressan
pratio = dabs(press2/(press1+pressvac+pressan))
c This is the ratio of the last contribution (i.e. at largest k) to the total
delep = (ep1+ep2)/enden
delp = (p1+p2)/press
Write(18,27) u,endenan,endenvac,enden2,enden1,enden
c Write(18,27) u,enden1,enden2,endenvac,endenan,enden,epratio
c Write(18,27) u,delep,press1,press2,press,pratio,delp
Write(24,27) u,pressan,pressvac,press2,press1,press
Enddo
27 format(1x,11(e14.7,1x))
STOP
END
!$$$This is the subroutine called by the differential equation solver to
!$$$step the equations forward in time.
SUBROUTINE DERIVS(u,Y,DYDU,NVAR)
IMPLICIT None
Integer i,ik,nvar,nn,nt,n
Real*8 u,Y, DYDU, k,Omega2,gamma2
PARAMETER(NN=10000,nt=100)
DIMENSION Y(NVAR),DYDU(NVAR),k(nn)
Common/param/gamma2,k
!$$$ Y(N) is Real(f(k)), Y(N+1) is Real(f'(k)), Y(N+2) IS IM(f(k))
!$$$ Y(N+3) IS IMP(f(k)) AND SIMILARLY FOR THE DYDT'S.
n = 0
DO i=1,NVAR-3,4
n = n+1
Omega2 = (k(n)**2-0.25d0)/dcosh(u)**2 + gamma2
DYDU(i) = Y(i+1)
DYDU(i+1) = - Omega2*Y(i)
DYDU(i+2) = Y(i+3)
DYDU(i+3) = - Omega2*Y(i+2)
end do
RETURN
END
Double Precision Function dsech(u)
Real*8 u
dsech = 1.d0/dcosh(u)
Return
End