/*
   FLRW background growth functions

   Copyright (C) 2013 Boud Roukema, Jan Ostrowski

   This program is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 2, or (at your option)
   any later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software Foundation,
   Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.

   See also http://www.gnu.org/licenses/gpl.html

*/

/*! \file growth_function.c

 * Calculates the growth function and its derivatives in FLRW
 * background.
 *
 * All values of the growth functions are scaled to match the \f$ a_0
 * = 1 \f$ (scale factor at the present time).
 */

#include <stdio.h>
#include <sys/types.h>
#include "config.h"
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_deriv.h>

/* for malloc_usable_size if available */
#ifdef __GNUC__
#include <malloc.h>
#endif

#include "lib/inhomog.h"

#define KASAI_GROWTH 1

/* *
   FLRW-background growth functions.

   All q_growth values are scaled to match
   a_0 = 1 for the scale factor. This should
   not affect external values of the scale
   factor.
*/


/* * growth function for GSL differentiation */

double growth_FLRW_func(
                        double t_background, void * params)
{

  struct background_cosm_params_s background_cosm_params =
    *(struct background_cosm_params_s *) params;
  int want_verbose = 0;

  return growth_FLRW(&background_cosm_params,
                     t_background,
                     want_verbose);
}

/* * growth function for direct use */
/*! \brief Calculates the growth function for the FLRW background.
 *
 * If the model chosen in the background_cosm_params_s structure is the
 * EdS model, calculates the growth function through the \ref a_EdS
 * function (for a detailed description, see FLRW_background.c).
 *
 * If the chosen model is a flat FLRW model (also chosen through the
 * background_cosm_params_s structure), calculates the growth function
 * according to Bildhauer et al. 1992. \latexonly
 * (\href{http://adsabs.harvard.edu/abs/1992A%26A...263...23B}
 * {1992A\&A...263...23B}). \endlatexonly
 *
 * If the choice of the model is not specified or if a flat FLRW model
 * is chosen with non-zero curvature, prints out an error message (only
 * if \a want_verbose is set to 1).
 *
 * If no error occurs, returns \a q_growth.
 *
 * \param [in] background_cosm_params pointer to the
 * background_cosm_params_s containing relevant cosmological parameters
 * \param [in] t_background time values matrix
 * \param [in] want_verbose control parameter; defined and explained in
 * biscale_partition.c
 */
double growth_FLRW(/* INPUTS: */
                  struct background_cosm_params_s * background_cosm_params,
                  double t_background,
                  int want_verbose
                  /* OUTPUTS: */
                  ){
  double a_scale_factor;
  double q_growth;

#ifdef KASAI_GROWTH
  /* coefficients in Eq (5) & (10)--(16) Kasai 2010 */
  const double b1 = 619226202351.0/527102715964.0; /* 1.175 */
  const double b2 = 12478731282519.0/40730664415400.0; /* 0.3064 */
  const double b3 = 232758215919527.0/43467765064114880.0; /* 0.005355 */
  const double c1 = 88964947071.0/47918428724.0;  /* 1.875 */
  const double c2 = 2445658735707.0/2395921436200.0; /* 1.021 */
  const double c3 = 17010766061223.0/111170754639680.0; /* 0.1530 */
  /* Eq (6) Kasai 2010 */
  double x_eq5_K;
  double x_eq5_K2;
  double x_eq5_K3;
#else
  double OmLam0_a_cubed;
  double x_eq30_BildhauerBuchert1992;
  const double two_thirds = 2.0/3.0;
  const double five_sixths = 5.0/6.0;
  /* Bildhauer & Buchert 1992 Eq (31) normalisation:
     5/6 beta(5/6,2/3) \equiv
     2 / sqrt(pi) * gamma(11/6) * gamma(2/3) */
  const double BildBuchEq31_norm = 1.43728308846099;
#endif

  if(1==background_cosm_params->EdS){
    q_growth = a_EdS( background_cosm_params,
                      t_background,
                      want_verbose ) /
      background_cosm_params->inhomog_a_scale_factor_now;

    return q_growth;

  }else if( (!(background_cosm_params->flatFLRW) || background_cosm_params->Omm_0 > 1.0) && want_verbose){
    printf("t_flatFLRW ERROR: called for invalid bg model.\n");
    exit(1);
  };

  /* internally to this routine must be scaled to a_0 = 1 */
  a_scale_factor = a_flatFLRW( background_cosm_params,
                               t_background,
                               want_verbose ) /
    background_cosm_params->inhomog_a_scale_factor_now;

#ifdef KASAI_GROWTH
  /* Eq (5) Kasai arXiv:1012.2671 */
  x_eq5_K = (1-background_cosm_params->Omm_0)/
    background_cosm_params->Omm_0 *
    a_scale_factor * a_scale_factor * a_scale_factor;
  x_eq5_K2 = x_eq5_K * x_eq5_K;
  x_eq5_K3 = x_eq5_K2 * x_eq5_K;

  q_growth = a_scale_factor * sqrt(1.0 + x_eq5_K) *
    (1.0 + b1*x_eq5_K + b2*x_eq5_K2 + b3*x_eq5_K3)/
    (1.0 + c1*x_eq5_K + c2*x_eq5_K2 + c3*x_eq5_K3);
#else
  OmLam0_a_cubed = background_cosm_params->OmLam_0 *
    exp(3.0 * log(a_scale_factor));
  x_eq30_BildhauerBuchert1992 = OmLam0_a_cubed /
    (background_cosm_params->Omm_0 + OmLam0_a_cubed);

  /* Eq (28) Bildhauer & Buchert 1992 */
  q_growth = BildBuchEq31_norm *
    gsl_sf_beta_inc (five_sixths, two_thirds,
                     x_eq30_BildhauerBuchert1992) *
    exp(log( background_cosm_params->Omm_0 /
             background_cosm_params->OmLam_0 )/3.0) *
    sqrt( 1.0 +
          background_cosm_params->Omm_0 /
          OmLam0_a_cubed );
#endif

  /* DEBUG ONLY */
  /*   q_growth = a_scale_factor; */


  return q_growth;
}


/* * growth function for GSL differentiation */
double dot_growth_FLRW_func(
                        double t_background, void * params)
{

  struct background_cosm_params_s background_cosm_params =
    *(struct background_cosm_params_s *) params;
  int want_verbose = 0;

  return dot_growth_FLRW(&background_cosm_params,
                         t_background,
                         want_verbose);
}


/* First derivative of flat FLRW model growth function with respect to
   the scale factor (not time); see notes on Kasai 2010 below.
   Not normally intended for external use.
*/

static double dD_da_func(double aa, double omr);

static double dD_da_func(double aa, double omr){
  double dD_da;
  double a3omr; /* aa^3 * omr */

  /* maxima - Eq (5) Kasai 2010 arXiv:1012.2671 :
     D : a * sqrt(1+x) * (1 + 1.175*x + 0.3064*x^2 + 0.005355*x^3)/ (1 + 1.857*x + 1.021*x^2 + 0.1530*x^3), x: omr*a^3;
     diff(D, a, 1), factor;
     string(%);
  */

  a3omr = pow(aa,3)*omr;
  dD_da =
    (819315.0e0 *pow(a3omr,7) +
     2980287.0e0 *pow(a3omr,6) +
     66204367.0e0 *pow(a3omr,5) +
     432092178.0e0 *pow(a3omr,4) +
     1276790680.0e0 *pow(a3omr,3) +
     1902310000.0e0 *pow(a3omr,2) +
     1394400000.0e0* a3omr +
     400000000.0e0)/
    (400.0e0*sqrt( a3omr + 1.0e0)*
     pow( (153.0e0 *pow(a3omr,3) +
           1021.0e0 *pow(a3omr,2) +
           1857.0e0 *a3omr +
           1000.0e0), 2) );
  return dD_da;
}



/* * growth function time derivative */
/*! \brief Calculates first time derivative of the FLRW growth factor.
 *
 * The calculation depends on the model used. For EdS, uses \ref
 * a_dot_EdS function (for a more detailed description, go to
 * FLRW_background.c).
 *
 * If the choice of the model is not specified or if a flat FLRW model
 * is chosen with non-zero curvature, prints out an error message (only
 * if \a want_verbose is set to 1).
 *
 * For flat FLRW model, uses predefined GSL library routines to
 * calculate the first derivative. The routine used depends on the
 * relative step size (in comparison with \a t_background values). If
 * \a t_background value is bigger than \f$ 5h \f$, where \f$ h \f$ is
 * the step value, then the derivation happens according to the \a
 * gsl_deriv_central routine; if not, the routine used is \a
 * gsl_deriv_forward.
 *
 * If no error occurs, returns \a dq_dt.
 *
 * \param [in] background_cosm_params pointer to the
 * background_cosm_params_s containing relevant cosmological parameters
 * \param [in] t_background time values matrix
 * \param [in] want_verbose control parameter; defined and explained in
 * biscale_partition.c
 */
double dot_growth_FLRW(/* INPUTS: */
                       struct background_cosm_params_s * background_cosm_params,
                       double t_background,
                       int want_verbose
                       /* OUTPUTS: */
                       ){
#define DELTA_T_DERIV_TOL 1e-5

  double dq_dt;

#ifdef KASAI_GROWTH
  double aa, a_dot;
  double omr; /* (1-Omega_{m0})/Omega_{m0} - eq 6 Kasai 2010 */
#else
  gsl_function F_gsl;
  double h_step;
  double abserr;
#endif

  if(1==background_cosm_params->EdS){
    dq_dt = a_dot_EdS( background_cosm_params,
                          t_background,
                          want_verbose )  /
      background_cosm_params->inhomog_a_scale_factor_now;
    return dq_dt;

  }else if( (!(background_cosm_params->flatFLRW) || background_cosm_params->Omm_0 > 1.0) && want_verbose){
    printf("t_flatFLRW ERROR: called for invalid bg model.\n");
    exit(1);
  };


#ifdef KASAI_GROWTH

  /* internally to this routine must be scaled to a_0 = 1 */
  aa = a_flatFLRW( background_cosm_params,
                               t_background,
                               want_verbose ) /
    background_cosm_params->inhomog_a_scale_factor_now;
  a_dot = a_dot_flatFLRW( background_cosm_params,
                          t_background,
                          want_verbose ) /
    background_cosm_params->inhomog_a_scale_factor_now;
  omr = (1-background_cosm_params->Omm_0)/
    background_cosm_params->Omm_0;

  dq_dt = dD_da_func(aa, omr) * a_dot;

#else
  F_gsl.function = &growth_FLRW_func;
  F_gsl.params = background_cosm_params;

  /* use DELTA_T_DERIV_TOL to make a relative step size guess */
  h_step = DELTA_T_DERIV_TOL * t_background;

  if(t_background > 5*h_step){
    gsl_deriv_central(&F_gsl, t_background,
                      h_step,
                      &dq_dt,&abserr);
  }else{
    gsl_deriv_forward(&F_gsl, t_background,
                      h_step,
                      &dq_dt,&abserr);
  };
#endif

  return dq_dt;
}


/* * growth function second derivative wrt time */
/*! \brief Calculates second time derivative of the FLRW growth factor.
 *
 * The calculation depends on the model used. For EdS, uses \ref
 * a_ddot_EdS function (for a more detailed description, go to
 * FLRW_background.c).
 *
 * If the choice of the model is not specified or if a flat FLRW model
 * is chosen with non-zero curvature, prints out an error message (only
 * if \a want_verbose is set to 1).
 *
 * For flat FLRW model, uses predefined GSL library routines to
 * calculate the first derivative. The routine used depends on the
 * relative step size (in comparison with \a t_background values). If
 * \a t_background value is bigger than \f$ 5h \f$, where \f$ h \f$ is
 * the step value, then the derivation happens according to the \a
 * gsl_deriv_central routine; if not, the routine used is \a
 * gsl_deriv_forward.
 *
 * If no error occurs, returns \a d2q_dt2.
 *
 * \param [in] background_cosm_params pointer to the
 * background_cosm_params_s containing relevant cosmological parameters
 * \param [in] t_background time values matrix
 * \param [in] want_verbose control parameter; defined and explained in
 * biscale_partition.c
 */
double ddot_growth_FLRW(/* INPUTS: */
                        struct background_cosm_params_s * background_cosm_params,
                        double t_background,
                        int want_verbose
                        /* OUTPUTS: */
                        ){
  /* non-integer factor greater than DELTA_T_DERIV_TOL */
#define DELTA_T_DDOT_DERIV_TOL 2.3e-5

  double d2q_dt2;

#ifdef KASAI_GROWTH
  double aa, a_dot, a_ddot;
  double omr; /* (1-Omega_{m0})/Omega_{m0} - eq 6 Kasai 2010 */
  double d2D_da2;
  double a3omr; /* aa^3 * omr */
#else
  gsl_function F_gsl;
  double h_step;
  double abserr;
#endif

  if(1==background_cosm_params->EdS){
    d2q_dt2 = a_ddot_EdS( background_cosm_params,
                           t_background,
                           want_verbose ) /
      background_cosm_params->inhomog_a_scale_factor_now;
    return d2q_dt2;

  }else if( (!(background_cosm_params->flatFLRW) || background_cosm_params->Omm_0 > 1.0) && want_verbose){
    printf("t_flatFLRW ERROR: called for invalid bg model.\n");
    exit(1);
  };

#ifdef KASAI_GROWTH
  /* internally to this routine must be scaled to a_0 = 1 */
  aa = a_flatFLRW( background_cosm_params,
                               t_background,
                               want_verbose ) /
    background_cosm_params->inhomog_a_scale_factor_now;
  a_dot = a_dot_flatFLRW( background_cosm_params,
                          t_background,
                          want_verbose ) /
    background_cosm_params->inhomog_a_scale_factor_now;
  a_ddot = a_ddot_flatFLRW( background_cosm_params,
                            t_background,
                            want_verbose ) /
    background_cosm_params->inhomog_a_scale_factor_now;

  omr = (1-background_cosm_params->Omm_0)/
    background_cosm_params->Omm_0;

  /* maxima - Eq (5) Kasai 2010:
     D : a * sqrt(1+x) * (1 + 1.175*x + 0.3064*x^2 + 0.005355*x^3)/ (1 + 1.857*x + 1.021*x^2 + 0.1530*x^3), x: omr*a^3;
     diff(D, a, 2), factor;
     string(%);
  */

  /*
  (3*a^2*omr*
    (125355195*a^30*omr^10 +  3977329554*a^27*omr^9 -  2546850087*a^24*omr^8 -  206436376394*a^21*omr^7 -  1234623769273*a^18*omr^6 -  3925947239328*a^15*omr^5 -  7726431851616*a^12*omr^4 -  9613741100480*a^9*omr^3 -  7352622720000*a^6*omr^2 -  3156361600000*a^3*omr -  582400000000))
    /(800*(a^3*omr +  1)^(3/2)*
    (153*a^9*omr^3 +  1021* a^6*omr^2 +  1857*a^3*omr +  1000)^3)
  */


  a3omr = pow(aa,3)*omr;
  d2D_da2 = (3.0e0 *pow(aa,2)*omr*
             (125355195.0e0 * pow(a3omr,10) +
              3977329554.0e0 *pow(a3omr,9) -
              2546850087.0e0 *pow(a3omr,8) -
              206436376394.0e0 *pow(a3omr,7) -
              1234623769273.0e0 *pow(a3omr,6) -
              3925947239328.0e0 *pow(a3omr,5) -
              7726431851616.0e0 * pow(a3omr,4) -
              9613741100480.0e0 * pow(a3omr,3) -
              7352622720000.0e0 * pow(a3omr,2) -
              3156361600000.0e0 * a3omr-
              582400000000.0e0))/
    (800.0e0*
     pow((a3omr + 1.0e0),1.5)*
     pow((153.0e0 * pow(a3omr,3) +
          1021.0e0 *pow(a3omr,2) +
          1857.0e0 *a3omr +
          1000.0e0),3)
     );

  /* \ddot{y} = y'' \dot{x}^2 + y' \ddot{x} where ' = d/dx */
  d2q_dt2 = d2D_da2 *a_dot *a_dot +
    dD_da_func(aa, omr) * a_ddot;

#else

  F_gsl.function = &dot_growth_FLRW_func;
  F_gsl.params = background_cosm_params;

    /* use DELTA_T_DDOT_DERIV_TOL as a relative step size guess */
  h_step = DELTA_T_DDOT_DERIV_TOL * t_background;

  if(t_background > 5*h_step){
    gsl_deriv_central(&F_gsl, t_background,
                      h_step,
                      &d2q_dt2,&abserr);
  }else{
    gsl_deriv_forward(&F_gsl, t_background,
                      h_step,
                      &d2q_dt2,&abserr);
  };
#endif

  return d2q_dt2;
}