!-----------------------------------------------------------------------------
!
! Linear Growth Factor, g(z)=D(z)(1+z), which is constant and equal to 
! unity during the matter dominated era. Also computes its derivative 
! with respect to lna, dg/dlna.
! 
! *** Radiation density is ignored in the calculation ***
!
! << sample code >>
!
! USE cosmo
! USE growth
! double precision :: g,dgdlna,redshift
! external g,dgdlna
! ! Specify three cosmological parameters in double precision.
! ! The data type has been defined in MODULE cosmo.
! ode0=0.723d0
! om0=0.277d0
! w   = -1d0
! CALL setup_growth
! redshift=10d0
! print*,'omega matter=',om0
! print*,'omega de=',ode0
! print*,'w=',w
! print*,'g(z) at z=',redshift,'is',g(redshift)
! print*,'dg/dlna(z) at z=',redshift,'is',dgdlna(redshift)
! end
!
! August 23, 2008
! E. Komatsu
!
!-----------------------------------------------------------------------------
MODULE Growth
  INTEGER, parameter :: nw=501
  DOUBLE PRECISION, dimension(nw) :: xa,ga,g2a,dga,dg2a
contains
  SUBROUTINE Setup_Growth
    ! tabulate g(x), where x=ln(a)
    IMPLICIT none
    INTEGER :: i,ind
    DOUBLE PRECISION :: x,xini,y(2),c(24),work(2,9),tol=1.d-8
    EXTERNAL fderiv
    ind=1 
    do i=1,nw
       x = (i-nw)*0.01 ! x=ln(a)=-5.00,-4.99,...,0
       xa(i) = x
       if(i.eq.1)then
          xini = xa(i)
          y(1)=1d0 ! initial condition, g=1, at z=-5.00 (z=147.413)
          y(2)=0d0 ! initial condition, g'=1, at z=-5.00 (z=147.413)
       else
          call dverk(2,fderiv,xini,y,x,tol,ind,c,2,work)
          xini = x
       endif
       ga(i) = y(1)
       dga(i)= y(2)
    enddo
    CALL spline(xa,ga,nw,1.d30,1.d30,g2a)  ! evaluate g2a(i)=d^2[g(i)]/dy^2
    CALL spline(xa,dga,nw,1.d30,1.d30,dg2a)! evaluate g2a(i)=d^2[dg/dlna(i)]/dy^2
    return
  END SUBROUTINE Setup_Growth
END MODULE Growth
!-----------------------------------------------------------------------------
SUBROUTINE fderiv(n,x,y,yprime)
  USE cosmo
! x = ln(a) = -ln(1+z)
! y(1) = g
! y(2) = g'
! See Eq.(76) of Komatsu et al. (2008) [WMAP 5-year]
  IMPLICIT none 
  INTEGER, intent(IN) :: n
  DOUBLE PRECISION :: x,y(n),yprime(n)
  DOUBLE PRECISION :: ok,ode,ok0,h2,z,a
  DOUBLE PRECISION, parameter :: or0 = 0d0 ! ignore radiation density
  z=dexp(-x)-1d0
  a=dexp(x)
  ok0=1d0-om0-ode0
  h2=om0/a**3d0+or0/a**4d0+ok0/a**2d0+ode0/a**(3d0*(1.+w))
  ok=ok0/h2*dexp(-2d0*x)
  ode=ode0/h2*dexp(-3d0*x*(1d0+w)) 
  yprime(1)= y(2) 
  yprime(2)= -(2.5d0+0.5d0*(ok-3d0*w*ode))*y(2) &
             -(2d0*ok+1.5d0*(1d0-w)*ode)*y(1)
  return
END SUBROUTINE fderiv
!-----------------------------------------------------------------------------
DOUBLE PRECISION FUNCTION g(z)
  USE growth
  ! Returns g(z) = D(z)(1+z) by interpolating the table generated by
  ! Setup_Growth. 
  IMPLICIT none
  DOUBLE PRECISION, intent(IN) :: z
  DOUBLE PRECISION :: a,b,hh,x
  INTEGER :: jlo
  x=-dlog(1d0+z) ! x=ln(a)=-ln(1+z)
  CALL hunt(xa,nw,x,jlo)
  hh=xa(jlo+1)-xa(jlo)
  a=(xa(jlo+1)-x)/hh
  b=(x-xa(jlo))/hh
  g=a*ga(jlo)+b*ga(jlo+1)+((a**3-a)*g2a(jlo)+(b**3-b)*g2a(jlo+1))*(hh**2)/6.
  return
END FUNCTION g
!-----------------------------------------------------------------------------
DOUBLE PRECISION FUNCTION dgdlna(z)
  USE growth
  ! Returns dg/dlna(z) by interpolating the table generated by
  ! Setup_Growth.
  IMPLICIT none
  DOUBLE PRECISION, intent(IN) :: z
  DOUBLE PRECISION :: a,b,hh,x
  INTEGER :: jlo
  x=-dlog(1d0+z) ! x=ln(a)=-ln(1+z)
  CALL hunt(xa,nw,x,jlo)
  hh=xa(jlo+1)-xa(jlo)
  a=(xa(jlo+1)-x)/hh
  b=(x-xa(jlo))/hh
  dgdlna=a*dga(jlo)+b*dga(jlo+1)+((a**3-a)*dg2a(jlo)+(b**3-b)*dg2a(jlo+1))*(hh**2)/6.
  return
END FUNCTION dgdlna