/*
   scale_factor_D - wrapper for Raychaudhuri or Hamiltonian ODE integration of a_D

   Copyright (C) 2013-2015 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 scale_factor_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>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>

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

#include "lib/inhomog.h"

#define DEBUG 1

/* #undef DEBUG */


/*! \brief Calculates the collapse time of the structure.
 *
 * Besides the collapse time calculation, this function returns the
 * state of the structure (collapsed/not collapsed). Calculation happens
 * by skipping through expansion phase to the beginning of the collapse,
 * and then monitoring the value of \a a_D and sign of \a dot_a_D.
 *
 * Depending on the \b DETECT_SUDDEN_COLLAPSE macro value, the function
 * may allow the collapse detection even though \a dot_a_D is still
 * positive. This is switched on by default; to turn this off at the
 * compilation level with gcc, use a compile option:
 * -UDETECT_SUDDEN_COLLAPSE .
 *
 * Prints out a programming error and sets the value of \a i_t_collapse
 * to -1 if something went wrong with the iteration variable \a i_t.
 *
 * \param [in] n_t_outputs number of time steps taken into consideration
 * \param [in] a_D pointer to the scale factor
 * \param [in] dot_a_D pointer to the first derivative of the scale factor
 * \param [out] i_t_collapse time step of the collapse
 * \param [out] collapsed boolean; 1 - collapsed, 0 - not collapsed
 *
 */
int find_collapse_time(/* INPUTS */
                       int  n_t_outputs,
                       double *a_D,
                       double *dot_a_D,
                       /* OUTPUTS */
                       int  *i_t_collapse,
                       int  *collapsed /* boolean */
                       ){

  int i_t;
#ifndef INHOM_TURNAROUND_IS_COLLAPSE
  int i_t_post_expansion; /* first step after expansion */
#endif

  /* go to end of expanding phase */
  for(i_t = 1;
      i_t < n_t_outputs &&
        dot_a_D[i_t] >= 0.0;
      i_t++){
  };

#ifndef INHOM_TURNAROUND_IS_COLLAPSE
  /* i_t_post_expansion should represent the first contracting step */
  i_t_post_expansion = i_t;

  /* continue to collapse */
  for(i_t = i_t_post_expansion;
      i_t < n_t_outputs &&
        /* check that a_D has not dropped to very low positive, zero, negative */
        isnormal(a_D[i_t]) &&
        a_D[i_t] > 1.01 * SCALE_FACTOR_A_D_NEARLY_ZERO && /* inhomog.h */
        /* check that dot_a_D is still negative */
        isnormal(dot_a_D[i_t]) &&
        dot_a_D[i_t] < 0.0 &&
        /* if already second collapsing step, then double-check that
           a_D really *is* decreasing, as claimed by dot_a_D */
        (i_t == i_t_post_expansion || a_D[i_t] < a_D[i_t-1] );
      i_t++){
  };
#endif
  /* The a_D, a_D_dot values after collapse may be numerically
     nonsensical: e.g. both positive, or NaN, or \pm infinity;
     otherwise, they are merely physically nonsensical. So avoid them
     by going back to the last valid i_t value. If a_D_dot was not
     negative by the last valid i_t value, then collapse is not
     considered to have been detected, except if the sudden
     collapse case is enabled with DETECT_SUDDEN_COLLAPSE. (To turn
     this off at the compile stage with gcc, use a compile option
     -UDETECT_SUDDEN_COLLAPSE .)
  */
  i_t--;

#define DETECT_SUDDEN_COLLAPSE 1

  if(0 <= i_t && i_t < n_t_outputs){
    if( dot_a_D[i_t] < 0.0 &&
        isnormal(a_D[i_t]) && /* drop the test if NaN or infinity */
        a_D[i_t] > 1.01 * SCALE_FACTOR_A_D_NEARLY_ZERO){
      *collapsed = 1;
      *i_t_collapse = i_t;
#ifdef DETECT_SUDDEN_COLLAPSE
      /* Sudden collapse: If latest i_t value for which
         both a_D and dot_a_D are reasonable still has dot_a_D > 0,
         then this *may* represent a valid collapse if it occurs
         before the end of the integration
      */
    }else if( i_t+1 < n_t_outputs ){
      *collapsed = 1;
      *i_t_collapse = i_t;
#endif
    }else{
      *collapsed = 0;
      *i_t_collapse = n_t_outputs-1;
    };
  }else{
    *collapsed = 0;
    *i_t_collapse = -1;
    printf("find_collapse_time: Warning: programming error. i_t = %d\n",i_t);
    return (-1);
  };
  return 0;
}

/*! \brief Calculates the scale factor.
 *
 * The calculation can be done by either of two integration methods.
 * The default is the Raychaudhuri integrator; to use the Hamiltonian one, switch
 * on \b INHOM_ENABLE_HAMILTON_INTEGRATOR in the file inhomog.c .
 *
 * The Hamiltonian calculation happens according to the scale_factor_D_Ham
 * function based on Buchert et al. 2013 \latexonly
 * (\href{https://arxiv.org/abs/1303.6193}{arXiv: 1303.6193})
 * \endlatexonly . For a more detailed description go to
 * documentation for the file scale_factor_D_Ham.c .
 *
 * The Raychaudhuri calculation happens according to the
 * scale_factor_D_Ray function based on Buchert et al. 2013 \latexonly
 * (\href{https://arxiv.org/abs/1303.6193}{arXiv: 1303.6193})
 * \endlatexonly . For a more detailed description go to
 * documentation for the file scale_factor_D_Ray.c .
 *
 * \param [in] rza_integrand_params_s structure containing parameters
 * necessary for the ODE integration
 * \param [in] background_cosm_params_s structure containing all
 * relevant cosmological parameters (which are defined outside this
 * file)
 * \param [in] t_background pointer to the time values matrix
 * \param [in] n_t_background number of entries in the \a t_background
 * matrix
 * \param [in] n_sigma[3] an array with maximal numbers of standard
 * deviations for each invariant
 * \param [in] n_calls_invariants integration parameter for functions
 * sigma_sq_invariant_I and sigma_sq_invariant_II
 * \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] a_D pointer to the scale factor
 * \param [out] dot_a_D pointer to the first derivative of a_D
 * \param [out] unphysical control parameter for checking if values are
 * physically reasonable
 */
int scale_factor_D(/* INPUTS: */
                   struct rza_integrand_params_s *rza_integrand_params,
                   struct background_cosm_params_s background_cosm_params,
                   double *t_background,
                   int  n_t_background,
                   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 *a_D,
                   double *dot_a_D,
                   int *unphysical
                   ){


#ifdef INHOM_ENABLE_HAMILTON_INTEGRATOR
    /* Hamilton version */
    scale_factor_D_Ham(/* INPUTS: */
                       rza_integrand_params,
                       background_cosm_params,
                       t_background,
                       n_t_background,
                       n_sigma,
                       n_calls_invariants,  /* for invariant I integration */
                       want_planar, /* cf RZA2 V.A */
                       want_spherical, /* cf RZA2 V.B.3 */
                       want_verbose,
                       /* OUTPUTS: */
                       a_D,
                       dot_a_D,
                       unphysical
                       );
#else
    /* Raychaudhuri version */
    scale_factor_D_Ray(/* INPUTS: */
                       rza_integrand_params,
                       background_cosm_params,
                       t_background,
                       n_t_background,
                       n_sigma,
                       n_calls_invariants,  /* for invariant I integration */
                       want_planar, /* cf RZA2 V.A */
                       want_spherical, /* cf RZA2 V.B.3 */
                       want_verbose,
                       /* OUTPUTS: */
                       a_D,
                       dot_a_D,
                       unphysical
                       );
#endif /* ifdef INHOM_ENABLE_HAMILTON_INTEGRATOR */

  return 0;
}