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,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,cfit
	Real*8 G,H3,f,af,ap1,ap2,ap3,ap4
	Real*8 u1,zeta3,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),H01(NN),H03(NN),Traan(NN)
        Dimension G0(NN),rhos2(NN),rhos3(NN),Td(NN,nt),a4(NN),conan(NN)
        Dimension Hu1(NN),Hu3(NN),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='modes2m2i1t4.dat')
        open(16,FILE='SETresult2m2i1t4.dat')
        open(14,FILE = 'scalefresult2m2i1t4.dat')
c        open(13,FILE='h1h3info0m2i1t3.dat')
        open(13,FILE='rhotest2m2i1t4.dat')
        open(12,FILE='Ttest2m2i1t4.dat')
        open(17,FILE='ap2fit2m2i1t4.dat')
        open(18,FILE='ap1fit2m2i1t4.dat')
        open(19,FILE='afit2m2i1t4.dat')
c        open(13,FILE='SETtest.dat')
c        Open(14,File='initialD4.dat')
c        Open(13,File='initialDY4.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(17,*)  cfit(j,3)
        enddo

	pi = 2.d0*dasin(1.d0)
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
c        Y(1) =dsqrt(3.d0/lamda)*a(dsqrt(lamda/3.d0)*uin)
        Y(1)=a(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)
c        a1(it) = -dsqrt(-Y(1)**2+lamda*Y(1)**4/3.d0)
c        a2(it) = -Y(1)+2.d0*lamda*Y(1)**3/3.d0
c        a3(it) = -a1(it)+2.d0*lamda*Y(1)**2*a1(it)
c        a4(it) = Y(1)-20.d0*lamda*Y(1)**3/3.d0
c     +  +8.d0*lamda**2*Y(1)**5/3.d0
c        Y(2) = a1(it)
        H01(it) =-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


        H03(it) = 3.d0*a1(it)**4/a(uin)**8+6.d0*a1(it)**2/a(uin)**6
     +  +3.d0/a(uin)**4

        Hu1(it)=-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

        Hu3(it)=12.d0*a2(it)*a1(it)**2/a(uin)**7-12.d0*a1(it)**4/a(uin)**8
     +  +12.d0*a2(it)/a(uin)**5-12.d0*a1(it)**2/a(uin)**6
        Y(2) = a1(it)
        addr = a2(it)
c        write(13,21) H01(it), H03(it),Hu1(it),Hu3(it)
c         write(14,21) a1(it),a2(it),a3(it),a4(it) 
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
c this is the initial condition for second order wkb       
c           Wd = -(k(n)**2-0.25d0)*adot(uin)/(W*a(uin)**3)
         k(n) = n
         Omega = dsqrt(k(n)**2 + m**2*Y(1)**2)
c this is the initial condition for second order wkb       
c         W = Omega+5.d0*m**4*Y(1)**2*Y(2)**2/(8.d0*Omega**5)
c     -   -m**2*Y(2)**2/(4.d0*Omega**3)-m**2*Y(1)*addr/(4.d0*Omega**3)
c        W = Omega        
         W = Omega+5.d0*m**4*Y(1)**2*Y(2)**2/(8.d0*Omega**5)
     -   -m**2*Y(2)**2/(4.d0*Omega**3)-m**2*Y(1)*addr/(4.d0*Omega**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))
c           f(i+1) = -m**2*Y(1)*Y(2)/(Omega*(2.d0*W)**(1.5d0))
c           -Wd/(2.d0*W)**1.5d0
           f(i+1) = -m**2*Y(1)*Y(2)/(Omega*(2.d0*W)**(1.5d0))
           f(i+2) = 0.d0
           f(i+3) = -dsqrt(W/2.D0)

c           -Wd/(2.d0*W)**1.5d0
        
   
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 
        n = 0
        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)<reldiff1) exit
        rhos(it) = rhos(it)+rho(it,n)-rhod(it,n)
c        write(16,23) n,rho(it,n), rhod(it,n),rho(it,n)-rhod(it,n)   
        Enddo
        G0(it) =-3.d0*a1(it)**2/a(uin)**4-3.d0/a(uin)**2



c        rhoan(it) =m**2*G0(it)/(288.d0*pi**2)-m**4*(1.d0/2.d0
         rhoan(it)=(-H01(it)/6.d0+H03(it))/(2880.d0*pi**2)
     +  +m**2*G0(it)/(288.d0*pi**2)-m**4*(1.d0/2.d0
     +  +DlOG(m**2*a(uin)**2/4.d0)+2.d0*Eulergamma)/(64.d0*pi**2)
c        write(16,*) rhos(it),rhoan(it),rhos(it)+rhoan(it)
c        Y(2) = -dsqrt(8.d0*pi*(rhos(it)+rhoan(it))*Y(1)**4/3.d0
c     -  -Y(1)**2+lamda*Y(1)**4/3.d0)
c        Y(2) = a1(it)/Y(1)
        Y(2) = -sqrt(8.d0*pi*(rhos(1)+rhoan(1))*Y(1)**2/3.d0-1.d0
     +  +lamda*Y(1)**2/3.d0)

c        ad = Y(2)/Y(1)
c        Y(2) = a1(it)/Y(1)
c        ad2 = (a2(it)-Y(1)*ad**2)/(Y(1)**2)
c        Y(2) = ad	
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(5(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        
        it = 1
         
        Do i = 3,ntot-3,4
        Y(i) = Sqrt(Y(1))*f(i)
        Y(i+1) = 0.5d0*ad1(uin)*f(i)/Sqrt(Y(1))+f(i+1)/Sqrt(Y(1))
        Y(i+2) = 0.d0
        Y(i+3) = f(i+3)/Sqrt(Y(1))

        enddo

        n = 0
        Do i =3,ntot-3,4
        n = n+1 
        T(it,n)=k(n)**2*(Y(i)**2+Y(i+2)**2)*m**2/(2.d0*pi**2*Y(1)**3)

        Td(it,n)=(k(n)*m**2*Y(1)**2-m**4*Y(1)**4/(2.d0*k(n)))
     /  /(4*pi**2*Y(1)**4)
        reldifft(it,n) = Abs((T(it,n)-Td(it,n))/Td(it,n))
        if(reldifft(it,n)<reldiff1) exit
c        write(16,23) n,T(it,n),Td(it,n),T(it,n) - Td(it,n)
        Tre(it) = Tre(it) + T(it,n) - Td(it,n)
        enddo
        Gu(it) = -6.d0*a2(it)/a(uin)**3-6.d0/a(uin)**2
        
        Traan(it) = (-Hu1(it)/6.d0+Hu3(it))/(2880.d0*pi**2) 
     +  +m**2*Gu(it)/(288.d0*pi**2)-m**4*(1.d0
c        Traan(it) = m**2*Gu(it)/(288.d0*pi**2)-m**4*(1.d0
     +  +DlOG(m**2*a(uin)**2/4.d0)+2.d0*Eulergamma)/(16.d0*pi**2) 

        addr =-1.d0/Y(1)-Y(2)**2/Y(1)+4.d0*pi*(Tre(1)+Traan(1))*Y(1)/3.d0
     +  +2.d0*lamda*Y(1)/3.d0
c        Y(2) = -dsqrt(8.d0*pi*(rhos(it)+rhoan(it))*Y(1)**4/3.d0
c     -  -Y(1)**2+lamda*Y(1)**4/3.d0)
c        addr = -Y(1)-4.d0*pi*(Tre(it)+Traan(it))*Y(1)**3/3.d0
c     +  +2.d0*lamda*Y(1)**3/3.d0
c        ad = Y(2)/Y(1)
c        ad2 = (addr-Y(1)*ad**2)/(Y(1)**2)
c        Y(2) = ad      
c        ad2= -1.d0/Y(1)-Y(2)**2/Y(1)-4.d0*pi*(Tre(it)+Traan(it))*Y(1)/3.d0
c     +  +2.d0*lamda*Y(1)/3.d0
c        ap1 = a(uin)*ad1(uin)
c        ap2 = a(uin)*ad1(uin)**2+a(uin)**2*ad2(uin)
c        ap3 = a(uin)*ad1(uin)**3+4.d0*a(uin)**2*ad1(uin)*ad2(uin)
c     +  +a(uin)**3*ad3(uin)
c        ap4 = a(uin)*ad1(uin)**4+11.d0*a(uin)**2*ad1(uin)**2*ad2(uin)
c     +  +4.d0*a(uin)**3*ad2(uin)**2+7.d0*a(uin)**3*ad1(uin)*ad3(uin)
c     +  +a(uin)**4*ad4(uin)
c        write(14,21) ap1,ap2,ap3,ap4
        write(14,21) Y(1),Y(2),addr
c        write(13,21) H01(it), H03(it),Hu1(it),Hu3(it),G0(it)
c        write(16,21) a(uin),ad1(uin),ad2(uin),ad3(uin),ad4(uin)
c        write(13,21) a1(it),a2(it),a3(it),a4(it)
        write(16,23) it,rhos(it),rhoan(it),Tre(it),Traan(it)
	u = uin
        j=1
	DO it = 2,nu+1
	   ui = u
	   u = uin+delu*dble(it-1)
c           tk = (kk-1)*4+1
c           write(16,20) Y(tk),Y(tk+1),Y(tk+2),Y(tk+3)
c  call the ode solver 
c  compile it like 
c  link stdS stdS_1mod.o /home/anderson/source/obj/odebs.o
           CALL ODEINT(Y,ntot,ui,u,EPS,H1,HMIN,NOK,NBAD,DERIVS,BSSTEP)
        rhos(it) = 0
        Tre(it) = 0
        rhos2(it) = 0 
        n = 0
        DO i = 3,ntot-3,4
        n = n+1
c        wronski(i)=Y(i+1)*Y(i+2)-Y(i+3)*Y(i)-0.5d0
c  this is the method in the paper PhysRevD.61.024003
c  and PhysRevD.36.2963
         rho(it,n)=k(n)**2*(a(u)*(Y(i+1)**2+Y(i+3)**2)-ad1(u)*(Y(i)*Y(i+1)
     +  + Y(i+2)*Y(i+3))
     +  + (Y(i)**2+Y(i+2)**2)*(k(n)**2+m**2*a(u)**2
     +  + ad1(u)**2*0.25d0)/a(u))
     /  /(4.d0*pi**2*a(u)**4) 
        
        rhod(it,n) = (k(n)**3+k(n)*m**2*a(u)**2/2-m**4*a(u)**4/(8.d0*k(n)))
c     -  -(adot(u)**2+1)*(3.d0*m**2*a(u)**2)/(2.d0*k(n)))
     /  /(4*pi**2*a(u)**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)<reldiff1) exit
c        if(reldiff(it,n)<reldiff1) cycle
        rhos(it) = rhos(it)+rho(it,n)-rhod(it,n)
c        write(16,21) k(n), rho(it,n),rhod(it,n),rho(it,n)-rhod(it,n)        
c        write(16,*) wronski(i)
c        rhos(it) = rhos(it)+rho(it,n) -rhod(it,n)

c+rhoan(it,n)

     
c        reldiff(it,n) = Abs((T(it,n)-Td(it,n))/Td(it,n))
c        if(reldiff(it,n)<reldiff1) exit

c        write(16,*) k(n),T(it,n) - Td(it,n)
c        drho(it,i) = (adot(u)*(12.d0*ci-6.d0)*(Y(i+1)**2+Y(i+3)**2)
c     +  + 2.d0*adot(u)**2*(-42.d0*ci+9.d0)*(Y(i)*Y(i+1)
c     +  + Y(i+2)*Y(i+3))/a(u)
c     +  + adot(u)*(Y(i)**2+Y(i+2)**2)*((-12.d0*ci-2.d0)*k(n)**2
c     +  + 1.d0*(-72.d0*ci**2+2.d0)+m**2*a(u)**2*(-12.d0*ci)
c     +  +adot(u)**2*(-72.d0*ci**2+87.d0*ci-27.d0/2.d0)
c     -  - 12.d0*ci*(6.d0*ci-1.d0)*a(u)**2)/a(u)**2)
c     /  /(4.d0*pi**2*a(u)**4)

c        rhos3(it) = rhos3(it) + rhoan(it,n) 
c        write(16,*) k(n),T(it,n)-Td(it,n)
c        write(16,*) rho(it,i),T(it,i),drho(it,i)
c	write(16,*) rho(it,n), rhod(it,n),rho(it,n)-rhod(it,n)
	Enddo
        n=0
        DO i = 3,ntot-3,4
        n = n+1
c        wronski(i)=Y(i+1)*Y(i+2)-Y(i+3)*Y(i)-0.5d0
c  this is the method in the paper PhysRevD.61.024003
c  and PhysRevD.36.2963


c     -  -(adot(u)**2+1)*(3.d0*m**2*a(u)**2)/(2.d0*k(n)))
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  
c        if(reldiff(it,n)<reldiff1) cycle
c        write(16,20) k(n), rho(it,n),rhod(it,n),rho(it,n)-rhod(it,n)        
c        write(16,*) wronski(i)
c        rhos(it) = rhos(it)+rho(it,n) -rhod(it,n)

c+rhoan(it,n)

        T(it,n)=k(n)**2*(Y(i)**2+Y(i+2)**2)*m**2/(2.d0*pi**2*a(u)**3)

        Td(it,n)=(k(n)*m**2*a(u)**2-m**4*a(u)**4/(2.d0*k(n)))
     /  /(4*pi**2*a(u)**4)
        reldifft(it,n) = Abs((T(it,n)-Td(it,n))/Td(it,n))
        if(reldifft(it,n)<reldiff1) exit

        Tre(it) = Tre(it) + T(it,n) - Td(it,n)
c        write(16,23) n,T(it,n),Td(it,n),T(it,n) - Td(it,n)
c        drho(it,i) = (adot(u)*(12.d0*ci-6.d0)*(Y(i+1)**2+Y(i+3)**2)
c     +  + 2.d0*adot(u)**2*(-42.d0*ci+9.d0)*(Y(i)*Y(i+1)
c     +  + Y(i+2)*Y(i+3))/a(u)
c     +  + adot(u)*(Y(i)**2+Y(i+2)**2)*((-12.d0*ci-2.d0)*k(n)**2
c     +  + 1.d0*(-72.d0*ci**2+2.d0)+m**2*a(u)**2*(-12.d0*ci)
c     +  +adot(u)**2*(-72.d0*ci**2+87.d0*ci-27.d0/2.d0)
c     -  - 12.d0*ci*(6.d0*ci-1.d0)*a(u)**2)/a(u)**2)
c     /  /(4.d0*pi**2*a(u)**4)

c        rhos3(it) = rhos3(it) + rhoan(it,n) 
c        write(16,*) k(n),T(it,n)-Td(it,n)
c        write(16,*) rho(it,i),T(it,i),drho(it,i)
c       write(16,*) rho(it,n), rhod(it,n),rho(it,n)-rhod(it,n)
        Enddo
c        Y(1)=a(u)
c        a1(it) = -dsqrt(-Y(1)**2+lamda*Y(1)**4/3.d0)
c        a2(it) = -Y(1)+2.d0*lamda*Y(1)**3/3.d0
c        a3(it) = -a1(it)+2.d0*lamda*Y(1)**2*a1(it)
c        a4(it) = Y(1)-20.d0*lamda*Y(1)**3/3.d0
c     +  +8.d0*lamda**2*Y(1)**5/3.d0
        
c        a1(it) = -Y(1)**2+lamda*Y(1)**4/3.d0
c this is actually aprime^2
c        a2(it) = -Y(1)+2.d0*lamda*Y(1)**3/3.d0
c        a3(it) = -a1(it)+2.d0*lamda*Y(1)**2*a1(it)
c this is actually ap1*ap3
c        a4(it) = Y(1)-20.d0*lamda*Y(1)**3/3.d0
c     +  +8.d0*lamda**2*Y(1)**5/3.d0
        a1(it) = a(u)*ad1(u)
        a2(it) = a(u)*ad1(u)**2+a(u)**2*ad2(u)
        a3(it) = a(u)*ad1(u)**3+4.d0*a(u)**2*ad1(u)*ad2(u)
     +  +a(u)**3*ad3(u)
        a4(it) = a(u)*ad1(u)**4+11.d0*a(u)**2*ad1(u)**2*ad2(u)
     +  +4.d0*a(u)**3*ad2(u)**2+7.d0*a(u)**3*ad1(u)*ad3(u)
     +  +a(u)**4*ad4(u)
c        H01(it) =-36.d0*a3(it)/Y(1)**6+72.d0*a2(it)*a1(it)
c     /  /Y(1)**7+18.d0*a2(it)**2/Y(1)**6+36.d0*a1(it)/Y(1)**6
c     -  -18.d0/Y(1)**4


c        H03(it) = 3.d0*a1(it)**2/Y(1)**8+6.d0*a1(it)/Y(1)**6
c     +  +3.d0/Y(1)**4

c        Hu1(it)=-36.d0*a4(it)/Y(1)**5+144.d0*a3(it)/Y(1)**6
c     -  -216.d0*a2(it)*a1(it)/Y(1)**7+108.d0*a2(it)**2/Y(1)**6
c     +  +72.d0*a2(it)/Y(1)**5-72.d0*a1(it)/Y(1)**6

c        Hu3(it)=12.d0*a2(it)*a1(it)/Y(1)**7-12.d0*a1(it)**2/Y(1)**8
c     +  +12.d0*a2(it)/Y(1)**5-12.d0*a1(it)/Y(1)**6
        H01(it) =-36.d0*a3(it)*a1(it)/a(u)**6+72.d0*a2(it)*a1(it)**2
     /  /a(u)**7+18.d0*a2(it)**2/a(u)**6+36.d0*a1(it)**2/a(u)**6
     -  -18.d0/a(u)**4


        H03(it) = 3.d0*a1(it)**4/a(u)**8+6.d0*a1(it)**2/a(u)**6
     +  +3.d0/a(u)**4

        Hu1(it)=-36.d0*a4(it)/a(u)**5+144.d0*a3(it)*a1(it)/a(u)**6
     -  -216.d0*a2(it)*a1(it)**2/a(u)**7+108.d0*a2(it)**2/a(u)**6
     +  +72.d0*a2(it)/a(u)**5-72.d0*a1(it)**2/a(u)**6

        Hu3(it)=12.d0*a2(it)*a1(it)**2/a(u)**7-12.d0*a1(it)**4/a(u)**8
     +  +12.d0*a2(it)/a(u)**5-12.d0*a1(it)**2/a(u)**6


c        write(13,21) H01(it), H03(it),Hu1(it),Hu3(it) 
        Gu(it) = -6.d0*a2(it)/a(u)**3-6.d0/a(u)**2

c        write(13,21) H01(it),H03(it),Hu1(it),Hu3(it)
        Traan(it) =(-Hu1(it)/6.d0+Hu3(it))/(2880.d0*pi**2) 
     +  +m**2*Gu(it)/(288.d0*pi**2)-m**4*(1.d0
     +  +DlOG(m**2*a(u)**2/4.d0)+2.d0*Eulergamma)/(16.d0*pi**2)



        addr =-1.d0/Y(1)-Y(2)**2/Y(1)+4.d0*pi*(Tre(it)+Traan(it))*Y(1)/3.d0
     +  +2.d0*lamda*Y(1)/3.d0
        

        write(14,21) Y(1),Y(2),addr


c        H01(it) = -18.d0/a(u)**4-72.d0*adot(u)**2/a(u)**2
c     +  +54.d0*adot(u)**4/a(u)**4+36.d0*adot(u)**2/a(u)**4
c     +  +18
c        G0(it) = -3.d0*adot(u)**2/a(u)**2-3.d0/a(u)**2
        G0(it) = -3.d0*a1(it)**2/a(u)**4-3.d0/a(u)**2
c        write(13,21) H01(it), H03(it),Hu1(it),Hu3(it),G0(it)
c        H03(it) = 3.d0*adot(u)**4/a(u)**4+3.d0*adot(u)**2/a(u)**4
c     +  +3/a(u)**4
c        rhoan(it) = (-H01(it)/6.d0+H03(it))/(2880.d0*pi**2)
c     +  +m**2*G0(it)/(288.d0*pi**2)-m**4*(1.d0/2.d0
c     +  +DlOG(m**2*Y(1)**2/4.d0)+2.d0*Eulergamma)/(64.d0*pi**2)
c   this is calculating the analytical term of stress energy tensor       
c        G0(it) = -3.d0/Y(1)**2-3.d0*Y(2)**2/Y(1)**2



        rhoan(it) =(-H01(it)/6.d0+H03(it))/(2880.d0*pi**2) 
     +  +m**2*G0(it)/(288.d0*pi**2)-m**4*(1.d0/2.d0
     +  +DlOG(m**2*a(u)**2/4.d0)+2.d0*Eulergamma)/(64.d0*pi**2)

c        rhos(it) = rhos(it)+rhoan(it)
c        write(16,20) it,n,Tre(it),Traan(it),Tre(it)+Traan(it)
c write the final result to the file we want
c        write(16,20) it,n,rhos(it),rhoan(it),rhos(it)+rhoan(it)
c        Gu(it) = -6.d0*Y(2)**2/Y(1)**2-6.d0/Y(1)**2
c        Traan(it) =(-Hu1(it)/6.d0+Hu3(it))/(2880.d0*pi**2) 
c     +  +m**2*Gu(it)/(288.d0*pi**2)-m**4*(1.d0
c     +  +DlOG(m**2*Y(1)**2/4.d0)+2.d0*Eulergamma)/(16.d0*pi**2)

c        Traan(it) = m**2*Gu(it)/(288.d0*pi**2)-m**4*(1.d0
c     +  +DlOG(m**2*Y(1)**2/4.d0)+2.d0*Eulergamma)/(16.d0*pi**2)



c        addr=(-1.d0/Y(1)-Y(2)**2/Y(1)-4.d0*pi*(Tre(it)+Traan(it))*Y(1)/3.d0
c     +  +2.d0*lamda*Y(1)/3.d0)/(1.d0-m**2/(36.d0*pi))


c        write(14,21) Y(1),Y(2),addr

c        Gu(it) = -6.d0*Y(2)**2/Y(1)**2-6.d0/Y(1)**2
c     -  -6.d0*addr/Y(1)
c        Traan(it) = m**2*Gu(it)/(288.d0*pi**2)-m**4*(1.d0
c     +  +DlOG(m**2*Y(1)**2/4.d0)+2.d0*Eulergamma)/(16.d0*pi**2)
c        Traan(it) =(-Hu1(it)/6.d0+Hu3(it))/(2880.d0*pi**2) 
c     +  +m**2*Gu(it)/(288.d0*pi**2)-m**4*(1.d0
c     +  +DlOG(m**2*Y(1)**2/4.d0)+2.d0*Eulergamma)/(16.d0*pi**2)
c        G0(it) = -3.d0/Y(1)**2-3.d0*Y(2)**2/Y(1)**2

c        rhoan(it) =(-H01(it)/6.d0+H03(it))/(2880.d0*pi**2) 
c     +  +m**2*G0(it)/(288.d0*pi**2)-m**4*(1.d0/2.d0
c     +  +DlOG(m**2*Y(1)**2/4.d0)+2.d0*Eulergamma)/(64.d0*pi**2)
c        ap1 = a(u)*ad1(u)
c        ap2 = a(u)*ad1(u)**2+a(u)**2*ad2(u)
c        ap3 = a(u)*ad1(u)**3+4.d0*a(u)**2*ad1(u)*ad2(u)
c     +  +a(u)**3*ad3(u)
c        ap4 = a(u)*ad1(u)**4+11.d0*a(u)**2*ad1(u)**2*ad2(u)
c     +  +4.d0*a(u)**3*ad2(u)**2+7.d0*a(u)**3*ad1(u)*ad3(u)
c     +  +a(u)**4*ad4(u)
c        write(14,21) ap1,ap2,ap3,ap4
c        write(14,21) Y(1),Y(2),addr
c        write(16,21) a(u),ad1(u),ad2(u),ad3(u),ad4(u)
c        write(13,21) a1(it),a2(it),a3(it),a4(it)
c        rhoan(it) = m**2*G0(it)/(288.d0*pi**2)-m**4*(1.d0/2.d0
c     +  +DlOG(m**2*Y(1)**2/4.d0)+2.d0*Eulergamma)/(64.d0*pi**2)
        write(16,23) it,rhos(it),rhoan(it),Tre(it),Traan(it) 
c  This is so I can use them later 
        if((it-1)==(1+(j-1)*tk)) then
          j=j+1
          n = 0
          Do i =3 ,ntot-3,4
          n = n+1
          write(12,23) n,T(it-1,n),Td(it-1,n),T(it-1,n)-Td(it-1,n)

          write(13,23) n,rho(it-1,n),rhod(it-1,n),rho(it-1,n)-rhod(it-1,n)
c        write(16,*) n,rhos(tk),Tre(tk)
          Enddo
        endif

        Enddo
c        write(13,23) n,T(tk,n),Td(tk,n),T(tk,n)-Td(tk,n)
c        write(16,23) n,rho(tk,n), rhod(tk,n),rho(tk,n)-rhod(tk,n)   
c        write(16,*) n,rhos(tk),Tre(tk)
c        write(16,*) n,rhos(tk),Tre(tk),ad
c        Do i =1,nu
c        write(16,*) ck(i)
c        Enddo
c        write(16,*) rhos(tk),rhoan(tk),rhos(tk)+rhoan(tk)
c this part is used to check the conservation of stress energy tensor
c        u = uin
c        Do i = 1,nu-1
c        ui = u
c        u = u+delu
c        conan(i) = (rhoan(i+1)-rhoan(i))/delu+adot(u)*(4.d0*rhoan(i)
c     -  -Traan(i))/a(u)       
c        conan(i) = (rhoan(i+1)-rhoan(i))/delu
c        conan2(i) = adot(u)*(4.d0*rhoan(i) - Traan(i))/a(u)
        n = 0
        Do i = 3,ntot-3,4
        n = n+1
        cons(tk,n) =( rho(tk+1,n)-rho(tk,n))/delu
        cons2(tk,n) = Y(2)*(4.d0*rho(tk,n) - T(tk,n))/Y(1)
c        cons2(i,j) = drho(i,j)+adot(u)*(4.d0*rho(i,j)
c     -  - T(i,j))/a(u)  
c        drho2(i,j) = (rho(i+1,j)-rho(i,j))/delu
c        write(16,21) n,cons(tk,n),cons2(tk,n),cons(tk,n)+cons2(tk,n)
         
        enddo
        conan(tk) = (rhos(tk+1)-rhos(tk))/delu
        conan2(tk) = Y(2)*(4.d0*rhos(tk)-Tre(tk))/Y(1)
c        write(14,*) conan(tk),conan2(tk),conan(tk)+conan2(tk) 
        STOP 
        END
c Now compute the energy density and the pressure

!$$$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,a,Y, DYDU,adot, k,Omega2,m,lamda,Eulergamma 
        Real*8 T,Td,TRE,reldifft,reldiff1,pi,Gu,Traan
        Real*8 a1,a2,a3,a4,Hu1,Hu3,ad1,ad2,ad3,ad4,cfit 
	PARAMETER(NN=60000,nt=15000)
	DIMENSION Y(NVAR),DYDU(NVAR),k(nt),T(nt),Td(nt)
        DIMENSION reldifft(nt),cfit(15,3)
	Common/param/m,k,reldiff1,pi,lamda,Eulergamma
        Common /cinfo/cfit
!$$$    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
c        a1 = -Y(1)**2+lamda*Y(1)**4/3.d0
c this is actually aprime^2
c        a2 = -Y(1)+2.d0*lamda*Y(1)**3/3.d0
c        a3 = -a1+2.d0*lamda*Y(1)**2*a1
c this is actually ap1*ap3
c        a4 = Y(1)-20.d0*lamda*Y(1)**3/3.d0
c     +  +8.d0*lamda**2*Y(1)**5/3.d0
        a1 = a(u)*ad1(u)
        a2 = a(u)*ad1(u)**2+a(u)**2*ad2(u)
        a3 = a(u)*ad1(u)**3+4.d0*a(u)**2*ad1(u)*ad2(u)
     +  +a(u)**3*ad3(u)
        a4 = a(u)*ad1(u)**4+11.d0*a(u)**2*ad1(u)**2*ad2(u)
     +  +4.d0*a(u)**3*ad2(u)**2+7.d0*a(u)**3*ad1(u)*ad3(u)
     +  +a(u)**4*ad4(u)
        Hu1=-36.d0*a4/a(u)**5+144.d0*a3*a1/a(u)**6
     -  -216.d0*a2*a1**2/a(u)**7+108.d0*a2**2/a(u)**6
     +  +72.d0*a2/a(u)**5-72.d0*a1**2/a(u)**6

        Hu3=12.d0*a2*a1**2/a(u)**7-12.d0*a1**4/a(u)**8
     +  +12.d0*a2/a(u)**5-12.d0*a1**2/a(u)**6
c        Hu1=-36.d0*a4/Y(1)**5+144.d0*a3/Y(1)**6
c     -  -216.d0*a2*a1/Y(1)**7+108.d0*a2**2/Y(1)**6
c     +  +72.d0*a2/Y(1)**5-72.d0*a1/Y(1)**6

c        Hu3=12.d0*a2*a1/Y(1)**7-12.d0*a1**2/Y(1)**8
c     +  +12.d0*a2/Y(1)**5-12.d0*a1/Y(1)**6
c        H01(it) =-36.d0*a3(it)/Y(1)**6+72.d0*a2(it)*a1(it)
c     /  /Y(1)**7+18.d0*a2(it)**2/Y(1)**6+36.d0*a1(it)/Y(1)**6
c     -  -18.d0/Y(1)**4


c        H03(it) = 3.d0*a1(it)**2/Y(1)**8+6.d0*a1(it)/Y(1)**6
c     +  +3.d0/Y(1)**4


        Tre = 0
        DO i = 3,nvar-3,4
        n = n+1


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  
        T(n)=k(n)**2*(Y(i)**2+Y(i+2)**2)*m**2/(2.d0*pi**2*a(u)**3)

        Td(n)=(k(n)*m**2*a(u)**2-m**4*a(u)**4/(2.d0*k(n)))
     /  /(4*pi**2*a(u)**4)
        reldifft(n) = Abs((T(n)-Td(n))/Td(n))
        if(reldifft(n)<reldiff1) exit

        Tre = Tre + T(n) - Td(n)

        Enddo

c        Gu = -6.d0*a2/Y(1)**3-6.d0/Y(1)**2
        Gu = -6.d0*a2/a(u)**3-6.d0/a(u)**2

c        Traan =m**2*Gu/(288.d0*pi**2)-m**4*(1.d0
        Traan=(-Hu1/6.d0+Hu3)/(2880.d0*pi**2) 
     +  +m**2*Gu/(288.d0*pi**2)-m**4*(1.d0
     +  +DlOG(m**2*a(u)**2/4.d0)+2.d0*Eulergamma)/(16.d0*pi**2)




        DYDU(1) = Y(2)
        DYDU(2) = -1.d0/Y(1)-Y(2)**2/Y(1)+4.d0*pi*(Tre+Traan)*Y(1)/3.d0
     +  +2.d0*lamda*Y(1)/3.d0

        n = 0
        DO  i=3,NVAR-3,4
           n = n+1
           Omega2=-0.5d0*ad2(u)/a(u)+0.25d0*(ad1(u)/a(u))**2+k(n)**2/a(u)**2+m**2
c (k(n)**2-0.25d0)/a(u)**2 + m**2-0.25d0
           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/a(u)
	Return
	End
        
        

c  define the scaler factor.
c in this case, it is de sitter space

        double precision function a(u)
        real*8 u,cfit
        Dimension cfit(15,3)
        Common /cinfo/ cfit 
        a=cfit(1,1)+cfit(2,1)*u+cfit(3,1)*u**2+cfit(4,1)*u**3
     +  +cfit(5,1)*u**4+cfit(6,1)*u**5+cfit(7,1)*u**6+cfit(8,1)*u**7
     +  +cfit(9,1)*u**8+cfit(10,1)*u**9+cfit(11,1)*u**10
     +  +cfit(12,1)*u**11+cfit(13,1)*u**12+cfit(14,1)*u**13
     +  +cfit(15,1)*u**14
        return
        end

        double precision function ad1(u)
        real*8 u,cfit
        Dimension cfit(15,3)
        Common /cinfo/ cfit
        ad1=cfit(1,2)+cfit(2,2)*u+cfit(3,2)*u**2+cfit(4,2)*u**3
     +  +cfit(5,2)*u**4+cfit(6,2)*u**5+cfit(7,2)*u**6+cfit(8,2)*u**7
     +  +cfit(9,2)*u**8+cfit(10,2)*u**9+cfit(11,2)*u**10
     +  +cfit(12,2)*u**11+cfit(13,2)*u**12+cfit(14,2)*u**13
     +  +cfit(15,2)*u**14        
        return 
        end

        double precision function ad2(u)
        real*8 u,cfit
        Dimension cfit(15,3)
        Common /cinfo/ cfit
        ad2 = cfit(1,3)+cfit(2,3)*u+cfit(3,3)*u**2+cfit(4,3)*u**3
     +  +cfit(5,3)*u**4+cfit(6,3)*u**5+cfit(7,3)*u**6+cfit(8,3)*u**7
     +  +cfit(9,3)*u**8+cfit(10,3)*u**9+cfit(11,3)*u**10 
     +  +cfit(12,3)*u**11+cfit(13,3)*u**12+cfit(14,3)*u**13
     +  +cfit(15,3)*u**14
        return
        end

        double precision function ad3(u)
        real*8 u,cfit
        Dimension cfit(15,3)
        Common /cinfo/ cfit 
        ad3 = cfit(2,3)+2.d0*cfit(3,3)*u+3.d0*cfit(4,3)*u**2
     +  +4.d0*cfit(5,3)*u**3+5.d0*cfit(6,3)*u**4
     +  +6.d0*cfit(7,3)*u**5+7.d0*cfit(8,3)*u**6
     +  +8.d0*cfit(9,3)*u**7+9.d0*cfit(10,3)*u**8+10.d0*cfit(11,3)*u**9
     +  +11.d0*cfit(12,3)*u**10+12.d0*cfit(13,3)*u**11
     +  +13.d0*cfit(14,3)*u**12+14.d0*cfit(15,3)*u**13
        return
        end

        double precision function ad4(u)
        
        real*8 u,cfit
        Dimension cfit(15,3)
        Common /cinfo/ cfit
        ad4 = 2.d0*cfit(3,3)+6.d0*cfit(4,3)*u
     +  +12.d0*cfit(5,3)*u**2+20.d0*cfit(6,3)*u**3+30.d0*cfit(7,3)*u**4
     +  +42.d0*cfit(8,3)*u**5
     +  +56.d0*cfit(9,3)*u**6+72.d0*cfit(10,3)*u**7
     +  +90.d0*cfit(11,3)*u**8
     +  +110.d0*cfit(12,3)*u**9+132.d0*cfit(13,3)*u**10
     +  +156.d0*cfit(14,3)*u**11+182.d0*cfit(15,3)*u**12

        return
        end