c The program now prints out everything in rfs.out to 14 digits.
c This program has commented out the print statement put in to print out the values of
c the modes for all values of k and for all values of t at which they
c are computed.
c This version has g -> g^2, i.e. the interaction has been changed
c from g phi^2 psi^2 to g^2 phi^2 psi^2
c
c To be careful I have everywhere changed g to gnew
c
!$$$c This program computes the backreaction of free scalar fields
!$$$c on an inflaton field in flat space in the large N limit which
!$$$c is equivalent to the one loop backreaction with just one
!$$$c free quantized field present. Quantum fluctuations of the
!$$$c inflaton field are ignored in this approximation.
!$$$c
!$$$c This version is for a flat space background. It is easily
!$$$c generalized to other expansions. A version to compute the
!$$$c backreaction on the spacetime geometry will come later.
!$$$c
!$$$c The self-coupling constant for the inflaton is lamphi and the
!$$$c coupling between it and the quantum field is g^2. If
!$$$c g is set to zero then the phi field obeys a classical
!$$$c wave equation. Also is zero in this case if the field
!$$$c is in the vacuum state which is what the program is assuming..
!$$$c
!$$$ The basic structure of the program is as follows:
!$$$ 1. Initial values and values for various parameters are read in.
!$$$ 2. A subroutine is called to compute some constants relating to
!$$$ intervals used in the numerical integrations of the modes.
!$$$ Simpson's rule is being used.
!$$$ 3. The initial values for the modes of the quantized fields are
!$$$ computed using a second order adiabatic approximation. This
!$$$ needs to be changed to fourth order when the backreaction on
!$$$ the spacetime geometry is included.
!$$$ 4. The differential equation solver is then called. It integrates
!$$$ the coupled equations for the inflaton and mode functions. It
!$$$ has also been set up to integrate the equations for the scale
!$$$ factor. Right now this is set up to be trivial so that the space
!$$$ is flat space. However I wrote it this way so that it would be
!$$$ easier to generalize the program to include the backreaction on
!$$$ the geometry.
IMPLICIT None
Integer nn,nm,n,i,ip,nk1,nk2,nk3,nt,ntot,nbad,it,nok,ntest,junk
Integer iprint,iprt,nprt,istate,iback,ifile,iw
Integer ii,iii,nmode
Real*8 y,dydt,c1,c2,c3,delk1,delk2,delk3,k,ti,t,delt,eps,h1,hmin
Real*8 mu,mphi,xi,phi,phip,phipp,a,ap,W,Wp,lamphi,gnew,gphi2
Real*8 pi,lamk,psi2,psi2a,rhophi,rhopsi,rhopsia,rhotot
Real*8 kadiabatic,omega,wm1,wterm,testphipp
Real*8 beta1,beta2,beta3,Lam,psi2cutoff,psi2exact,psi2lamdif,psi2modedif
Real*8 psi2dif,psi21,psi22,psi23,dif1,dif2,dif3
Real*8 epsilon,psi21exact,psi22exact,psi23exact,T1,T2,T3,T4,T5,psi2int
PARAMETER(NN=100000,nm=24000)
DIMENSION Y(NN),DYDT(NN),c1(nm),c2(nm),c3(nm),k(3*nm)
common/param/pi,mu,mphi,xi,lamphi,gnew
Common/in1/lamk,c1,c2,c3,k
Common/in2/nk1,nk2,nk3
Common/der/psi2,psi2a,rhophi,rhopsi,rhopsia,psi21,psi22,psi23,psi2int
Common/back/iback
Common/fil/ifile
EXTERNAL DERIVS
EXTERNAL BSSTEP
!$$$
!$$$rfs.dat is the file for all the values of parameters to be read into the
!$$$program. rfs.tst is a file that includes diagnostic output to determine
!$$$how accurately the mode integrations are being done. rfs.out is the main
!$$$output file.
!$$$
OPEN(15,FILE='rfs.dat')
OPEN(16,FILE='rfsint.tst')
OPEN(17,FILE='rfs.out')
c OPEN(18,FILE='rfsmodesum.tst')
OPEN(19,FILE='rfstmn.out')
c This file is for the modes which will be followed in time.
c Open(20,file='modes.out')
!$$$reading the variable junk is just a trick to allow for the labeling of the
!$$$data file.
!$$$
!$$$mphi is the mass of the inflaton, xi is the coupling of the inflaton to
!$$$the scalar curvature, lamphi is the self-coupling of the inflaton, gnew
!$$$is the coupling between the inflaton and the quantum fields.
!$$$
!$$$The program integrates the mode equation in a two step process. First it
!$$$goes from k = 0 to a cutoff lamk. In this stage only the terms in the
!$$$WKBdiv approximation that are of order k and k^3 are subtracted off. In
!$$$the second stage the integral goes from lamk to an upper limit cutoff.
!$$$Here the terms in WKBdiv that go like 1/k are also subtracted off. This
!$$$is consistent with the way in which we define the analytic approximation.
!$$$Here what is happening is the WKBdiv terms are being subtracted off and
!$$$then added back on as part of the analytic approximation. So nk1 is the
!$$$number of modes in the first integration and nk2 the number in the second
!$$$integration. However there is an overlap at k = lamk since that is the
!$$$end and beginning of the two integration ranges. delk1 and delk2 are the
!$$$intervals for the two integration regions.
!$$$
!$$$phi is the value of inflaton field, phip the value of its first derivative.
!$$$The quantity "a" is the scale factor. In this version of course it makes
!$$$sense to just set a = 1.
!$$$
!$$$eps1 is the target accuracy for the differential equation solver for the
!$$$quantities a, ap, phi, phip. These can be less accurate usually than the
!$$$modes. eps2 is the target accuracy for the modes. h1 is an initial value
!$$$for the differential equation solver to try and hmin is the minimum step
!$$$size. The program stops if the step size has to go below this value.
c iprt is the spacing between times when you print out the modes spectrum.
READ(15,*) junk
READ(15,*) mphi,xi,lamphi,gnew
READ(15,*) junk
READ(15,*) DELk1,nk1,delk2,nk2,delk3,nk3,DELT,NT
READ(15,*) junk
c READ(15,*) phi,phip,a
READ(15,*) gphi2,phip,a
phi = dsqrt(gphi2/gnew**2)
READ(15,*) junk
c READ(15,*) EPS1,EPS2,H1,HMIN
READ(15,*) EPS,H1,HMIN
Read(15,*) junk
c istate is 2 for second order and 4 for 4th order.
c iback is 0 for no backreaction on the phi field and 1 for backreaction
c Note that even when there is no backreaction I will use the same 4th order adiabiabatic
c state as if there was so that the comparison can be direct.
Read(15,*) nprt,istate,iback
c This is the n for one of the modes to be printed out as a function of time
Read(15,*) junk
Read(15,*) nmode
Write(17,*) ' 0 gnew = ',gnew,' (g*phi)^2 = ',gphi2
If(iback.eq.0) then
If(istate.eq.2) then
Write(17,*) '0 Background field with adiabatic order 2 states'
Else
Write(17,*) '0 Background field with adiabatic order 4 states'
Endif
Else
If(istate.eq.2) then
Write(17,*) ' 0 Backreaction with adiabatic order 2 states'
Else
Write(17,*) ' 0 Backreaction with adiabatic order 4 states'
Endif
Endif
pi = 2.d0*dasin(1.d0)
Print *,'pi = ',pi
!$$$This is the lower limit cutoff for the integrals that go like 1/k in
!$$$the analytic approximation.
lamk = nk1*delk1
Print *,'lower limit cutoff = ',lamk
!$$$c We'll set the parameter mu in the integrals of the analytic approximation
!$$$c to be equal to 1.
mu = 1.d0
call Intconst(delk1,delk2,delk3)
!$$$This is just for debugging purposes.
Print *,'c1(1) = ',c1(1),'c1(2) = ',c1(2),' c1(nk1) = ',c1(nk1)
Print *,'c2(1) = ',c2(1),'c2(2) = ',c2(2),' c2(nk2) = ',c2(nk2)
Print *,'c3(1) = ',c3(1),'c3(2) = ',c3(2),' c3(nk3) = ',c3(nk3)
Print *,'k(1) = ',k(1),' k(nk1) = ',k(nk1),' k(nk1+1) = ',k(nk1+1)
Print *,'k(nk1+nk2) = ',k(nk1+nk2),'k(nk1+nk2+1) = ',k(nk1+nk2+1),' k(nk1+nk2+nk3) = ',k(nk1+nk2+nk3)
!$$$ntest is some parameter in the differential equation solver. This solver
!$$$comes from Numerical Recipes.
ntest = 0
!$$$c There are 4 variables that we start with. Y(1) = a, Y(2) = ap,
!$$$c Y(3) = phi, Y(4) = phip.
!$$$c In the sums over modes there is an extra value of n that one must
!$$$c be careful of because there is an integral from k = 0 to lamk and
!$$$c then one from lamk to the upper limit cutoff for k. Now the value
!$$$c of k at the lower limit cutoff is the same for both integrals but
!$$$c the integration factors c1 and c2 are different.
!$$$c
!$$$c There are 4 modes for each k; the real part, the imaginary part
!$$$c and their derivatives.
!$$$c
!$$$c Thus the total number of variables for the differential equation
!$$$c solver is 4 + 4*nk1 + 4*(nk2 - 1)+4*(nk3-1) again because of the redundancy
!$$$c in the value of k. Thus
!$$$c
NTOT = 4*nk1 + 4*nk2 + 4*nk3
Y(1) = a
!$$$Let's work in a flat space background.
ap = 0.d0
Y(2) = ap
Y(3) = phi
Y(4) = phip
!$$$c
!$$$c Note that Y(5) = Real psi, Y(6) = Real psip, Y(7) = Im psi,
!$$$c Y(8) = Im psip.
!$$$ We have the same value of k at n = nk1 and n = nk1+1.
!$$$c This is because we want different values for c1(nk1)
!$$$c and c2(nk1+1) even though the value of k is the same
!$$$c for these since they are at the end of the first integration
!$$$c interval and the beginning of the second respectively.
c We can have a problem with initial values for the modes if g^2 is too small as
c it can be for small g^2 phi_0^2. This checks for that and then sets iw = 1 so that
c later on the second order WKB approximation can be used for small values of k.
c
c Now we keep the value of phipp the same for simplicity. It is too complicated to do anything else.
c This is not unreasonable so long as the correction to phipp is small and that is checked
c farther below. Here are some reasons why
c i) you have a UV cutoff anyway.
c ii) you have numerical error anyway
c
c
c This is for a fourth order adiabatic state and the value of phipp is determined
c self-consistently unless the starting value for Wm1 < 0 for one of the k's used.
c
c I do it this way to eliminate infrared divergences.
c
If(istate.eq.4) then
phipp = 1.d0/(1.d0+gnew**2/(48.d0*pi**2))*(-phi + gnew**4*phi**3/(16.d0*pi**2)
- - gnew**4*phi**3*dlog(gnew**2*phi**2/mu**2)/(16.*pi**2))
c Now test to see if Wm1 is going to be positive. If it is for the smallest k used then there is
c no problem.
Omega = dsqrt(k(1)**2 + gnew**2*a**2*phi**2)
Wm1 = 1.d0/Omega + gnew**2*(phip**2+phi*phipp)/(4.*Omega**5)
If(Wm1.lt.0) then
iw = 1
Else
iw = 0
Endif
Print*,'If iw = 0 then 4th order adiabatic was used everywhere which is what we want. iw = ',iw
c Now see if at nk1 it is OK.
Omega = dsqrt(k(nk1)**2 + gnew**2*a**2*phi**2)
Wm1 = 1.d0/Omega + gnew**2*(phip**2+phi*phipp)/(4.d0*Omega**5)
If(Wm1.le.0.d0) then
Print *,' Wm1 < 0 even at k(nk1)'
Stop
Endif
Else
phipp = -phi
Endif
!$$$ The Do loop sets the intial values for the modes.
c If istate = 2 it is 2nd order adiabatic and if istate = 4 it is fourth order
c In both cases it is an infrared finite state.
n = 0
DO i=5,NTOT-3,4
If(n.ne.nk1) then
n = n+1
else
n = n+2
endif
Omega = dsqrt(k(n)**2 + gnew**2*a**2*phi**2)
If(istate.eq.2) then
W = Omega
Wm1 = 1.d0/Omega
else
c Note to get the normalization correct I will use Wm1 everywhere so W = 1/Wm1
c This is useful because I used Wm1 to fourth adiabtic order to compute phipp
c Also as you can see in the commented statement, if I compute W directly I will
c sometimes get a negative number because of the minus sign. There is no such
c problem with Wm1 at large k but there can be at small k.
c W = Omega - gnew**2*(phip**2+phi*phipp)/(4.d0*Omega**3)
Wm1 = 1.d0/Omega + gnew**2*(phip**2+phi*phipp)/(4.d0*Omega**5)
c Use second order adiabatic up to and including k(nk1). There will be a sharp transition but
c there already is one because of the difference in subtraction terms on either side this value of k
c And of course phipp will be wrong but hopefully still close to the correct value
If (iw.eq.1.and.n.le.nk1) then
Wm1 = 1.d0/Omega
Endif
W = 1.d0/Wm1
Endif
c This value of Wp is 3rd order but in flat space with phip=0 as a starting value it is zero
c anyway.
Wp = gnew**2/(Omega)*(a**2*phi*phip + a*ap*phi**2)
Y(i) = dsqrt(Wm1/2.D0)
Y(i+1) = -Wp*(Wm1/2.D0)**1.5D0
Y(i+2) = 0.D0
Y(i+3) = -(W/2.D0)**.5D0
c I want to write the initial values of the modes as well
c t = 0 so it is easiest just to set it to that here
t = 0.d0
c Write(20,*) t,k(n),Y(i),Y(i+1),Y(i+2),Y(i+3)
c We need to print this out twice if n = nk1 because the do loop skips n = nk1+1 but it is easier
c to read that value in the correlation function routine.
c If(n.eq.nk1) then
c Write(20,*) t,k(n),Y(i),Y(i+1),Y(i+2),Y(i+3)
c Endif
end do
!$$$c This is to print out the initial values of app and phipp as well
!$$$c as the test integral and psi2 and psi2a. The latter are printed
!$$$c out through a command in the test subroutine.
TI = 0.d0
T = TI
!$$$DERIVS is a subroutine that steps the equations forward in time.
!$$$It is written by me and is used for the ODE solver. Here it is
!$$$being used to get initial values for certain quantities.
!$$$TEST is a subroutine that integrates a simple function that has the same
!$$$type of k behavior as does the actual mode integral for the quantity
!$$$<\psi>^2. It is therefore useful for determining whether the mode integration
!$$$is likely to be accurate or not, and by how much it might deviate
!$$$in its accuracty.
CALL DERIVS(T,Y,DYDT,NTOT)
WRITE(16,*) ' Time dif int 1 dif int 2 dif int 3 diff total'
CALL TEST(T,Y,ntot)
c This prints out the integrand for the mode sum. This is useful
c if gnew = 0 since then the integrand should always vanish.
ifile = 0
CALL TestModeSum(T,Y,DYDT,NTOT)
rhotot = rhophi+rhopsi
c This is just for debugging purposes
CALL DERIVS(T,Y,DYDT,NTOT)
c
c We can also compare the initial value of with its correct value. There is an
c error from both the cutoff and the density of modes. I'll compute them separately.
c Here psi2cutoff is the exact value we should get with with cutoff. psi2exact is the exact
c value we should get if the cutoff was infinity.
Lam = k(nk1+nk2+nk3)
c This is for the first integration interval which has no 1/k subtraction term.
epsilon = k(nk1)
psi21exact = (1.d0/(4.d0))*(epsilon*dsqrt(epsilon**2+gphi2)-epsilon**2
1 +gphi2*(-dlog((epsilon+dsqrt(epsilon**2+gphi2))/dsqrt(gphi2)))
2 + (phipp/(6.d0*phi))*epsilon**3/(epsilon**2+gphi2)**1.5d0)
dif1 = (psi21exact-psi21)/psi21exact
c This is for the second integration interval
T1 = 0.5d0*(k(nk1+nk2)*dsqrt(k(nk1+nk2)**2+gphi2)
1 - gphi2*dlog((k(nk1+nk2)+dsqrt(k(nk1+nk2)**2+gphi2))))
T2 = -0.5d0*(k(nk1)*dsqrt(k(nk1)**2+gphi2) - gphi2*dlog((k(nk1)+dsqrt(k(nk1)**2+gphi2))))
T3 = -0.5d0*(k(nk1+nk2)**2-k(nk1)**2)
T4 = dlog(k(nk1+nk2)/k(nk1))
T5 = (1.d0/(3.d0*gphi2))*(k(nk1+nk2)**3/(k(nk1+nk2)**2+gphi2)**1.5d0
1 - k(nk1)**3/(k(nk1)**2+gphi2)**1.5d0)
psi22exact = (0.5d0*T1+0.5d0*T2+0.5d0*T3+0.25d0*gphi2*T4+(gphi2*phipp/(8.d0*phi))*T5)
dif2 = (psi22exact-psi22)/psi22exact
T1 = 0.5d0*(Lam*dsqrt(Lam**2+gphi2)
1 - gphi2*dlog((Lam+dsqrt(Lam**2+gphi2))))
T2 = -0.5d0*(k(nk1+nk2)*dsqrt(k(nk1+nk2)**2+gphi2) - gphi2*dlog((k(nk1+nk2)+dsqrt(k(nk1+nk2)**2+gphi2))))
T3 = -0.5d0*(Lam**2-k(nk1+nk2)**2)
T4 = dlog(Lam/k(nk1+nk2))
T5 = (1.d0/(3.d0*gphi2))*(Lam**3/(Lam**2+gphi2)**1.5d0
1 - k(nk1+nk2)**3/(k(nk1+nk2)**2+gphi2)**1.5d0)
psi23exact = (0.5d0*T1+0.5d0*T2+0.5d0*T3+0.25d0*gphi2*T4+(gphi2*phipp/(8.d0*phi))*T5)
dif3 = (psi23exact-psi23)/psi23exact
Print *,'psi21 = ',psi21,' psi21exact = ',psi21exact,' dif1 = ',dif1
Print *,'psi22 = ',psi22,' psi22exact = ',psi22exact,' dif2 = ',dif2
Print *,'psi23 = ',psi23,' psi23exact = ',psi23exact,' dif3 = ',dif3
psi2cutoff = (1.d0/(8.d0*Pi**2))*(Lam*dsqrt(Lam**2+gphi2)-Lam**2+gphi2*(-dlog((Lam+dsqrt(Lam**2+gphi2))/Lam)
1 +dlog(dsqrt(gphi2)*2.d0/mu)) + (phipp/(6.d0*phi))*Lam**3/(Lam**2+gphi2)**1.5d0 - gphi2 )
Print*,' Original value of psi2cutoff = ',psi2cutoff
If(iback.eq.1) then
psi2cutoff = psi2cutoff + (1.d0/(4.d0*pi**2*k(nk1+nk2+nk3)**2))*(3.d0*gnew**4*phi**4/16.d0
1 + gnew**2*phip**2/8.d0 + gnew**2*phi*phipp/8.d0)
Print*,' Value of psi2cutoff when 4th order WKB approximation is added on = ',psi2cutoff
Endif
psi2exact = phipp/(48.d0*pi**2*phi) + (gphi2/(16.d0*pi**2))*(-1.d0 + dlog(gphi2/mu**2))
psi2dif = (psi2exact-psi2)/psi2exact
psi2lamdif = (psi2exact-psi2cutoff)/psi2exact
psi2modedif = (psi2cutoff-psi2)/psi2cutoff
Print*,'psi2 = ',psi2,'psi2a = ',psi2a
Print*,'psi2exact = ',psi2exact,' psi2cutoff = ',psi2cutoff
Print*,' Total relative difference =',psi2dif
Print*,' Cutoff = ',Lam
Print*, ' The relative error due to the existence of a cutoff = ', psi2lamdif
Print*,' The relative error in the mode integral = ',psi2modedif
Write(17,*),'0 Cutoff = ',Lam
Write(17,*) ' 0 Relative error due to the cutoff = ', psi2lamdif
Write(17,*),'0 Relative error for mode integral = ',psi2modedif
c We can check the accuracy of the mode sum by comparing the initial value of phipp
c above which is computed exactly analytically with DYDT(4) which is the numerical
c calculation of the initial value of phipp.
c This only makes sense to do for a 4th order adiabatic state when backreaction is being
c taken into account
If(istate.eq.4.and.iback.eq.1) then
testphipp = (phipp-DYDT(4))/phipp
Print *, 'Initial phipp analytic = ',phipp,' Initial phipp numeric = ',DYDT(4),
- ' percent difference = ', testphipp
Print *, ' If iw = 1 the difference needs to be small since phipp is wrong. iw = ',iw
Write(17,*) '0 Initial phipp analytic = ',phipp
Write(17,*) '0 Initial phipp numeric = ',DYDT(4)
Write(17,*) '0 Percent difference = ', testphipp
Endif
WRITE(17,*) ' 0 Time phi phip phipp psi2 psi2a'
WRITE(19,*) ' Time rhophi rhopsi rhopsia rho_total'
c beta2 = |beta|^2 = Carmen and Emil's definition of produced particles
WRITE(17,27) T,Y(3),Y(4),DYDT(4),psi2,psi2a
WRITE(19,27) T,rhophi,rhopsi,rhopsia,rhotot
c
!$$$
!$$$Now start the do loop that integrates the equations forward in time.
!$$$Ti is the starting time interval for the solver and t is the ending
!$$$interval. Delt is the spacing in time between intervals. nt is the
!$$$number of intervals.
ip = 0
iprint = 1
ifile = 1
iprt = NT/nprt
DO IT = 1,NT
ip = ip+1
print *,ip
TI = T
T = T+DELT
c CALL ODEINT(Y,NTOT,TI,T,EPS1,EPS2,H1,HMIN,NOK,NBAD,DERIVS,BSSTEP)
CALL ODEINT(Y,NTOT,TI,T,EPS,H1,HMIN,NOK,NBAD,DERIVS,BSSTEP)
CALL DERIVS(T,Y,DYDT,NTOT)
CALL TEST(T,Y,ntot)
c This prints out the integrand for the mode sum. This is useful
c if gnew = 0 since then the integrand should always vanish.
If(iprint.eq.iprt) then
CALL TestModeSum(T,Y,DYDT,NTOT)
iprint = 1
else
iprint = iprint+1
endif
rhotot = rhophi+rhopsi
WRITE(17,27) T,Y(3),Y(4),DYDT(4),psi2,psi2a
WRITE(19,27) T,rhophi,rhopsi,rhopsia,rhotot
c beta1 = (1.d0/(2.D0*k(1)))*(Y(6)**2+Y(8)**2)+(-Y(6)*Y(7)+Y(8)*Y(5))
c - +(k(1)/2.d0)*(Y(1)**2+Y(3)**2)
c If(nmode.le.nk1) then
c ii = 4*(nmode)+1
c Else
c ii = 4*(nmode-1)+1
c Endif
c beta2 = (1.d0/(2.D0*k(nmode)))*(Y(ii+1)**2+Y(ii+3)**2)+(-Y(ii+1)*Y(ii+2)+Y(ii+3)*Y(ii))
c - +(k(nmode)/2.d0)*(Y(ii)**2+Y(ii+2)**2)
c iii = 4*(nk1+nk2+nk3-1)+1
c beta3 = (1.d0/(2.D0*k(nk1+nk2+nk3)))*(Y(iii+1)**2+Y(iii+3)**2)+(-Y(iii+1)*Y(iii+2)+Y(iii+3)*Y(iii))
c - +(k(nk1+nk2+nk3)/2.d0)*(Y(iii)**2+Y(iii+2)**2)
c Now print out the real and imaginary parts of the modes
n = 0
DO i=5,NTOT-3,4
If(n.ne.nk1) then
n = n+1
else
n = n+2
endif
c Note that this way I have all the modes at each value of t
c Write(20,*) t,k(n),Y(i),Y(i+1),Y(i+2),Y(i+3)
c I need to print it out twice if n = nk1.
c If(n.eq.nk1) then
c Write(20,*) t,k(n),Y(i),Y(i+1),Y(i+2),Y(i+3)
c Endif
Enddo
If(iback.eq.0.and.dabs(psi2).gt.1.d20) then
Print *,'phi2 has reached 1.d20'
Stop
Endif
27 format(10(e21.14,1x))
end do
STOP
END
Subroutine Intconst(delk1,delk2,delk3)
!$$$c This subroutine computes the constants needed for the mode integrations.
!$$$It uses Simpson's rule.
!$$$c The integrand is assumed to be equal to zero for k = 0.
!$$$c If the number of intervals is wrong use the trapezoidal rule.
Implicit None
Integer n, nm, nk1,nk2,nk3, j
Real*8 c1,c2,c3,c12,c22,c32,k,kk,delk1,delk2,delk3,lamk
Parameter (nm = 24000)
Dimension c1(nm),c2(nm),c3(nm),k(3*nm),c12(nm),c22(nm),c32(nm)
Common/in1/lamk,c1,c2,c3,k
Common/in2/nk1,nk2,nk3
Common/in3/c12,c22,c32
!$$$This subroutine computes the intervals of integration.
!$$$There are two intervals because we have a lower limit cutoff and so want
!$$$to integrate a certain set up terms up to that and then all other terms
!$$$starting at the cutoff.
!$$$Note that for the last interval we want to use the trapezoidal rule
!$$$if the number of intervals does not work out right. This changes
!$$$the value of c(n)for the last two values of n. The same thing happens
!$$$at the end of the second interval of integration.
!$$$
!$$$
c
c I have now split the second interval into two parts to make the integration more accurate.
c
kk = 0.d0
j = 0
do n=1,nk1
kk = kk+delk1
k(n) = kk
if(n.eq.nk1) then
If(j.eq.0) then
c1(n) = delk1/2.d0
else
c1(n) = delk1/3.d0
endif
else if(n.eq.nk1-1) then
if(j.eq.1) then
c1(n) = delk1*(1.d0/3.d0 + 0.5d0)
j = 0
else
c1(n) = 4.d0*delk1/3.d0
j = 1
endif
else if(j.eq.0) then
c1(n) = 4.d0*delk1/3.d0
j = 1
else
c1(n) = 2.d0*delk1/3.d0
j = 0
endif
enddo
!$$$
!$$$Here we start at the lower limit cutoff so the sequence is different.
!$$$Also k(nk1+1) = k(nk1) because of the way the limits of integration are.
j = 0
kk = kk - delk2
do n=1,nk2
kk = kk + delk2
k(n+nk1) = kk
if(n.eq.nk2) then
If(j.eq.0) then
c2(n) = delk2/2.d0
else
c2(n) = delk2/3.d0
endif
else if(n.eq.nk2-1) then
if(j.eq.1) then
c2(n) = delk2*(1.d0/3.d0 + 0.5d0)
j = 0
else
c2(n) = 4.d0*delk2/3.d0
j = 1
endif
else if (n.eq.1) then
c2(n) = delk2/3.d0
else if(j.eq.0) then
c2(n) = 4.d0*delk2/3.d0
j = 1
else
c2(n) = 2.d0*delk2/3.d0
j = 0
endif
enddo
!$$$
!$$$Here we start at the lower limit cutoff so the sequence is different.
!$$$Also k(nk1+nk2+1) = k(nk1+nk2) because of the way the limits of integration are.
j = 0
kk = kk - delk3
do n=1,nk3
kk = kk + delk3
k(n+nk1+nk2) = kk
if(n.eq.nk3) then
If(j.eq.0) then
c3(n) = delk3/2.d0
else
c3(n) = delk3/3.d0
endif
else if(n.eq.nk3-1) then
if(j.eq.1) then
c3(n) = delk3*(1.d0/3.d0 + 0.5d0)
j = 0
else
c3(n) = 4.d0*delk3/3.d0
j = 1
endif
else if (n.eq.1) then
c3(n) = delk3/3.d0
else if(j.eq.0) then
c3(n) = 4.d0*delk3/3.d0
j = 1
else
c3(n) = 2.d0*delk3/3.d0
j = 0
endif
enddo
c Now do the intervals when I have 1/2 the number of points
c I am going to simply assume there is an odd number of points
c and that the intervals work out correctly
c You need (nk-1) to be divisible by 4 for the second two bands.
c For the first band you need nk1 to be divisible by 4 since you
c start with n = 0 (i.e. k = 0) but that doesn't contribute.
c Throughout delk is multiplied by 2.
j = 0
c We start with n = 2 because of the fact that we are skipping every other step and n = 0 gives
c no contribution because k = 0
do n=2,nk1,2
if(n.eq.nk1) then
c12(n) = 2.d0*delk1/3.d0
else if(j.eq.0) then
c12(n) = 4.d0*2.d0*delk1/3.d0
j = 1
else
c12(n) = 2.d0*2.d0*delk1/3.d0
j = 0
endif
enddo
!$$$
!$$$Here we start at the lower limit cutoff so the sequence is different.
!$$$Also k(nk1+1) = k(nk1) because of the way the limits of integration are.
c Here we start with n = 1 because the beginning of the interval does contribute
j = 0
do n=1,nk2,2
if(n.eq.nk2.or.n.eq.1) then
c22(n) = 2.d0*delk2/3.d0
else if(j.eq.0) then
c22(n) = 4.d0*2.d0*delk2/3.d0
j = 1
else
c22(n) = 2.d0*2.d0*delk2/3.d0
j = 0
endif
enddo
!$$$
!$$$Here we start at the lower limit cutoff so the sequence is different.
!$$$Also k(nk1+nk2+1) = k(nk1+nk2) because of the way the limits of integration are.
j = 0
do n=1,nk3,2
if(n.eq.nk3.or.n.eq.1) then
c32(n) = 2.d0*delk3/3.d0
else if(j.eq.0) then
c32(n) = 4.d0*2.d0*delk3/3.d0
j = 1
else
c32(n) = 2.d0*2.d0*delk3/3.d0
j = 0
endif
enddo
Return
End
!$$$This is the subroutine called by the differential equation solver to
!$$$step the equations forward in time.
SUBROUTINE DERIVS(T,Y,DYDT,NVAR)
IMPLICIT None
Integer n,nk1,nk2,nk3,nvar,nm,i,ntest,nt,iback
Real*8 t,Y, DYDT, k, om2, c1, c2, c3,psi2, lamphi, gnew, f2, mphi
Real*8 mu,capom2,a,ap,app,phi,phip,xi,pi,psi2a,lamk
Real*8 f2p,int1,int2,int3,rhopsi,rhopsia,rhophi,psi21,psi22,psi23
Real*8 int11,psi2int
Parameter (nm = 24000)
DIMENSION Y(NVAR),DYDT(NVAR),c1(nm),c2(nm),c3(nm),k(3*nm)
common/param/pi,mu,mphi,xi,lamphi,gnew
Common/in1/lamk,c1,c2,c3,k
Common/in2/nk1,nk2,nk3
Common/der/psi2,psi2a,rhophi,rhopsi,rhopsia,psi21,psi22,psi23,psi2int
Common/back/iback
!$$$ Y(1) is a, Y(2) is ap, Y(3) is phi, y(4) is phip. DYDT(1) is ap,
!$$$ DYDT(2) is ap, DYDT(3) is phip, DYDT(4) is phipp.
!$$$ 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.
!$$$ This version of the program is for flat space. That's why
!$$$I explicitly set a = 1, ap = app = 0.
a = 1.d0
ap = 0.d0
app = 0.d0
phi = Y(3)
phip = Y(4)
DYDT(1) = ap
DYDT(3) = phip
rhophi = (phip**2 + mphi**2*phi**2)/2.d0
!$$$ We have the same value of k at n = nk1 and n = nk1+1.
!$$$c This is because we want different values for c1(nk1)
!$$$c and c2(nk1+1) even though the value of k is the same
!$$$c for these since they are at the end of the first integration
!$$$c interval and the beginning of the second respectively.
!$$$We first update the values of the derivatives of the modes.
n = 0
DO i=5,NVAR-3,4
If(n.ne.nk1) then
n = n+1
else
n = n+2
endif
capom2 = k(n)**2+gnew**2*a**2*phi**2
DYDT(i) = Y(i+1)
DYDT(i+1) = -capom2*Y(i)
DYDT(i+2) = Y(i+3)
DYDT(i+3) = -capom2*Y(i+2)
end do
!$$$psi2 is <\psi^2>. Recall that c1(n) = const. delk1 and
!$$$c2(n) = const. delk2.
psi2 = 0.D0
rhopsi = 0.d0
DO n = 1,nk1
i = 4*n + 1
f2 = (Y(i)**2+Y(i+2)**2)
f2p = (Y(i+1)**2+Y(i+3)**2)
INT1 = k(n)**2*f2 - k(n)/2.d0
int2 = k(n)**2*(f2p+(k(n)**2+a**2*gnew**2*phi**2)*f2)
- - k(n)**3 - k(n)*gnew**2*phi**2*a**2/2.
psi2 = psi2 + int1*c1(n)
rhopsi = rhopsi + int2*c1(n)
Enddo
psi21 = psi2
!$$$Now n is for c2(n) and recall k(nk1) = k(nk1+1)
psi22 = 0.d0
DO n = 1,nk2
i = 4*(n-1+nk1) + 1
f2 = (Y(i)**2+Y(i+2)**2)
f2p = (Y(i+1)**2+Y(i+3)**2)
int1 = k(n+nk1)**2*f2 - k(n+nk1)/2.d0 +
- gnew**2*phi**2*a**2/(4.d0*k(n+nk1))
int2 = k(n+nk1)**2*(f2p+(k(n+nk1)**2+a**2*gnew**2*phi**2)*f2)
- - k(n+nk1)**3 - k(n+nk1)*gnew**2*phi**2*a**2/2.d0 +
- gnew**4*phi**4*a**4/(8.d0*k(n+nk1))
psi2 = psi2 + int1*c2(n)
psi22 = psi22 + int1*c2(n)
rhopsi = rhopsi + int2*c2(n)
Enddo
psi23 = 0.d0
DO n = 1,nk3
i = 4*(n-1+nk1+nk2) + 1
f2 = (Y(i)**2+Y(i+2)**2)
f2p = (Y(i+1)**2+Y(i+3)**2)
c int11 = k(n+nk1+nk2)**2*f2
int1 = k(n+nk1+nk2)**2*f2 - k(n+nk1+nk2)/2.d0 +
- gnew**2*phi**2*a**2/(4.d0*k(n+nk1+nk2))
int2 = k(n+nk1+nk2)**2*(f2p+(k(n+nk1+nk2)**2+a**2*gnew**2*phi**2)*f2)
- - k(n+nk1+nk2)**3 - k(n+nk1+nk2)*gnew**2*phi**2*a**2/2.d0 +
- gnew**4*phi**4*a**4/(8.d0*k(n+nk1+nk2))
psi2 = psi2 + int1*c3(n)
psi23 = psi23 + int1*c3(n)
rhopsi = rhopsi + int3*c3(n)
Enddo
!$$$c Now put in the normalization factor.
psi2 = psi2/(2*pi**2*a**2)
rhopsi = rhopsi/(4*pi**2*a**4)
!$$$c Now add on the analytic approximation.
!$$$Note that lamk is the lower limit cutoff. We called it epsilon in our paper.
!$$$c Note also that when we do backreaction of the geometry the app will
!$$$c be a problem that will have to be dealt with. It will have to be
!$$$c computed from the backreaction equations before psi2a is computed.
psi2a = -app/(48*pi**2*a**3)-gnew**2*phi**2/(8*pi**2)*(1 - dlog(2*lamk/(a*mu)))
psi2 = psi2 + psi2a
c rhopsia = - gnew**4*phi**4/(32.d0*pi**2) + gnew**2*phi**4/(32.d0*pi**2) *dlog(2*lamk/(a*mu))
rhopsia = - gnew**4*phi**4/(32.d0*pi**2) + gnew**4*phi**4/(32.d0*pi**2) *dlog(2*lamk/(a*mu))
rhopsi = rhopsi + rhopsia
DYDT(4) = -2.d0*ap/a*phip - a**2*(mphi**2 + xi * 6.d0*DYDT(2)/a**3
- + lamphi/6.d0*phi**2 + iback*gnew**2*psi2)*phi
c Finally put in the terms from the 4th order WKB approximation for modes above the cutoff
If(iback.eq.1) then
DYDT(4) = (1.d0/(1.d0+gnew**4*phi**2/(32.d0*pi**2*k(nk1+nk2+nk3)**2)))*(DYDT(4)
1 - 3.d0*gnew**6*phi**5/(64.d0*pi**2*k(nk1+nk2+nk3)**2)
2 - gnew**4*phip**2*phi/(32.d0*pi**2*k(nk1+nk2+nk3)**2))
psi2 = psi2 + (1.d0/(4.d0*pi**2*k(nk1+nk2+nk3)**2))*(3.d0*gnew**4*phi**4/16.d0
1 + gnew**2*phip**2/8.d0 + gnew**2*phi*DYDT(4)/8.d0)
Endif
c
c This is enough if we work in a fixed background. Then I just need
c specify DYDT(2) which is a''.
DYDT(2) = 0.d0
RETURN
END
!$$$This subroutine computes the modes sums when half the modes are used.
SUBROUTINE TEST(t,Y,nvar)
IMPLICIT None
Integer nm,nk1,nk2,nk3,j,n,i,nvar
Real*8 t,c1,c2,c3,k,w,pi,mphi,xi,lamphi,gnew,int,dif1,dif2,dif3,a
Real*8 lamk,phi,psi2a,psi2,mu,rhophi,rhopsi,rhopsia
Real*8 c12,c22,c32,psi21,psi22,psi23,psi2t1,psi2t2,psi2t3,psi2t
Real*8 psi2int,diftot,f2,Y,int1
Parameter(nm = 24000)
DIMENSION Y(NVAR)
Dimension c1(nm),c2(nm),c3(nm),k(3*nm)
Dimension c12(nm),c22(nm),c32(nm)
Common/der/psi2,psi2a,rhophi,rhopsi,rhopsia,psi21,psi22,psi23,psi2int
Common/in1/lamk,c1,c2,c3,k
Common/in2/nk1,nk2,nk3
Common/in3/c12,c22,c32
common/param/pi,mu,mphi,xi,lamphi,gnew
a = 1.d0
phi = Y(3)
J = 0
psi2t1 = 0.d0
DO n = 2,nk1,2
i = 4*n + 1
f2 = (Y(i)**2+Y(i+2)**2)
INT1 = k(n)**2*f2 - k(n)/2.d0
psi2t1 = psi2t1 + int1*c12(n)
Enddo
psi2t2 = 0.d0
c We are starting with n = 0 really although that gives zero contribution.
c Since we skip every other mode we should start thus with n = 2.
c For the other bands we start with n = 1 because the beginning of the band does contribute
DO n = 1,nk2,2
i = 4*(n-1+nk1) + 1
f2 = (Y(i)**2+Y(i+2)**2)
int1 = k(n+nk1)**2*f2 - k(n+nk1)/2.d0 +
- gnew**2*phi**2*a**2/(4.d0*k(n+nk1))
psi2t2 = psi2t2 + int1*c22(n)
Enddo
psi2t3 = 0.d0
DO n = 1,nk3
i = 4*(n-1+nk1+nk2) + 1
f2 = (Y(i)**2+Y(i+2)**2)
int1 = k(n+nk1+nk2)**2*f2 - k(n+nk1+nk2)/2.d0 +
- gnew**2*phi**2*a**2/(4.d0*k(n+nk1+nk2))
psi2t3 = psi2t3 + int1*c32(n)
Enddo
psi2t = (psi2t1+psi2t2+psi2t3)/(2*pi**2*a**2) + psi2a
dif1 = (psi21-psi2t1)/psi21
dif2 = (psi22-psi2t2)/psi22
dif3 = (psi23-psi2t3)/psi23
diftot = (psi2-psi2t)/psi2
Write(16,17) t,psi2t1,psi2t2,psi2t3
WRITE(16,17) t,dif1,dif2,dif3,diftot
17 format(8(e12.5,2x))
RETURN
END
c
c This subroutine is used to check the convergence of the mode sum.
c
c In version 2 there was a program that was checking to see what
c happens to the integrand when you should have the modes canceling
c with the counterterms. The upshot seemed to be that the cancelation
c gets worse at large k because the integrand goes like k^2 or something
c like that.
c
SUBROUTINE TestModeSum(T,Y,DYDT,NVAR)
IMPLICIT None
Integer n,nk1,nk2,nk3,nvar,nm,i,ntest,nt,ifile,nfile,iback
Integer nn,j,nnk2,nnk3
Real*8 t,Y, DYDT, k, om2, c1, c2, c3, psi2, lamphi, gnew, f2, f2p,int1,int2, mphi
Real*8 mu,capom2,a,ap,app,phi,phip,phipp,xi,pi,psi2a,lamk
Real*8 difk,difk1,psi2dif,psi2old,inttest,inttest1,cc,rhophi,rhopsi,rhopsia
Real*8 beta2,psi21,psi22,psi23,psi2int
Parameter (nm = 24000)
DIMENSION Y(NVAR),DYDT(NVAR),c1(nm),c2(nm),c3(nm),k(3*nm),cc(3)
common/param/pi,mu,mphi,xi,lamphi,gnew
Common/in1/lamk,c1,c2,c3,k
Common/in2/nk1,nk2,nk3
c Common/der/psi2,psi2a,rhophi,rhopsi,rhopsia
Common/der/psi2,psi2a,rhophi,rhopsi,rhopsia,psi21,psi22,psi23,psi2int
Common/fil/ifile
Common/back/iback
ifile = ifile+1
nfile = ifile+30
If(ifile.eq.1) then
Open(nfile,file='integrand.1')
else if(ifile.eq.2) then
Open(nfile,file='integrand.2')
else if(ifile.eq.3) then
Open(nfile,file='integrand.3')
else if(ifile.eq.4) then
Open(nfile,file='integrand.4')
else if(ifile.eq.5) then
Open(nfile,file='integrand.5')
else if(ifile.eq.6) then
Open(nfile,file='integrand.6')
else if(ifile.eq.7) then
Open(nfile,file='integrand.7')
else if(ifile.eq.8) then
Open(nfile,file='integrand.8')
else if(ifile.eq.9) then
Open(nfile,file='integrand.9')
else if(ifile.eq.10) then
Open(nfile,file='integrand.10')
else if(ifile.eq.11) then
Open(nfile,file='integrand.11')
else if(ifile.eq.12) then
Open(nfile,file='integrand.12')
else if(ifile.eq.13) then
Open(nfile,file='integrand.13')
else if(ifile.eq.14) then
Open(nfile,file='integrand.14')
else if(ifile.eq.15) then
Open(nfile,file='integrand.15')
else if(ifile.eq.16) then
Open(nfile,file='integrand.16')
else if(ifile.eq.17) then
Open(nfile,file='integrand.17')
else if(ifile.eq.18) then
Open(nfile,file='integrand.18')
else if(ifile.eq.19) then
Open(nfile,file='integrand.19')
else if(ifile.eq.20) then
Open(nfile,file='integrand.20')
else if(ifile.eq.21) then
Open(nfile,file='integrand.21')
Endif
!$$$ Y(1) is a, Y(2) is ap, Y(3) is phi, y(4) is phip. DYDT(1) is ap,
!$$$ DYDT(2) is ap, DYDT(3) is phip, DYDT(4) is phipp.
!$$$ 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.
phi = Y(3)
phip = Y(4)
phipp = DYDT(4)
a = 1.d0
!$$$ We have the same value of k at n = nk1 and n = nk1+1.
!$$$c This is because we want different values for c1(nk1)
!$$$c and c2(nk1+1) even though the value of k is the same
!$$$c for these since they are at the end of the first integration
!$$$c interval and the beginning of the second respectively.
c Write(18,*) 't = ',t
c Write(18,*) ' n k(n) Real f Imag f int1 int2 psi2 psi2dif'
!$$$ psi2 is <\psi^2>. Recall that c1(n) = const. delk1 and
!$$$ c2(n) = const. delk2.
cc(1) = c2(1)
cc(2) = c2(2)
cc(3) = cc(1)
psi2 = 0.D0
psi2old = psi2
DO n = 1,nk1
i = 4*n + 1
f2 = (Y(i)**2+Y(i+2)**2)
f2p = (Y(i+1)**2+Y(i+3)**2)
INT1 = k(n)**2*f2 - k(n)/2.d0
c This was used to check that the sums
c This is the % difference.
c inttest = f2 - 1.d0/(2.d0*k(n))
c - + gnew**2*phi**2*a**2/(4.d0*k(n)**3)
c difk = inttest/f2
c difk1 = -1.d0
c This is the integrand for the energy density of the psi field
int2 = k(n)**2*(f2p+(k(n)**2+a**2*gnew**2*phi**2)*f2)
- - k(n)**3 - k(n)*gnew**2*phi**2*a**2/2.
psi2 = psi2 + int1*c1(n)
c
c This is to check the convergence of the modes.
c
psi2dif = (psi2-psi2old)/psi2
psi2old = psi2
c beta2 = |beta|^2
beta2 = (1.d0/(2.D0*k(n)))*(Y(i+1)**2+Y(i+3)**2)+(-Y(1)*Y(2)+Y(3)*Y(1))
- +(k(n)/2.d0)*(Y(i)**2+Y(i+2)**2)
c Write(nfile,37) t,n,k(n),Y(i),Y(i+2),int1,int2,psi2,psi2dif
Write(nfile,37) t,n,k(n),Y(i),beta2,int1,int2,psi2,psi2dif
Enddo
c psi21 = psi2
37 Format(1x,d12.5,1x,i5,1x,7(d12.5,1x))
!$$$Now n is for c2(n) and recall k(nk1) = k(nk1+1)
c This is to make the intervals work better.
c
c We need an odd number of n's because there are 3 n's in each interval.
c We integrate in each interval from nn-1 to nn+1.
nnk2 = ((nk2-1)/2)*2
c psi22 = 0.d0
DO nn = 2,nnk2,2
Do j=1,3
n = nn-2 + j
i = 4*(n-1+nk1) + 1
f2 = (Y(i)**2+Y(i+2)**2)
f2p = (Y(i+1)**2+Y(i+3)**2)
int1 = k(n+nk1)**2*f2 - k(n+nk1)/2.d0 +
- gnew**2*phi**2*a**2/(4.d0*k(n+nk1))
c inttest = f2 - 1.d0/(2.d0*k(n+nk1)) +
c - gnew**2*phi**2*a**2/(4.d0*k(n+nk1)**3)
c difk = inttest/f2
c This is to test whether a different subtraction would lead to
c a faster convergence of the mode sums.
c inttest1 = f2 - 1.d0/(2.d0*dsqrt(k(n+nk1)**2 +
c 1 gnew**2*phi**2*a**2))
c difk1 = inttest1/f2
int2 = k(n+nk1)**2*(f2p+(k(n+nk1)**2+a**2*gnew**2*phi**2)*f2)
- - k(n+nk1)**3 - k(n+nk1)*gnew**2*phi**2*a**2/2. +
- gnew**4*phi**4*a**4/(8.*k(n+nk1))
psi2 = psi2 + int1*cc(j)
c psi22 = psi22 + int1*cc(j)
End Do
psi2dif = (psi2-psi2old)/psi2
psi2old = psi2
beta2 = (1.d0/(2.D0*k(n+nk1)))*(Y(i+1)**2+Y(i+3)**2)+(-Y(1)*Y(2)+Y(3)*Y(1))
- +(k(n+nk1)/2.d0)*(Y(i)**2+Y(i+2)**2)
c Write(nfile,37) t,n,k(n+nk1),Y(i),Y(i+2),int1,int2,psi2,psi2dif
Write(nfile,37) t,n,k(n+nk1),Y(i),beta2,int1,int2,psi2,psi2dif
Enddo
!$$$Now n is for c3(n) and recall k(nk1+nk2) = k(nk1+nk2+1)
!$$$ psi2 is <\psi^2>. Recall that c1(n) = const. delk1 and
!$$$ c2(n) = const. delk2, c3(n) = const*delk3
cc(1) = c3(1)
cc(2) = c3(2)
cc(3) = c3(1)
c This is to make the intervals work better.
c
c We need an odd number of n's because there are 3 n's in each interval.
c We integrate in each interval from nn-1 to nn+1.
nnk3 = ((nk3-1)/2)*2
c psi23 = 0.d0
DO nn = 2,nnk3,2
Do j=1,3
n = nn-2 + j
i = 4*(n-1+nk1+nk2) + 1
f2 = (Y(i)**2+Y(i+2)**2)
f2p = (Y(i+1)**2+Y(i+3)**2)
int1 = k(n+nk1+nk2)**2*f2 - k(n+nk1+nk2)/2.d0 +
- gnew**2*phi**2*a**2/(4.d0*k(n+nk1+nk2))
c inttest = f2 - 1.d0/(2.d0*k(n+nk1+nk2)) +
c - gnew**2*phi**2*a**2/(4.d0*k(n+nk1+nk2)**3)
c difk = inttest/f2
c This is to test whether a different subtraction would lead to
c a faster convergence of the mode sums.
c inttest1 = f2 - 1.d0/(2.d0*dsqrt(k(n+nk1+nk2)**2 +
c 1 gnew**2*phi**2*a**2))
c difk1 = inttest1/f2
int2 = k(n+nk1+nk2)**2*(f2p+(k(n+nk1+nk2)**2+a**2*gnew**2*phi**2)*f2)
- - k(n+nk1+nk2)**3 - k(n+nk1+nk2)*gnew**2*phi**2*a**2/2. +
- gnew**4*phi**4*a**4/(8.*k(n+nk1+nk2))
psi2 = psi2 + int1*cc(j)
c psi23 = psi23 + int1*cc(j)
End Do
psi2dif = (psi2-psi2old)/psi2
psi2old = psi2
beta2 = (1.d0/(2.D0*k(n+nk1+nk2)))*(Y(i+1)**2+Y(i+3)**2)+(-Y(i+1)*Y(i+2)+Y(i+3)*Y(i))
- +(k(n+nk1+nk2)/2.d0)*(Y(i)**2+Y(i+2)**2)
c Write(nfile,37) t,n,k(n+nk1+nk2),Y(i),Y(i+2),int1,int2,psi2,psi2dif
Write(nfile,37) t,n,k(n+nk1+nk2),Y(i),beta2,int1,int2,psi2,psi2dif
Enddo
!$$$c Now put in the normalization factor.
psi2 = psi2/(2*pi**2*a**2)
!$$$c Now add on the analytic approximation.
!$$$Note that lamk is the lower limit cutoff. We called it epsilon in our paper.
!$$$c Note also that when we do backreaction of the geometry the app will
!$$$c be a problem that will have to be dealt with. It will have to be
!$$$c computed from the backreaction equations before psi2a is computed.
psi2a = -app/(48*pi**2*a**3)-gnew**2*phi**2/(8*pi**2)*(1 - dlog(2*lamk/(a*mu)))
psi2 = psi2 + psi2a
c Now add on the 4'th order wkb approximation for modes above the cutoff
If(iback.eq.1) then
psi2 = psi2 + (1.d0/(4.d0*pi**2*k(nk1+nk2+nk3)**2))*(3.d0*gnew**4*phi**4/16.d0
1 + gnew**2*phip**2/8.d0 + gnew**2*phi*phipp/8.d0)
Endif
RETURN
END