!!The latest version is now maintained at 
!!   http://www.cita.utoronto.ca/~zqhuang/QCDM

!!================= F90 Module to calculate quintessence + CDM cosmology (when quintessence is minimally coupled to gravity);
!! written by Zhiqi Huang (zhiqihuang.cea.fr; zqhuang@cita.utoronto.ca)
!! Please cite Huang, Bond, Kofman, 2011 (arXiv: 1007.5297) if you use this code
!! Assuming 3 species of massless neutrinos and Omega_k = 0; Please contact me for an extension of this code that includes Omega_k.
!! This module does not evaluate D(z) (growth factor); Please contact me if you want an extension of this code that includes D(z) calculation.
!! **********************************************************************************************************
!! Conventions:
!! z: redshift
!! V(\phi): the quintessence potential;  
!! G_N: Newton's gravitational constant;
!! \rho_{m}: the total matter density (including cold dark matter and baryon); it is a function of redshift;
!! \dot\phi: d\phi/dt, where t is the cosmological time;
!! \rho_{\phi}: the energy density of quintessence; ( = 1/2 \dot\phi ^2 + V(\phi))
!! p_{\phi}: the pressure of quintessence; (= 1/2 \dot\phi ^2 - V(\phi) )
!!  w_{\phi}: the equation of state of quintessence; ( = p_{\phi} / \rho_{\phi} )
!!  a: scale factor 1/(1+z) (i.e. normalized to 1 today);
!!  H_0 = 100 h km/s/Mpc: current Hubble parameter
!!**********************************************************************************************************
!! The parameters:
!!==========================================================================================================
!! Omega_m0:
!! $\Omega_{m0}$ is the fraction of matter at redshift zero (including cold dark matter and baryon).
!!    \Omega_{m0}   =    \rho_m  / ( 3 H^2 / 8\pi G_N)    evaluated at z=0
!!===========================================================================================================
!! epsilon_s:
!!  $\epsilon_s$ controls the slope of the quintessence potential at low redshift where DE is relevant; 
!!  Its exact definition is:
!!    \epsilon_s = 1/(16\pi G_N) (d \ln V / d\phi )^2  evaluated where \rho_{\phi} = \rho_m
!!  For quintessence \epsilon_s >= 0 (by definition); A recommended prior is 0 <= \epsilon_s <= 2  (8-sigma upperbound using current SN+CMB data); 
!!  If you want to choose reasonable fiducial \epsilon_s to run Nbody simulation, you can use the current one-sigma bound epsilon_s = 0.25.
!!  Note: using a negative \epsilon_s won't cause the code crash; In this case you are actually doing phantom models (see HBK2011).
!!==========================================================================================================
!! epsilon_infty
!!  $\epsilon_\infty$ controls the high-redshift attractor for tracking models.
!!  Its exact definition is
!!    \epsilon_\infty = | 1 + w_{\phi) | / (1 + w_m) evaluated at z>>1; 
!!    where w_m is 1/3 in the radiation dominated regime and 0 in the matter dominated regime.
!!  This limit exists and is nonzero for tracking models (by definition);
!!  For thawing models, this parameter is zero (w_{\phi} = -1 at high redshift, because the quintessence field does not move in early universe).
!!  In order not to spoil the BBN and CMB physics, you can use a prior 0 <= \epsilon_\infty < 1; 
!!  If you want to choose a reasonable fiducial value, use \epsilon_\infty ~ 0.1 (negative power-law) or \epsilon_\infty = 0 (all thawing models);
!!============================================================================================================
!! zeta_s
!!   For thawing models, this parameter controls the second derivative of $\ln V(\phi)$ at low redshift; For tracking models, it is a derived parameter that depends on both the second derivative of $\ln V(\phi)$ and $\dot \phi$ at low redshift. See HBK2011 for the exact definition.
!!   For most quintessence models, |\zeta_s|<1; A recommended prior is -1<\zeta_s<1;
!!   This parameter will only give ~0.01 corrections to w(a); If you are choosing a fiducial model, you can just use \zeta_s = 0. 
!!*****************************************************************************************************************************

!! Instruction continues on line 240.

!! This module does the initialization of QCDM cosmology; w_\phi, \rho_\phi, and comoving distance are tabulated here.
!! If you do not care about the details (how these quantities are calculated), just jump to line 240 and read the instruction.
Module QCDMInit 
  implicit none
  public
  private::dl,TBL_SIZE, wphiCosmology, QCDM_initialized

  integer,parameter::dl=Kind(1.d0) 
  integer,parameter::TBL_SIZE = 8192

  Type WphiCosmology
     real(dl)  h   !! Hubble constant H_0 = 100 h km/s/Mpc
     real(dl)  Omega_m0 !! nonrealativistic matter energy fraction at redshift zero
     !! w parameters ==>
     real(dl)  epsilon_s !!\epsilon_s, negative for phantom, positive for quintessence; 
     real(dl)  epsilon_infty !!|\epsilon_{\infty\phi}|/\epsilon_m, always positive; 0<=epsilon_infty<1
     real(dl)  zeta_s !! correction from the second derivative of lnV at low redshift; -1<zeta_s<1
     !!derived parameters ==>
     real(dl) qpsign !! 1 for quintessence (epsilon_s>0); -1 for phantom (epsilon_s<0)
     real(dl) a_eq !!the scale factor at DE-matter equality
     real(dl) Omega_r !! effective radiation energy fraction at readshift zero
     real(dl) Omega_L  !! dark energy energy fraction at redshift zero
     real(dl),dimension(TBL_SIZE)::x, ln_rho_DE, chi
     !!the table of \ln (\rho_{DE}(x)/\rho_{DE,0})
     !! where x:= -\ln a = \ln (1+z)
     !! \chi is in unit Mpc
     real(dl) xstep  !! x(2) - x(1)
  End type WphiCosmology


  Type(WphiCosmology) WPhiCosmo
  Logical:: QCDM_initialized = .false.


contains

  !!Initialization, this subroutine is not thread safe, please protect it with omp critical when using openmp
  subroutine QCDM_initialize(h, Omega_m0, epsilon_s, epsilon_infty, zeta_s)
    real(dl) h, Omega_m0, epsilon_s, epsilon_infty, zeta_s !!input parameters
    real(dl) delta, l0,l1,l2, eps_inf,zmax,xmax, w(TBL_SIZE), a(TBL_SIZE)
    real(dl),parameter::sqrt2 = 1.414213562373_dl
    integer i


    !! now prepare for w(z) calculation:
    if(Omega_m0 .lt. 0.1 .or. Omega_m0 .gt. 0.5 .or. h.lt.0.4 .or. h.gt. 1.) &
         stop "bad input in QCDM: I am assuming 0.1 < Omega_m0<0.5 and 0.4<h<1."
    if(epsilon_s.gt.2. .or. epsilon_s<-2.) &
         stop "bad input in QCDM: |epsilon_s| > 2; the w_phi approximation will be inaccurate"
    if(epsilon_infty.lt.0. .or. epsilon_infty .gt.1) &
         stop "bad input in QCDM: this module only works for 0<=epsilon_infty<1"
    if(zeta_s .gt. 2 .or. zeta_s .lt. -2) &
         stop "bad input in QCDM: |zeta_s|>2; The correction parameter should be < O(1) (light quintessence field)"
    if(QCDM_initialized)then
       if(       abs(WphiCosmo%Omega_m0 - Omega_m0) .lt. 1.e-5 &
            .and. abs(wphiCosmo%epsilon_s - epsilon_s) .lt. 1.e-5 &
            .and. abs(wphiCosmo%epsilon_infty - epsilon_infty) .lt. 1.e-5 &
            .and. abs(wphiCosmo%zeta_s - zeta_s) .lt. 1.e-5 ) then
          if(abs(wphiCosmo%h -h) .lt. 1.e-5 )then
             return    !!do not need to do the initialization
          else
             wphiCosmo%chi = wphiCosmo%chi*(wphiCosmo%h/h)
             wphiCosmo%h = h
             return
          endif
       endif
    else
       QCDM_initialized = .true.
    endif
    !!copy parameters over
    WPhiCosmo%h = h
    WPhiCosmo%Omega_m0 = Omega_m0
    WPhiCosmo%epsilon_s = epsilon_s
    WPhiCosmo%epsilon_infty = epsilon_infty
    WPhiCosmo%zeta_s = zeta_s
    if(epsilon_s.lt.0.)then
       WPhiCosmo%qpsign = - 1.
    else
       WPhiCosmo%qpsign = 1.
    endif

    !!calculate a_{eq}
    WPhiCosmo%Omega_r = 4.1814d-5/h**2 !!assuming three spieces of light neutrinos
    WPhiCosmo%Omega_L = 1. - WPhiCosmo%Omega_m0 - WPhiCosmo%Omega_r
    eps_inf = 1.5_dl * epsilon_infty
    l0 = sqrt(eps_inf)
    l1 = sqrt(abs(epsilon_s))-sqrt2 * l0

    delta=(l0 + (0.533 - 0.1 * zeta_s)*l1 )**2 &
         + (l0 + ((0.91-0.78*Omega_m0) + (0.236-0.76*Omega_m0)*zeta_s)*l1 )**2

    WPhiCosmo%a_eq = (Omega_m0/WPhiCosmo%Omega_L)**(1./(3.- WphiCosmo%qpsign * delta))

    !!define the x table:
    zmax = 20._dl * WPhiCosmo%Omega_m0/WPhiCosmo%Omega_r !!for z>zmax, w~const
    xmax = log ( 1. + zmax )  !! ~ 10
    WPhiCosmo%xstep = xmax/(TBL_SIZE-1._dl)
    do i=1,TBL_SIZE
       WPhiCosmo%x(i) = xmax * (i-1._dl)/(TBL_SIZE - 1._dl)
    enddo
    a =  exp( - WPhiCosmo%x)
    do i=1, TBL_SIZE
       w(i) = wphi_of_a ( a(i) )
    enddo

    !! now tabulate ln_rho_DE(x) := \ln [\rho_{DE}(z) / \rho_{DE,0}]
    WPhiCosmo%ln_rho_DE(1) = 0. !!
    do i=2, TBL_SIZE
       WPhiCosmo%ln_rho_DE(i) = WPhiCosmo%ln_rho_DE(i-1) + (6. + w(i-1)+w(i)+4.*wphi_of_a(exp(- (WPhiCosmo%x(i)+WPhiCosmo%x(i-1))/2.)))
    enddo
    WPhiCosmo%ln_rho_DE = WPhiCosmo%ln_rho_DE * WPhiCosmo%xstep/2.

    !! now tabulate chi := \chi(x) in Mpc
    WPhiCosmo%chi(1) = 0.
    do i=2, TBL_SIZE
       WPhiCOsmo%chi(i) = WPhiCosmo%chi(i-1) + 1./Ha_ratio_of_a( a(i) ) + 1./Ha_ratio_of_a( a(i-1) ) + 4./Ha_ratio_of_a( exp(- (WPhiCosmo%x(i)+WPhiCosmo%x(i-1))/2.) )
    enddo
    WPhiCosmo%chi = WPhiCosmo%chi/6.*WPhiCosmo%xstep * 2997.92458/WPhiCosmo%h
    return
  contains
    function Ha_ratio_of_a(a) !!return  a * H(a) / H_0
      real(dl) a, Ha_ratio_of_a
      Ha_ratio_of_a = sqrt((WPhiCosmo%Omega_m0+WPhiCosmo%Omega_r/a) /a + WphiCosmo%Omega_L*rho_DE_ratio_of_a(a)*a**2)
    end function Ha_ratio_of_a

  end subroutine QCDM_initialize

  Function wphi_of_a (a)
    !! return w_\phi as a function of a
    !! thread safe
    real(dl) a, wphi_of_a
    real(dl) mu, mu3, s0, s1, Fx, F2x, eps_inf
    real(dl),parameter::sqrt2 = 1.414213562373_dl
    eps_inf = WPhiCosmo%epsilon_infty * (1. + WPhiCosmo%Omega_r/3./(WPhiCosmo%Omega_r + WPhiCosmo%Omega_m0 * a))*1.5
    mu=a/WPhiCosmo%a_eq
    if(mu .lt. 1.e-4 )then  !!high redshift  
       wphi_of_a = -1. + WPhiCosmo%qpsign/1.5 * eps_inf
       return
    endif
    mu3 = mu**3
    s0 = sqrt(mu3)
    s1 = sqrt(1.+mu3)
    Fx = s1/s0 - log( s0 + s1 )/mu3
    F2x = sqrt2 * (1. - log(1.+mu3)/mu3) - Fx
    s0 = sqrt(eps_inf)
    wphi_of_a = -1. + WPhiCosmo%qpsign/1.5 * ( &
         s0 + (sqrt(abs(WPhiCosmo%epsilon_s)) - sqrt2 * s0 )*(Fx + WPhiCosmo%zeta_s * F2x) &
         )**2
    return
  end Function wphi_of_a

  Function rho_DE_ratio_of_a(a) 
    !!return \rho_\phi (a) / \rho_\phi (z=0)
    !!thread safe
    real(dl) a, rho_DE_ratio_of_a, x, pos
    integer i
    x = -log(max(a, 1.e-10_dl))
    pos = x/WPhiCosmo%xstep + 1._dl
    i =  floor(pos) 
    If(i .ge. TBL_SIZE)then
       rho_DE_ratio_of_a = exp(WPhiCosmo%ln_rho_DE(TBL_SIZE) + 4. *(x - WPhiCosmo%x(TBL_SIZE)) * WPhiCosmo%epsilon_infty * WPhiCosmo%qpsign )
       return
    endif
    x = pos - i
    rho_DE_ratio_of_a = exp(WPhiCosmo%ln_rho_DE(i) * (1.-x) + WPhiCosmo%ln_rho_DE(i+1)*x)
    return
  end Function rho_DE_ratio_of_a

  Function comoving_distance_of_a(a) 
    !!return comoving distance in unit Mpc; it equals to the comoving angular diameter distance 
    !!thread safe
    real(dl) a, comoving_distance_of_a , x , pos
    integer i
    x = -log(max(a, 1.e-10_dl))
    pos = x/WPhiCosmo%xstep + 1._dl
    i =  floor(pos) 
    if(i.ge.TBL_SIZE)then
       comoving_distance_of_a = WPhiCosmo%chi(TBL_SIZE)
       return
    endif
    x = pos - i
    comoving_distance_of_a = WPhiCosmo%chi(i) * (1.-x) + WPhiCosmo%chi(i+1)*x
    return
  end function comoving_distance_of_a

End Module QCDMInit


!!********************************************************************************************************************
!!This module is a wrapped blackbox for you to use. 
!!It is NOT thread safe. Do NOT calculate different cosmology simultaneously in multiple threads.

!!All functions have six arguments: 
!!   the first five define the cosmology: h, Omega_m0, epsilon_s, epsilon_infty, zeta_s
!!   the last argument is redshift z
!!   all variables are double precision
!! To understand the meaning of the parameters, read the instructions at the beginning of this file.

!!The list of functions are:
!!QCDM_w_phi (h, Omega_m0, epsilon_s, epsilon_infty, zeta_s, z)   
!!          returns the equation of state of quintessence (or phantom if epsilon_s<0) at redshift z
!!QCDM_rho_phi_ratio (h, Omega_m0, epsilon_s, epsilon_infty, zeta_s, z)   
!!          returns \rho_\phi (z) /  \rho_\phi(0); i.e. the ratio between quintessence energy densities at redshift z and redshift 0
!!QCDM_Hubble (h, Omega_m0, epsilon_s, epsilon_infty, zeta_s, z)
!!          returns the Hubble expansion rate in unit of km/s/Mpc at redshift z 
!!QCDM_angular_diameter_distance (h, Omega_m0, epsilon_s, epsilon_infty, zeta_s, z)
!!          returns the physical angular diameter distance in unit of Mpc
!!QCDM_luminosity_distance (h, Omega_m0, epsilon_s, epsilon_infty, zeta_s, z)
!!          returns the luminosity distance in unit of Mpc

Module QCDM
  use QCDMInit
  implicit none
  public
  private:: dl
  integer,parameter::dl=Kind(1.d0) !!change to Kind(1.0) for single precision

contains
  
  function QCDM_w_phi (h, Omega_m0, epsilon_s, epsilon_infty, zeta_s, z)
    real(dl) QCDM_w_phi
    real(dl) h, Omega_m0, epsilon_s, epsilon_infty, zeta_s, z
    call QCDM_Initialize(h, Omega_m0, epsilon_s, epsilon_infty, zeta_s)
    QCDM_w_phi = wphi_of_a( 1./(1.+z))
  end function QCDM_w_phi

  function QCDM_rho_phi_ratio(h, Omega_m0, epsilon_s, epsilon_infty, zeta_s, z)
    real(dl) QCDM_rho_phi_ratio
    real(dl) h, Omega_m0, epsilon_s, epsilon_infty, zeta_s, z
    call QCDM_Initialize(h, Omega_m0, epsilon_s, epsilon_infty, zeta_s)
    QCDM_rho_phi_ratio =  rho_DE_ratio_of_a(1./(1.+z))
  end function QCDM_rho_phi_ratio

  function QCDM_Hubble (h, Omega_m0, epsilon_s, epsilon_infty, zeta_s, z)
    real(dl) QCDM_Hubble
    real(dl) h, Omega_m0, epsilon_s, epsilon_infty, zeta_s, z
    call QCDM_Initialize(h, Omega_m0, epsilon_s, epsilon_infty, zeta_s)
    QCDM_Hubble = wphiCosmo%h*100.*sqrt(wphiCosmo%Omega_L * rho_DE_ratio_of_a(1./(1.+z)) & 
         + (wphiCosmo%Omega_r*(1.+z)+wphiCosmo%Omega_m0)*(1.+z)**3)
  end function QCDM_Hubble

  function QCDM_angular_diameter_distance (h, Omega_m0, epsilon_s, epsilon_infty, zeta_s, z)
    real(dl) QCDM_angular_diameter_distance
    real(dl) h, Omega_m0, epsilon_s, epsilon_infty, zeta_s, z
    call QCDM_Initialize(h, Omega_m0, epsilon_s, epsilon_infty, zeta_s)
    QCDM_angular_diameter_distance = Comoving_Distance_of_a(1./(1.+z))/(1.+z)
  end function QCDM_angular_diameter_distance

  function QCDM_luminosity_distance (h, Omega_m0, epsilon_s, epsilon_infty, zeta_s, z)
    real(dl) QCDM_luminosity_distance
    real(dl) h, Omega_m0, epsilon_s, epsilon_infty, zeta_s, z
    call QCDM_Initialize(h, Omega_m0, epsilon_s, epsilon_infty, zeta_s)
    QCDM_luminosity_distance = Comoving_Distance_of_a(1./(1.+z))*(1.+z)
  end function QCDM_LUMINOSITY_DISTANCE


End Module QCDM

