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