! This program takes the values of the Bunch-Davies modes at u = -1 and
! uses them as initial conditions for the mode equation. The mode equation is 
! then solved in terms of proper time u to get the exact Bunch-Davies modes for all times.

c  I've decided for now the best way to do it is to have different versions.  This version
c  is for m = 10 and k = 1.  There will be 6 combinations total.
	Implicit None
	Integer i,im, nu, it, nk, nbad,nok,junk
	real*8 ustart,uin,u,delu,m,uout,y(4)
	Real*8 pi,eps,h1,hmin,k
	Real*8 gamma
	Real*8 B2k,Ckr,Cki,vr,vpr,vi,vpi
	common/intp/k
	common/mass/m
        EXTERNAL DERIVS
        EXTERNAL BSSTEP
	Open(unit = 9, file = "BDmodes_v.dat")
	open(unit = 10, file = "BDmodes_v.out")
	open(unit = 11, file = "C_m_10_k_1.out")
	open(unit = 12, file = "B_m_10_k_1.out")
	read(9,*) junk
	read(9,*) uin,delu,nu
	read(9,*) junk
	read(9,*),eps,h1,hmin
	pi = 2.d0*dasin(1.d0)
	Print *,'pi = ',pi
c  Start with m = k = 1
	m = 1.d1
        gamma = dsqrt(m**2 - 0.25d0)
	k = 1.d0
c  We start the integration at u = -1
	  ustart = -1.d0
c  Y(1) is vr, Y(2) is vpr, Y(3) is vi, Y(4) is vpi
	  y(1) = -0.09683147691231392d0
	  y(2) = -0.7609130766547556d0
	  y(3) = -0.06502026338505628d0
	  y(4) = 0.8944231520632057d0
          CALL ODEINT(Y,4,ustart,uin,EPS,H1,HMIN,NOK,NBAD,DERIVS,BSSTEP)
	  u = uin
	  uout = uin + delu
	  vr = Y(1)
	  vpr = Y(2)
	  vi = Y(3)
	  vpi = Y(4)
	  write(10,*) k, u, y(1), y(2), y(3), y(4)  
c  Now compute  the coefficients of Ckr and Cki
	  Ckr = -(k**2*vi**2)/(8.d0*Pi**2) - (gamma**2*k**2*vi**2)/(2.d0*Pi**2) - 
     -  (k**2*vpi**2)/(2.d0*Pi**2) + (k**2*vpr**2)/(2.d0*Pi**2) + 
     -  (k**2*vr**2)/(8.d0*Pi**2) + (gamma**2*k**2*vr**2)/(2.d0*Pi**2) - 
     -  (k**4*vi**2)/(2.d0*Pi**2*dcosh(u)**2) + (k**4*vr**2)/(2.d0*Pi**2*dcosh(u)**2) - 
     -  (k**2*vi*vpi*dtanh(u))/Pi**2 + (k**2*vpr*vr*dtanh(u))/Pi**2 - 
     -  (k**2*vi**2*dtanh(u)**2)/(2.d0*Pi**2) + (k**2*vr**2*dtanh(u)**2)/(2.d0*Pi**2)
	  Cki = -((k**2*vpi*vpr)/Pi**2) - (k**2*vi*vr)/(4.d0*Pi**2) - 
     -  (gamma**2*k**2*vi*vr)/Pi**2 - (k**4*vi*vr)/(Pi**2*dcosh(u)**2) - 
     -  (k**2*vi*vpr*dtanh(u))/Pi**2 - (k**2*vpi*vr*dtanh(u))/Pi**2 - 
     -  (k**2*vi*vr*dtanh(u)**2)/Pi**2
	  Write(11,29) u,k,Ckr,Cki
          B2k = (k**2*vi**2)/(8.d0*Pi**2) + (gamma**2*k**2*vi**2)/(2.d0*Pi**2) +
     -  (k**2*vr**2)/(8.d0*Pi**2) + (gamma**2*k**2*vr**2)/(2.d0*Pi**2) +
     -  (k**4*vi**2/dcosh(u)**2)/(2.d0*Pi**2) +
     -  (k**4*vr**2/dcosh(u)**2)/(2.d0*Pi**2) +
     -  (k**2*vi**2*dtanh(u)**2)/(2.d0*Pi**2) + (k**2*vr**2*dtanh(u)**2)/(2.d0*Pi**2)
     -  + k**2*(vpr**2+vpi**2 + 2.d0*dtanh(u)*(vr*vpr+vi*vpi))/(2.d0*Pi**2)

          Write(12,29) u,k,B2k
	  do it = 1,nu
            CALL ODEINT(Y,4,u,uout,EPS,H1,HMIN,NOK,NBAD,DERIVS,BSSTEP)
	    write(10,*) k, uout, y(1), y(2), y(3), y(4)
	    u = uout
	    uout = uout + delu
	    vr = Y(1)
	    vpr = Y(2)
	    vi = Y(3)
	    vpi = Y(4)
c  Now compute  the coefficients of Ckr and Cki
	    Ckr = -(k**2*vi**2)/(8.d0*Pi**2) - (gamma**2*k**2*vi**2)/(2.d0*Pi**2) - 
     -  (k**2*vpi**2)/(2.d0*Pi**2) + (k**2*vpr**2)/(2.d0*Pi**2) + 
     -  (k**2*vr**2)/(8.d0*Pi**2) + (gamma**2*k**2*vr**2)/(2.d0*Pi**2) - 
     -  (k**4*vi**2)/(2.d0*Pi**2*dcosh(u)**2) + (k**4*vr**2)/(2.d0*Pi**2*dcosh(u)**2) - 
     -  (k**2*vi*vpi*dtanh(u))/Pi**2 + (k**2*vpr*vr*dtanh(u))/Pi**2 - 
     -  (k**2*vi**2*dtanh(u)**2)/(2.d0*Pi**2) + (k**2*vr**2*dtanh(u)**2)/(2.d0*Pi**2)
	    Cki = -((k**2*vpi*vpr)/Pi**2) - (k**2*vi*vr)/(4.d0*Pi**2) - 
     -  (gamma**2*k**2*vi*vr)/Pi**2 - (k**4*vi*vr)/(Pi**2*dcosh(u)**2) - 
     -  (k**2*vi*vpr*dtanh(u))/Pi**2 - (k**2*vpi*vr*dtanh(u))/Pi**2 - 
     -  (k**2*vi*vr*dtanh(u)**2)/Pi**2
	    Write(11,29) u,k,Ckr,Cki
	    B2k = (k**2*vi**2)/(8.d0*Pi**2) + (gamma**2*k**2*vi**2)/(2.d0*Pi**2) +
     -  (k**2*vr**2)/(8.d0*Pi**2) + (gamma**2*k**2*vr**2)/(2.d0*Pi**2) +
     -  (k**4*vi**2/dcosh(u)**2)/(2.d0*Pi**2) +
     -  (k**4*vr**2/dcosh(u)**2)/(2.d0*Pi**2) +
     -  (k**2*vi**2*dtanh(u)**2)/(2.d0*Pi**2) + (k**2*vr**2*dtanh(u)**2)/(2.d0*Pi**2)
     -  + k**2*(vpr**2+vpi**2 + 2.d0*dtanh(u)*(vr*vpr+vi*vpi))/(2.d0*Pi**2)
            Write(12,29) u,k,B2k

	  end do
   29   Format(1x,20(d12.5,2x))
	stop
	end

 	SUBROUTINE DERIVS(u,Y,DYDT,NVAR)
        INTEGER nvar
        REAL*8 u, Y(4), DYDT(4), k, m, a, h
        COMMON/IntP/k
        COMMON/mass/m
        a = dcosh(u)
        h = dtanh(u)
        DYDT(1) = Y(2)
        DYDT(2) = -3.d0*h*Y(2)-((k**2-1.d0)/dcosh(u)**2 + (m**2+2.d0))*Y(1)
        DYDT(3) = Y(4)
        DYDT(4) = -3.d0*h*Y(4)-((k**2-1.d0)/dcosh(u)**2 + (m**2+2.d0))*Y(3)
        RETURN
        END