/*
   Omega_D - averaged Omega parameters - RZA2 arXiv:1303.6193v2

   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 Omega_D.c */

#include <stdio.h>
#include <sys/types.h>
#include "config.h"
#include <math.h>

#include <gsl/gsl_rng.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_statistics.h>

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

#include "lib/inhomog.h"

#define DEBUG 1

/* #undef DEBUG */

/* Correct the P(k) normalisation */
int correct_Pk_normalisation(/* INPUTS: */
                             struct rza_integrand_params_s *rza_integrand_params, /* +OUT */
                             struct background_cosm_params_s *background_cosm_params, /* +OUT */
                             int want_verbose,
                             long   n_calls_invariants /* for invariant I integration */
                             ){

#define EIGHT_MPC_H100 8.0

  struct rza_integrand_params_s rza_integrand_params_sigma8;

  double sigma8sq_calc, sigma8sq_err_calc;

  /* if not yet normalised */
  if(0 == background_cosm_params->correct_Pk_norm_known){

    /* calculate sigma_8^2 using the first invariant, i.e. same
       as for FLRW linear theory */


    /* Copy the domain shape and power spectrum but
       override the domain size and declare the Ist invariant
       to be unknown. */
    /* TODO: do any rza_integrand_params really need to be inherited? */
    rza_integrand_params_sigma8 = *rza_integrand_params;
    rza_integrand_params_sigma8.w_type = 1; /* spherical */
    rza_integrand_params_sigma8.sigma_sq_inv_triple.I_known = 0;

    rza_integrand_params_sigma8.R_domain = EIGHT_MPC_H100/
      (background_cosm_params->H_0/100.0) *
      background_cosm_params->inhomog_a_scale_factor_initial
      / background_cosm_params->inhomog_a_scale_factor_now;

    rza_integrand_params_sigma8.background_cosm_params =
      *background_cosm_params;

    /* integrate \sigma^2 at 8 Mpc/h0eff at the initial
       time */
    sigma_sq_invariant_I(rza_integrand_params_sigma8,
                         n_calls_invariants,
                         want_verbose,
                         &sigma8sq_calc,
                         &sigma8sq_err_calc
                         );

    background_cosm_params->correct_Pk_norm =
      background_cosm_params->inhomog_a_scale_factor_initial /
      background_cosm_params->inhomog_a_scale_factor_now *
      background_cosm_params->sigma8 /
      sqrt(fmax(sigma8sq_calc,TOL_LENGTH_SQUARED));
    background_cosm_params->correct_Pk_norm *=
      background_cosm_params->correct_Pk_norm;

    background_cosm_params->correct_Pk_norm_known = 1;

  };
  return 0;
}


/*! \brief Calculates \f$ \Omega's \f$ and other important cosmological
 * parameters.
 *
 * Calls \ref kinematical_backreaction, \ref curvature_backreaction and
 * \ref scale_factor_D functions before main calculation to get values
 * of \f$ ^{RZA}Q_D \f$, \f$ ^{RZA}R_D \f$ and the scale factor \f$ a_D
 * \f$. If necessary, calculates the \f$ <I_i> \f$ invariant using the
 * \ref sigma_sq_invariant_I function.
 *
 * Initial time depends on the chosen model. The choice is controlled by
 * the \a t_EdS and \a t_flatFLRW parameters (for a more detailed
 * description see biscale_partition.c). Checks with the \a t_background
 * values to prevent them from being smaller than \a t_initial.
 *
 * To evaluate \f$ \Omega_D \f$, uses (18) and (19) of Buchert et al.
 * 2013; \latexonly \href{https://arxiv.org/abs/1303.6193}{arXiv:
 * 1303.6193} \endlatexonly ).
 *
 *
 * \param [in,out] rza_integrand_params_s pointer to the structure
 * containing parameters necessary for the ODE integration
 * \param [in] background_cosm_params pointer to the
 * background_cosm_params_s containing relevant cosmological parameters
 * \param [in] t_background pointer to the matrix of time values
 * \param [in] n_t_background size of the \a t_background_in matrix
 * \param [in] n_calls_invariants for invariant I integration
 * \param [in] want_planar control parameter; defined in
 * biscale_partition.c
 * \param [in] want_spherical control parameter; defined in
 * biscale_partition.c
 * \param [in] want_verbose control parameter; defined in
 * biscale_partition.c
 * \param [out] rza_Q_D pointer to the RZA kinematical backreaction
 * parameter
 * \param [out] rza_R_D pointer to the RZA Ricci scalar curvature
 * parameter
 * \param [out] a_D pointer to the scale factor \f$ a_D \f$
 * \param [out] dot_a_D pointer to the first derivative of the scale
 * factor \f$ \dot{a_D} \f$
 * \param [out] unphysical pointer to the control parameter for checking
 * if values are physically reasonable
 * \param [out] H_D pointer to the Hubble constant
 * \param [out] Omm_D pointer to the total \f$ \Omega_m \f$
 * \param [out] OmQ_D pointer to \f$ \Omega_Q \f$ (backreaction term)
 * \param [out] OmLam_D pointer to \f$ \Omega_{\Lambda} \f$
 * \param [out] OmR_D pointer to the \f$ \Omega_R \f$ (curvature term)
 */
int Omega_D(/* INPUTS: */
            struct rza_integrand_params_s *rza_integrand_params, /* +OUT */
            struct background_cosm_params_s background_cosm_params,
            double *t_background_in,
            int  n_t_background_in,
            double n_sigma[3],
            long   n_calls_invariants,  /* for invariant I integration */
            int want_planar, /* cf RZA2 V.A */
            int want_spherical, /* cf RZA2 V.B.3 */
            int want_verbose,
            /* OUTPUTS: */
            double *rza_Q_D,
            double *rza_R_D,
            double *a_D,  /* all of size n_t_background_in */
            double *dot_a_D,
            int *unphysical,
            double *H_D,
            double *Omm_D,
            double *OmQ_D,
            double *OmLam_D,
            double *OmR_D
            ){

  /*
    double *rza_Q_D;
    double *rza_R_D;
    int *unphysical;
  */
  double t_initial, a_initial, a_dot_initial;
  double a_FLRW, a_dot_FLRW;
  /* double t_0;  */   /* locally calculate initial and final times */
  int i_t;

  /* Aconst_bis is 2/3 of Aconst in scale_factor_d.c, i.e.
     Aconst_bis = H_H^2(t_i) a_i^3 Omega_m(t_i) (1 - inv_I) */
  double Aconst_bis = 0.0;
  double inv_I, inv_I_err;

  /*   int want_spherical=0; */

  /* prepare for calculating Q_D estimate */
  /*
  rza_Q_D = malloc((size_t)n_t_background_in * sizeof(double));
  rza_R_D = malloc((size_t)n_t_background_in * sizeof(double));
  unphysical = malloc((size_t)n_t_background_in*sizeof(int));
  */
  /* initial time; check that input times are later */
  if(1 == background_cosm_params.EdS){
    t_initial = t_EdS(&background_cosm_params,
                      background_cosm_params.inhomog_a_scale_factor_initial,
                      want_verbose);
  }else if(1 == background_cosm_params.flatFLRW){
    t_initial = t_flatFLRW(&background_cosm_params,
                           background_cosm_params.inhomog_a_scale_factor_initial,
                            want_verbose);
  }else{
      printf("No other options for background_cosm_params so far in program.\n");
      return 1;
  };

  if( gsl_stats_min(t_background_in, 1, (size_t)n_t_background_in)
      < t_initial ){
    printf("scale_factor_D ERROR: smallest t_background_in = %g = too small < %g\n",
           gsl_stats_min(t_background_in, 1, (size_t)n_t_background_in),
           t_initial);
    exit(1);
  };

  /* initial, final times, normally should be first and last */
  /* t_0 = gsl_stats_max(t_background_in, 1, (size_t)n_t_background_in); */ /* not used */



  /* calculate Q_D */
  kinematical_backreaction(rza_integrand_params,
                           background_cosm_params,
                           t_background_in, n_t_background_in,
                           n_sigma,
                           n_calls_invariants,
                           want_planar, /* cf RZA2 V.A */
                           want_spherical,
                           want_verbose,
                           rza_Q_D
                           );

  /* calculate R_D */
  curvature_backreaction(rza_integrand_params,
                         background_cosm_params,
                         t_background_in, n_t_background_in,
                         n_sigma,
                         n_calls_invariants,
                         want_planar, /* cf RZA2 V.A */
                         want_spherical,
                         want_verbose,
                         rza_R_D
                         );

  /* calculate a_D, dot_a_D */
  scale_factor_D(rza_integrand_params,
                 background_cosm_params,
                 t_background_in, n_t_background_in,
                 n_sigma,
                 n_calls_invariants,
                 want_planar, /* cf RZA2 V.A */
                 want_spherical,
                 want_verbose,
                 a_D,
                 dot_a_D,
                 unphysical
                 );

  /* TODO: fix programming risk - the programmer might set
     rza_integrand_params.background_cosm_params
     at a higher level and not expect it to be modified here */
  rza_integrand_params->background_cosm_params =
    background_cosm_params; /* needed by the invariant integrators for P(k)
                             */


  if(rza_integrand_params->precalculated_invariants.enabled){
    inv_I = rza_integrand_params->precalculated_invariants.inv_I;
  }else /* (re-)calculate inv_I if needed */ {
    if(rza_integrand_params->sigma_sq_inv_triple.I_known){
      /* set the scale of inv_I */
      inv_I = n_sigma[0] *
        rza_integrand_params->sigma_sq_inv_triple.sqrt_E_sigma_sq_I;
    }else{
      sigma_sq_invariant_I( *rza_integrand_params,
                            n_calls_invariants,
                            want_verbose,
                            &inv_I, &inv_I_err );
      /* set the scale of inv_I */
      rza_integrand_params->sigma_sq_inv_triple.sqrt_E_sigma_sq_I =
        sqrt(inv_I);
      rza_integrand_params->sigma_sq_inv_triple.I_known = 1;
      inv_I = n_sigma[0] *
        rza_integrand_params->sigma_sq_inv_triple.sqrt_E_sigma_sq_I;
    };
  };


  if(1==background_cosm_params.EdS){
    a_initial =
      a_EdS(&background_cosm_params, t_initial, want_verbose);
    a_dot_initial =
      a_dot_EdS(&background_cosm_params, t_initial, want_verbose);
    Aconst_bis =
      a_dot_initial*a_dot_initial /
      (a_initial*a_initial) *
      background_cosm_params.inhomog_a_d_scale_factor_initial *
      background_cosm_params.inhomog_a_d_scale_factor_initial *
      background_cosm_params.inhomog_a_d_scale_factor_initial *
      (1.0 - inv_I);
  };

  /* Use (18) and (19) of RZA2 to evaluate Omega^D's */
  for(i_t=0; i_t<n_t_background_in; i_t++){
    if(!unphysical[i_t]){
      H_D[i_t] = dot_a_D[i_t] / a_D[i_t];
      OmQ_D[i_t] = -rza_Q_D[i_t] /(6.0* H_D[i_t]*H_D[i_t]);
      OmR_D[i_t] = -rza_R_D[i_t] /(6.0* H_D[i_t]*H_D[i_t]);
      if(1==background_cosm_params.EdS){
        Omm_D[i_t] = Aconst_bis
          / (fmax(TOL_ADOT_SQ_A_OMD,
                  (dot_a_D[i_t]*dot_a_D[i_t] * a_D[i_t])));
        OmLam_D[i_t] = 0.0;
      }else if(1==background_cosm_params.flatFLRW){
        a_FLRW = a_flatFLRW(&background_cosm_params,
                            t_background_in[i_t],
                            want_verbose);
        a_dot_FLRW = a_dot_flatFLRW(&background_cosm_params,
                                    t_background_in[i_t],
                                    want_verbose);
        /* normalise to a_0 = 1*/
        /* TODO: more general formula valid for non-flat FLRW bg */
        /* this formula: e.g. (1.1) arXiv:1303.4444 and
           Omm_D/Omm_{FLRW}
        */
        Omm_D[i_t] = background_cosm_params.Omm_0 /
          (exp(3.0*log(a_FLRW
                       /background_cosm_params.inhomog_a_scale_factor_now)) *
           background_cosm_params.OmLam_0 +
           background_cosm_params.Omm_0 )
          * (a_dot_FLRW*a_dot_FLRW * a_FLRW)
          / (dot_a_D[i_t]*dot_a_D[i_t] * a_D[i_t]);

        OmLam_D[i_t] = background_cosm_params.OmLam_0 *
          background_cosm_params.H_0 *
          background_cosm_params.H_0 *
          (COSM_H_0_INV_GYR*COSM_H_0_INV_GYR) /
          (H_D[i_t]*H_D[i_t]);
      }; /*  if(1==background_cosm_params.EdS) */
    }else{
      /* TODO: feed up unphysical[i_t] to a higher level */
      H_D[i_t] = -9e9;
      OmQ_D[i_t] = -9e9;
      OmR_D[i_t] = -9e9;
      Omm_D[i_t] = -9e9;
      OmLam_D[i_t] = -9e9;
    }; /*(!unphysical[i_t]) */
  };  /* for(i_t=0; i_t<n_t_background_in; i_t++) */

  /*
      free(unphysical);
      free(rza_R_D);
      free(rza_Q_D);
  */
  return 0;
}