text:
/* integration/qk21.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
 *
 * 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 3 of the License, 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <gsl/gsl_integration.h>

/* Gauss quadrature weights and kronrod quadrature abscissae and
   weights as evaluated with 80 decimal digit arithmetic by
   L. W. Fullerton, Bell Labs, Nov. 1981. */

static const double xgk[11] =   /* abscissae of the 21-point kronrod rule */
{
  0.995657163025808080735527280689003,
  0.973906528517171720077964012084452,
  0.930157491355708226001207180059508,
  0.865063366688984510732096688423493,
  0.780817726586416897063717578345042,
  0.679409568299024406234327365114874,
  0.562757134668604683339000099272694,
  0.433395394129247190799265943165784,
  0.294392862701460198131126603103866,
  0.148874338981631210884826001129720,
  0.000000000000000000000000000000000
};

/* xgk[1], xgk[3], ... abscissae of the 10-point gauss rule.
   xgk[0], xgk[2], ... abscissae to optimally extend the 10-point gauss rule */

static const double wg[5] =     /* weights of the 10-point gauss rule */
{
  0.066671344308688137593568809893332,
  0.149451349150580593145776339657697,
  0.219086362515982043995534934228163,
  0.269266719309996355091226921569469,
  0.295524224714752870173892994651338
};

static const double wgk[11] =   /* weights of the 21-point kronrod rule */
{
  0.011694638867371874278064396062192,
  0.032558162307964727478818972459390,
  0.054755896574351996031381300244580,
  0.075039674810919952767043140916190,
  0.093125454583697605535065465083366,
  0.109387158802297641899210590325805,
  0.123491976262065851077958109831074,
  0.134709217311473325928054001771707,
  0.142775938577060080797094273138717,
  0.147739104901338491374841515972068,
  0.149445554002916905664936468389821
};

void
gsl_integration_qk21 (const gsl_function * f, double a, double b,
                      double *result, double *abserr,
                      double *resabs, double *resasc)
{
  double fv1[11], fv2[11];
  gsl_integration_qk (11, xgk, wg, wgk, fv1, fv2, f, a, b, result, abserr, resabs, resasc);
}


text:
#ifndef __GSL_VECTOR_H__
#define __GSL_VECTOR_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <gsl/gsl_vector_complex_long_double.h>
#include <gsl/gsl_vector_complex_double.h>
#include <gsl/gsl_vector_complex_float.h>

#include <gsl/gsl_vector_long_double.h>
#include <gsl/gsl_vector_double.h>
#include <gsl/gsl_vector_float.h>

#include <gsl/gsl_vector_ulong.h>
#include <gsl/gsl_vector_long.h>

#include <gsl/gsl_vector_uint.h>
#include <gsl/gsl_vector_int.h>

#include <gsl/gsl_vector_ushort.h>
#include <gsl/gsl_vector_short.h>

#include <gsl/gsl_vector_uchar.h>
#include <gsl/gsl_vector_char.h>

#endif /* __GSL_VECTOR_H__ */

ext: h

max_stars_repo_name: berkus/music-cs

max_issues_repo_name: berkus/music-cs

max_forks_repo_name: berkus/music-cs


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv2.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_roots.h>

#include "ccl.h"

/* --------- ROUTINE: h_over_h0 ---------
INPUT: scale factor, cosmology
TASK: Compute E(a)=H(a)/H0
*/
static double h_over_h0(double a, ccl_cosmology * cosmo, int *status)
{
  // Check if massive neutrinos are present - if not, we don't need to
  // compute their contribution
  double Om_mass_nu;
  if ((cosmo->params.N_nu_mass)>1e-12) {
    Om_mass_nu = ccl_Omeganuh2(
      a, cosmo->params.N_nu_mass, cosmo->params.m_nu, cosmo->params.T_CMB,
      status) / (cosmo->params.h) / (cosmo->params.h);
  }
  else {
    Om_mass_nu = 0;
  }

  /* Calculate h^2 using the formula (eqn 2 in the CCL paper):
    E(a)^2 = Omega_m a^-3 +
             Omega_l a^(-3*(1+w0+wa)) exp(3*wa*(a-1)) +
             Omega_k a^-2 +
             (Omega_g + Omega_nu_rel) a^-4 +
             Om_mass_nu
  */
  return sqrt(
    (cosmo->params.Omega_c + cosmo->params.Omega_b +
     cosmo->params.Omega_l *
       pow(a,-3*(cosmo->params.w0+cosmo->params.wa)) *
       exp(3*cosmo->params.wa*(a-1)) +
     cosmo->params.Omega_k * a +
     (cosmo->params.Omega_g + cosmo->params.Omega_nu_rel) / a +
     Om_mass_nu * a*a*a) / (a*a*a));
}

/* --------- ROUTINE: ccl_omega_x ---------
INPUT: cosmology object, scale factor, species label
TASK: Compute the density relative to critical, Omega(a) for a given species.
Possible values for "label":
ccl_species_crit_label <- critical (physical)
ccl_species_m_label <- matter
ccl_species_l_label <- DE
ccl_species_g_label <- radiation
ccl_species_k_label <- curvature
ccl_species_ur_label <- massless neutrinos
ccl_species_nu_label <- massive neutrinos
*/
double ccl_omega_x(ccl_cosmology * cosmo, double a, ccl_species_x_label label, int *status)
{
  // If massive neutrinos are present, compute the phase-space integral and
  // get OmegaNuh2. If not, set OmegaNuh2 to zero.
  double OmNuh2;
  if ((cosmo->params.N_nu_mass) > 0.0001) {
    // Call the massive neutrino density function just once at this redshift.
    OmNuh2 = ccl_Omeganuh2(a, cosmo->params.N_nu_mass, cosmo->params.m_nu,
                           cosmo->params.T_CMB, status);
  }
  else {
    OmNuh2 = 0.;
  }

  double hnorm = h_over_h0(a, cosmo, status);

  switch(label) {
    case ccl_species_crit_label :
      return 1.;
    case ccl_species_m_label :
      return (cosmo->params.Omega_c + cosmo->params.Omega_b) / (a*a*a) / hnorm / hnorm +
          OmNuh2 / (cosmo->params.h) / (cosmo->params.h) / hnorm / hnorm;
    case ccl_species_l_label :
      return
        cosmo->params.Omega_l *
        pow(a,-3 * (1 + cosmo->params.w0 + cosmo->params.wa)) *
        exp(3 * cosmo->params.wa * (a-1)) / hnorm / hnorm;
    case ccl_species_g_label :
      return cosmo->params.Omega_g / (a*a*a*a) / hnorm / hnorm;
    case ccl_species_k_label :
      return cosmo->params.Omega_k / (a*a) / hnorm / hnorm;
    case ccl_species_ur_label :
      return cosmo->params.Omega_nu_rel / (a*a*a*a) / hnorm / hnorm;
    case ccl_species_nu_label :
      return OmNuh2 / (cosmo->params.h) / (cosmo->params.h) / hnorm / hnorm;
    default:
      *status = CCL_ERROR_PARAMETERS;
      ccl_cosmology_set_status_message(
        cosmo, "ccl_background.c: ccl_omega_x(): Species %d not supported\n", label);
      return NAN;
  }
}

/* --------- ROUTINE: ccl_rho_x ---------
INPUT: cosmology object, scale factor, species label
TASK: Compute rho_x(a), with x defined by species label.
Possible values for "label":
ccl_species_crit_label <- critical (physical)
ccl_species_m_label <- matter (physical)
ccl_species_l_label <- DE (physical)
ccl_species_g_label <- radiation (physical)
ccl_species_k_label <- curvature (physical)
ccl_species_ur_label <- massless neutrinos (physical)
ccl_species_nu_label <- massive neutrinos (physical)
*/
double ccl_rho_x(ccl_cosmology * cosmo, double a, ccl_species_x_label label, int is_comoving, int *status)
{
  double comfac;
  if (is_comoving) {
     comfac = a*a*a;
  } else {
      comfac = 1.0;
  }
  double hnorm = h_over_h0(a, cosmo, status);
  double rhocrit =
    ccl_constants.RHO_CRITICAL *
    (cosmo->params.h) *
    (cosmo->params.h) * hnorm * hnorm * comfac;

  return rhocrit * ccl_omega_x(cosmo, a, label, status);
}

/* --------- ROUTINE: ccl_mu_MG ---------
INPUT: cosmology object, scale factor
TASK: Compute mu(a) where mu is one of the the parameterizating functions
of modifications to GR in the quasistatic approximation.
*/

double ccl_mu_MG(ccl_cosmology * cosmo, double a, int *status)
{
    // This function can be extended to include other
    // z-dependences for mu in the future.
    return cosmo->params.mu_0 * ccl_omega_x(cosmo, a, ccl_species_l_label, status) / cosmo->params.Omega_l;
}

/* --------- ROUTINE: ccl_Sig_MG ---------
INPUT: cosmology object, scale factor
TASK: Compute Sigma(a) where Sigma is one of the the parameterizating functions
of modifications to GR in the quasistatic approximation.
*/

double ccl_Sig_MG(ccl_cosmology * cosmo, double a, int *status)
{
    // This function can be extended to include other
    // z-dependences for Sigma in the future.
    return cosmo->params.sigma_0 * ccl_omega_x(cosmo, a, ccl_species_l_label, status) / cosmo->params.Omega_l;
}

// Structure to hold parameters of chi_integrand
typedef struct {
  ccl_cosmology *cosmo;
  int * status;
} chipar;

/* --------- ROUTINE: chi_integrand ---------
INPUT: scale factor
TASK: compute the integrand of the comoving distance
*/
static double chi_integrand(double a, void * params_void)
{
  ccl_cosmology * cosmo = ((chipar *)params_void)->cosmo;
  int *status = ((chipar *)params_void)->status;

  return ccl_constants.CLIGHT_HMPC/(a*a*h_over_h0(a, cosmo, status));
}

/* --------- ROUTINE: growth_ode_system ---------
INPUT: scale factor
TASK: Define the ODE system to be solved in order to compute the growth (of the density)
*/
static int growth_ode_system(double a,const double y[],double dydt[],void *params)
{
  int status = 0;
  ccl_cosmology * cosmo = params;

  double hnorm=h_over_h0(a,cosmo, &status);
  double om=ccl_omega_x(cosmo, a, ccl_species_m_label, &status);

  dydt[1]=1.5*hnorm*a*om*y[0];
  dydt[0]=y[1]/(a*a*a*hnorm);

  return status;
}

/* --------- ROUTINE: growth_ode_system_muSig ---------
INPUT: scale factor
TASK: Define the ODE system to be solved in order to compute the growth (of the density)
* in the case in which we use the mu / Sigma quasistatic parameterisation of modified gravity
*/
static int growth_ode_system_muSig(double a,const double y[],double dydt[],void *params)
{
  int status = 0;
  ccl_cosmology * cosmo = params;

  double hnorm=h_over_h0(a,cosmo, &status);
  double om=ccl_omega_x(cosmo, a, ccl_species_m_label, &status);

  double mu = ccl_mu_MG(cosmo, a, &status);
  dydt[1]=1.5*hnorm*a*om*y[0]*(1. + mu);

  dydt[0]=y[1]/(a*a*a*hnorm);

  return status;
}

/* --------- ROUTINE: df_integrand ---------
INPUT: scale factor, spline object
TASK: Compute integrand from modified growth function
*/
static double df_integrand(double a,void * spline_void)
{
  if(a<=0)
    return 0;
  else {
    gsl_spline *df_a_spline=(gsl_spline *)spline_void;

    return gsl_spline_eval(df_a_spline,a,NULL)/a;
  }
}

/* --------- ROUTINE: growth_factor_and_growth_rate ---------
INPUT: scale factor, cosmology
TASK: compute the growth (D(z)) and the growth rate, logarithmic derivative (f?)
*/

static int growth_factor_and_growth_rate(double a, double *gf, double *fg, ccl_cosmology *cosmo, int *stat)
{
  if(a < cosmo->gsl_params.EPS_SCALEFAC_GROWTH) {
    *gf = a;
    *fg = 1;
    return 0;
  }
  else {
    int gslstatus;
    double y[2];
    double ainit = cosmo->gsl_params.EPS_SCALEFAC_GROWTH;

    // if mu0 == 0, call normal growth_ode_system, otherwise call growth_ode_system_muSig
    if (cosmo->params.mu_0 > 1e-12 || cosmo->params.mu_0 < -1e-12) {
      gsl_odeiv2_system sys = {growth_ode_system_muSig, NULL, 2, cosmo};

      gsl_odeiv2_driver *d = gsl_odeiv2_driver_alloc_y_new(
        &sys, gsl_odeiv2_step_rkck,
        0.1*cosmo->gsl_params.EPS_SCALEFAC_GROWTH, 0, cosmo->gsl_params.ODE_GROWTH_EPSREL);

      if (d == NULL) {
        return CCL_ERROR_MEMORY;
      }

      y[0] = cosmo->gsl_params.EPS_SCALEFAC_GROWTH;
      y[1] = (
        cosmo->gsl_params.EPS_SCALEFAC_GROWTH *
        cosmo->gsl_params.EPS_SCALEFAC_GROWTH *
        cosmo->gsl_params.EPS_SCALEFAC_GROWTH *
        h_over_h0(cosmo->gsl_params.EPS_SCALEFAC_GROWTH, cosmo, stat));

      gslstatus = gsl_odeiv2_driver_apply(d, &ainit, a, y);
      gsl_odeiv2_driver_free(d);

      if(gslstatus != GSL_SUCCESS) {
        ccl_raise_gsl_warning(gslstatus, "ccl_background.c: growth_factor_and_growth_rate():");
        return 0;
      }

      *gf = y[0];
      *fg = y[1]/(a*a*h_over_h0(a, cosmo, stat)*y[0]);
      return 0;
    }
    else {
      gsl_odeiv2_system sys = {growth_ode_system, NULL, 2, cosmo};

      gsl_odeiv2_driver *d = gsl_odeiv2_driver_alloc_y_new(
        &sys, gsl_odeiv2_step_rkck,
        0.1*cosmo->gsl_params.EPS_SCALEFAC_GROWTH, 0, cosmo->gsl_params.ODE_GROWTH_EPSREL);

      if (d == NULL) {
        return CCL_ERROR_MEMORY;
      }

      y[0] = cosmo->gsl_params.EPS_SCALEFAC_GROWTH;
      y[1] = (
        cosmo->gsl_params.EPS_SCALEFAC_GROWTH *
        cosmo->gsl_params.EPS_SCALEFAC_GROWTH *
        cosmo->gsl_params.EPS_SCALEFAC_GROWTH*
        h_over_h0(cosmo->gsl_params.EPS_SCALEFAC_GROWTH, cosmo, stat));

      gslstatus = gsl_odeiv2_driver_apply(d, &ainit, a, y);
      gsl_odeiv2_driver_free(d);

      if(gslstatus != GSL_SUCCESS) {
        ccl_raise_gsl_warning(gslstatus, "ccl_background.c: growth_factor_and_growth_rate():");
        return 0;
      }

      *gf = y[0];
      *fg = y[1]/(a*a*h_over_h0(a, cosmo, stat)*y[0]);

      return 0;
    }
  }
}

/* --------- ROUTINE: compute_chi ---------
INPUT: scale factor, cosmology
OUTPUT: chi -> radial comoving distance
TASK: compute radial comoving distance at a
*/
void compute_chi(double a, ccl_cosmology *cosmo, double * chi, int * stat)
{
  int gslstatus;
  double result;
  chipar p;

  p.cosmo=cosmo;
  p.status=stat;

  gsl_integration_cquad_workspace * workspace = NULL;

  gsl_function F;
  F.function = &chi_integrand;
  F.params = &p;

  workspace = gsl_integration_cquad_workspace_alloc(cosmo->gsl_params.N_ITERATION);

  if (workspace == NULL) {
      *stat = CCL_ERROR_MEMORY;
  } else {
    //TODO: CQUAD is great, but slower than other methods. This could be sped up if it becomes an issue.
    gslstatus=gsl_integration_cquad(
      &F, a, 1.0, 0.0, cosmo->gsl_params.INTEGRATION_DISTANCE_EPSREL, workspace, &result, NULL, NULL);
    *chi=result/cosmo->params.h;

    if (gslstatus != GSL_SUCCESS) {
      ccl_raise_gsl_warning(gslstatus, "ccl_background.c: compute_chi():");
      *stat = CCL_ERROR_COMPUTECHI;
    }
  }

  gsl_integration_cquad_workspace_free(workspace);
}

//Root finding for a(chi)
typedef struct {
  double chi;
  ccl_cosmology *cosmo;
  int * status;
} Fpar;

static double fzero(double a,void *params)
{
  double chi,chia,a_use=a;

  chi=((Fpar *)params)->chi;
  compute_chi(a_use,((Fpar *)params)->cosmo,&chia, ((Fpar *)params)->status);

  return chi-chia;
}

static double dfzero(double a,void *params)
{
  ccl_cosmology *cosmo=((Fpar *)params)->cosmo;
  int *stat = ((Fpar *)params)->status;

  chipar p;
  p.cosmo=cosmo;
  p.status=stat;

  return chi_integrand(a,&p)/cosmo->params.h;
}

static void fdfzero(double a,void *params,double *f,double *df)
{
  *f=fzero(a,params);
  *df=dfzero(a,params);
}

/* --------- ROUTINE: a_of_chi ---------
INPUT: comoving distance chi, cosmology, stat, a_old, gsl_root_fdfsolver
OUTPUT: scale factor
TASK: compute the scale factor that corresponds to a given comoving distance chi
Note: This routine uses a root solver to find an a such that compute_chi(a) = chi.
The root solver uses the derivative of compute_chi (which is chi_integrand) and
the value itself.
*/
static void a_of_chi(double chi, ccl_cosmology *cosmo, int* stat, double *a_old, gsl_root_fdfsolver *s)
{
  if(chi==0) {
    *a_old=1;
  }
  else {
    Fpar p;
    gsl_function_fdf FDF;
    double a_previous,a_current=*a_old;

    p.cosmo=cosmo;
    p.chi=chi;
    p.status=stat;
    FDF.f=&fzero;
    FDF.df=&dfzero;
    FDF.fdf=&fdfzero;
    FDF.params=&p;
    gsl_root_fdfsolver_set(s,&FDF,a_current);

    int iter=0, gslstatus;
    do {
      iter++;
      gslstatus=gsl_root_fdfsolver_iterate(s);
      if(gslstatus!=GSL_SUCCESS) ccl_raise_gsl_warning(gslstatus, "ccl_background.c: a_of_chi():");
      a_previous=a_current;
      a_current=gsl_root_fdfsolver_root(s);
      gslstatus=gsl_root_test_delta(a_current, a_previous, 0, cosmo->gsl_params.ROOT_EPSREL);
    } while(gslstatus==GSL_CONTINUE && iter <= cosmo->gsl_params.ROOT_N_ITERATION);

    *a_old=a_current;

    // Allows us to pass a status to h_over_h0 for the neutrino integral calculation.
    if(gslstatus==GSL_SUCCESS) {
      *stat = *(p.status);
    }
    else {
      ccl_raise_gsl_warning(gslstatus, "ccl_background.c: a_of_chi():");
      *stat = CCL_ERROR_COMPUTECHI;
    }
  }
}

/* ----- ROUTINE: ccl_cosmology_compute_distances ------
INPUT: cosmology
TASK: if not already there, make a table of comoving distances and of E(a)
*/

void ccl_cosmology_compute_distances(ccl_cosmology * cosmo, int *status)
{
  //Do nothing if everything is computed already
  if(cosmo->computed_distances)
    return;

  // Create logarithmically and then linearly-spaced values of the scale factor
  int na = cosmo->spline_params.A_SPLINE_NA+cosmo->spline_params.A_SPLINE_NLOG-1;
  double * a = ccl_linlog_spacing(
    cosmo->spline_params.A_SPLINE_MINLOG, cosmo->spline_params.A_SPLINE_MIN,
    cosmo->spline_params.A_SPLINE_MAX, cosmo->spline_params.A_SPLINE_NLOG,
    cosmo->spline_params.A_SPLINE_NA);
  // Allocate arrays for all three of E(a), chi(a), and a(chi)
  double *E_a = malloc(sizeof(double)*na);
  double *chi_a = malloc(sizeof(double)*na);
  // Allocate E(a) and chi(a) splines
  gsl_spline * E = gsl_spline_alloc(cosmo->spline_params.A_SPLINE_TYPE, na);
  gsl_spline * chi = gsl_spline_alloc(cosmo->spline_params.A_SPLINE_TYPE, na);
  // a(chi) spline allocated below

  //Check for too little memory
  if (a == NULL || E_a == NULL || chi_a == NULL || E == NULL || chi == NULL) {
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(
      cosmo, "ccl_background.c: ccl_cosmology_compute_distances(): ran out of memory\n");
  }

  //Check for messed up scale factor conditions
  if (!*status){
    if ((fabs(a[0]-cosmo->spline_params.A_SPLINE_MINLOG)>1e-5) ||
        (fabs(a[na-1]-cosmo->spline_params.A_SPLINE_MAX)>1e-5) ||
        (a[na-1]>1.0)) {
      *status = CCL_ERROR_LINSPACE;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_background.c: ccl_cosmology_compute_distances(): Error "
        "creating first logarithmic and then linear spacing in a\n");
    }
  }

  // Fill in E(a) - note, this step cannot change the status variable
  if (!*status)
    for (int i=0; i<na; i++)
      E_a[i] = h_over_h0(a[i], cosmo, status);

  // Create a E(a) spline
  if (!*status){
    if (gsl_spline_init(E, a, E_a, na)){
      *status = CCL_ERROR_SPLINE;
      ccl_cosmology_set_status_message(
        cosmo, "ccl_background.c: ccl_cosmology_compute_distances(): Error creating  E(a) spline\n");
    }
  }

  // Compute chi(a)
  if (!*status){
    for (int i=0; i<na; i++)
      compute_chi(a[i], cosmo, &chi_a[i], status);
    if (*status){
      *status = CCL_ERROR_INTEG;
      ccl_cosmology_set_status_message(
        cosmo, "ccl_background.c: ccl_cosmology_compute_distances(): chi(a) integration error \n");
    }
  }

  // Initialize chi(a) spline
  if (!*status){
    if (gsl_spline_init(chi, a, chi_a, na)){//in Mpc
      *status = CCL_ERROR_SPLINE;
      ccl_cosmology_set_status_message(
        cosmo, "ccl_background.c: ccl_cosmology_compute_distances(): Error creating  chi(a) spline\n");
    }
  }

  if (*status){ //If there was an error, free the GSL splines and return
    gsl_spline_free(E); // Note: you are allowed to call gsl_free() on NULL
    gsl_spline_free(chi);
    E = NULL;
    chi = NULL;
  }

  // Set up the boundaries for the a(chi) spline
  double dchi, chi0, chif, a0, af;
  if(!*status){
    dchi=5.;
    chi0=chi_a[na-1];
    chif=chi_a[0];
    a0=a[na-1];
    af=a[0];
  }

  //TODO: The interval in chi (5. Mpc) should be made a macro
  free(a); //Free these, in preparation for making a(chi) splines
  free(E_a);
  free(chi_a);
  a = NULL;
  E_a = NULL;
  chi_a = NULL;
  //Note: you are allowed to call free() on NULL

  ////////////////////////////////////////////////////////////////////////////////////////////////////
  //Below here na (length of some arrays) changes, so this function has to be split at this point.
  ////////////////////////////////////////////////////////////////////////////////////////////////////

  na = (int)((chif-chi0)/dchi);
  dchi  = (chif-chi0)/na; // <=5, since na is an integer
  //Allocate new arrays for a and chi(a)
  chi_a = ccl_linear_spacing(chi0, chif, na);
  a     = malloc(sizeof(double)*na);
  //Allocate space for GSL root finders
  const gsl_root_fdfsolver_type *T=gsl_root_fdfsolver_newton;
  gsl_root_fdfsolver *s = NULL;
  gsl_spline *achi;

  s = gsl_root_fdfsolver_alloc(T);
  achi = gsl_spline_alloc(cosmo->spline_params.A_SPLINE_TYPE, na);

  //Check for too little memory
  if (!*status){
    if (a == NULL || chi_a == NULL || s == NULL || achi == NULL) {
      *status=CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(
        cosmo, "ccl_background.c: ccl_cosmology_compute_distances(): ran out of memory\n");
    }
    else if (fabs(chi_a[0]-chi0) > 1e-5 || fabs(chi_a[na-1]-chif) > 1e-5) { //Check for messed up chi conditions
      *status = CCL_ERROR_LINSPACE;
      ccl_cosmology_set_status_message(
        cosmo, "ccl_background.c: ccl_cosmology_compute_distances(): Error creating linear spacing in chi\n");
    }
  }

  // Calculate a(chi)
  if (!*status){
    a[0]=a0; a[na-1]=af;
    for(int i=1;i<na-1;i++) {
      // we are using the previous value as a guess here to help the root finder
      // as long as we use small steps in a this should be fine
      a_of_chi(chi_a[i],cosmo, status, &a0, s);
      a[i]=a0;
    }
    if(*status) {
      *status = CCL_ERROR_ROOT;
      ccl_cosmology_set_status_message(
        cosmo, "ccl_background.c: ccl_cosmology_compute_distances(): a(chi) root-finding error \n");
    }
  }

  // Initialize the a(chi) spline
  if (!*status){
    if(gsl_spline_init(achi, chi_a, a, na)){
      *status = CCL_ERROR_SPLINE;
      ccl_cosmology_set_status_message(
        cosmo, "ccl_background.c: ccl_cosmology_compute_distances(): Error creating  a(chi) spline\n");
    }
  }

  free(a);
  free(chi_a); //Note: you are allowed to call free() on NULL
  gsl_root_fdfsolver_free(s);
  if (*status){//If there was an error, free the GSL splines and return
    gsl_spline_free(E); //Note: you are allowed to call gsl_free() on NULL
    gsl_spline_free(chi);
    gsl_spline_free(achi);
  }

  if (*status == 0) {
    //If there were no errors, attach the splines to the cosmo struct and end the function.
    cosmo->data.E             = E;
    cosmo->data.chi           = chi;
    cosmo->data.achi          = achi;
    cosmo->computed_distances = true;
  }
}

/* ----- ROUTINE: ccl_cosmology_distances_from_input ------
INPUT: cosmology, scale factor array, comoving distance chi computed at the scale factor array values
TASK: if not already there, make a table of comoving distances from an input array
*/
void ccl_cosmology_distances_from_input(ccl_cosmology * cosmo, int na, double a[], double chi_a[], double E_a[],
   int *status)
{
  double *chi_a_reversed = NULL;
  double *a_reversed = NULL;

  //Do nothing if everything is computed already
  if(cosmo->computed_distances)
    return;

  // Allocate E(a), chi(a) and a(chi) splines
  gsl_spline * E = gsl_spline_alloc(cosmo->spline_params.A_SPLINE_TYPE, na);
  gsl_spline * chi = gsl_spline_alloc(cosmo->spline_params.A_SPLINE_TYPE, na);
  gsl_spline * achi = gsl_spline_alloc(cosmo->spline_params.A_SPLINE_TYPE, na);

  //Check for too little memory
  if (a == NULL || E_a == NULL || chi_a == NULL || E == NULL || chi == NULL || achi == NULL) {
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(
      cosmo, "ccl_background.c: ccl_cosmology_distances_from_input(): ran out of memory\n");
  }

  // Update the minimum scale factor value used by the user, which is necessary because it is used in ccl_tracers.c
  cosmo->spline_params.A_SPLINE_MIN = a[0];

  // Initialize a E(a) spline
  if (!*status){
    if (gsl_spline_init(E, a, E_a, na)){
      *status = CCL_ERROR_SPLINE;
      ccl_cosmology_set_status_message(
        cosmo, "ccl_background.c: ccl_cosmology_distances_from_input(): Error creating  E(a) spline\n");
    }
  }

  // Initialize chi(a) spline
  if (!*status){
    if (gsl_spline_init(chi, a, chi_a, na)){//in Mpc
      *status = CCL_ERROR_SPLINE;
      ccl_cosmology_set_status_message(
        cosmo, "ccl_background.c: ccl_cosmology_distances_from_input(): Error creating  chi(a) spline\n");
    }
  }

  // Reverse the order of chi(a) so that we can initialize the a(chi) spline, which needs monotonically decreasing x-array.
  if (!*status) {
    chi_a_reversed = malloc(na*sizeof(double));
    a_reversed = malloc(na*sizeof(double));
    if (chi_a_reversed == NULL || a_reversed == NULL) {
      *status=CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(
        cosmo, "ccl_background.c: ccl_cosmology_distances_from_input(): ran out of memory\n");
    }
    else {
      for (int i=0; i<na; i++) {
        chi_a_reversed[i] = chi_a[na-1-i];
        a_reversed[i] = a[na-1-i];
      }
    }
  }

  // Initialize a(chi) spline
  if (!*status){
    if (gsl_spline_init(achi, chi_a_reversed, a_reversed, na)){
      *status = CCL_ERROR_SPLINE;
      ccl_cosmology_set_status_message(
        cosmo, "ccl_background.c: ccl_cosmology_distances_from_input(): Error creating  chi(a) spline\n");
    }
  }
  free(chi_a_reversed);
  free(a_reversed);

  if (*status){ //If there was an error, free the GSL splines and return
    gsl_spline_free(E); // Note: you are allowed to call gsl_free() on NULL
    gsl_spline_free(chi);
    gsl_spline_free(achi);
    E = NULL;
    chi = NULL;
    achi = NULL;
  }

  if (*status == 0) {
    //If there were no errors, attach the splines to the cosmo struct and end the function.
    cosmo->data.E             = E;
    cosmo->data.chi           = chi;
    cosmo->data.achi          = achi;
    cosmo->computed_distances = true;
  }

  return;
}

/* ----- ROUTINE: ccl_cosmology_growth_from_input ------
INPUT: cosmology, scale factor array, growth array, growth rate array
TASK: if not already there, create growth splines with the input arrays and store them.
*/
void ccl_cosmology_growth_from_input(ccl_cosmology* cosmo, int na, double a[],
                                     double growth_arr[], double fgrowth_arr[], int* status)
{
  int chistatus;
  if (cosmo->computed_growth)
    return;

  double *growth_normed = NULL;
  gsl_spline * growth = gsl_spline_alloc(cosmo->spline_params.A_SPLINE_TYPE, na);
  gsl_spline * fgrowth = gsl_spline_alloc(cosmo->spline_params.A_SPLINE_TYPE, na);
  //The last element corresponds to a=1 (which is checked for in python).
  double growth0 = growth_arr[na-1];

  if (growth == NULL || fgrowth == NULL) {
    *status = CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(
        cosmo, "ccl_background.c: ccl_cosmology_growth_from_input(): ran out of memory\n");
  }

  // Need to normalize the growth factor
  if (*status == 0) {
    growth_normed = malloc(na*sizeof(double));
    if (growth_normed == NULL) {
      *status = CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(
           cosmo, "ccl_background.c: ccl_cosmology_growth_from_input(): ran out of memory\n");
    }
  }

  if (*status == 0) {
    for(int ia=0; ia<na; ia++)
      growth_normed[ia] = growth_arr[ia]/growth0;

    chistatus = gsl_spline_init(growth, a, growth_normed, na);
    if (chistatus) {
      *status = CCL_ERROR_SPLINE;
      ccl_cosmology_set_status_message(
        cosmo, "ccl_background.c: ccl_cosmology_growth_from_input(): Error creating D(a) spline\n");
    }
  }

  if (*status == 0) {
    chistatus = gsl_spline_init(fgrowth, a, fgrowth_arr, na);
    if (chistatus) {
      *status = CCL_ERROR_SPLINE;
      ccl_cosmology_set_status_message(
        cosmo, "ccl_background.c: ccl_cosmology_growth_from_input(): Error creating f(a) spline\n");
    }
  }

  if (*status){ //If there was an error, free the GSL splines and return
    gsl_spline_free(growth);
    gsl_spline_free(fgrowth);
    growth = NULL;
    fgrowth = NULL;
  }

  if (*status == 0) {
    // assign all the splines we've just made to the structure.
    cosmo->data.growth = growth;
    cosmo->data.fgrowth = fgrowth;
    cosmo->data.growth0 = growth0;
    cosmo->computed_growth = true;
  }

  free(growth_normed);
}

/* ----- ROUTINE: ccl_cosmology_compute_growth ------
INPUT: cosmology
TASK: if not already there, make a table of growth function and growth rate
      normalize growth to input parameter growth0
*/
void ccl_cosmology_compute_growth(ccl_cosmology* cosmo, int* status)
{
  if (cosmo->computed_growth)
    return;

  // Create logarithmically and then linearly-spaced values of the scale factor
  int chistatus = 0, na = cosmo->spline_params.A_SPLINE_NA+cosmo->spline_params.A_SPLINE_NLOG-1;
  double *a = NULL;
  gsl_integration_cquad_workspace * workspace = NULL;
  gsl_function F;
  gsl_spline *df_a_spline = NULL;
  double *df_arr = NULL;
  gsl_spline *df_z_spline = NULL;
  int status_mg = 0, gslstatus;
  double growth0, fgrowth0;
  double *y = NULL;
  double *y2 = NULL;
  double df, integ;

  if (*status == 0) {
    a = ccl_linlog_spacing(
      cosmo->spline_params.A_SPLINE_MINLOG, cosmo->spline_params.A_SPLINE_MIN,
      cosmo->spline_params.A_SPLINE_MAX, cosmo->spline_params.A_SPLINE_NLOG,
      cosmo->spline_params.A_SPLINE_NA);

    if (a == NULL ||
        (fabs(a[0]-cosmo->spline_params.A_SPLINE_MINLOG) > 1e-5) ||
        (fabs(a[na-1]-cosmo->spline_params.A_SPLINE_MAX) > 1e-5) ||
        (a[na-1] > 1.0)) {
      *status = CCL_ERROR_LINSPACE;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_background.c: ccl_cosmology_compute_growth(): "
        "Error creating logarithmically and then linear spacing in a\n");
    }
  }

  if(cosmo->params.has_mgrowth) {
    if (*status == 0) {
      df_arr = malloc(na*sizeof(double));

      if(df_arr == NULL) {
        *status = CCL_ERROR_MEMORY;
        ccl_cosmology_set_status_message(
          cosmo, "ccl_background.c: ccl_cosmology_compute_growth(): ran out of memory\n");
      }
    }

    // Generate spline for Delta f(z) that we will then interpolate into an array of a
    if (*status == 0) {
      df_z_spline = gsl_spline_alloc(cosmo->spline_params.A_SPLINE_TYPE, cosmo->params.nz_mgrowth);
      if (df_z_spline == NULL) {
        *status = CCL_ERROR_MEMORY;
        ccl_cosmology_set_status_message(
          cosmo, "ccl_background.c: ccl_cosmology_compute_growth(): ran out of memory\n");
      }
    }

    if (*status == 0) {
      chistatus = gsl_spline_init(df_z_spline,cosmo->params.z_mgrowth,
                                  cosmo->params.df_mgrowth,
                                  cosmo->params.nz_mgrowth);

      if(chistatus) {
        *status = CCL_ERROR_SPLINE;
        ccl_cosmology_set_status_message(
          cosmo, "ccl_background.c: ccl_cosmology_compute_growth(): Error creating Delta f(z) spline\n");
      }
    }

    if (*status == 0) {
      for (int i=0; i<na; i++) {
        if(a[i]>0) {
          double z=1./a[i]-1.;

          if(z<=cosmo->params.z_mgrowth[0])
            df_arr[i]=cosmo->params.df_mgrowth[0];
          else if(z>cosmo->params.z_mgrowth[cosmo->params.nz_mgrowth-1])
            df_arr[i]=cosmo->params.df_mgrowth[cosmo->params.nz_mgrowth-1];
          else
              chistatus |= gsl_spline_eval_e(df_z_spline,z,NULL,&df_arr[i]);
        } else {
            df_arr[i]=0;
        }
      }
      if(chistatus) {
        *status = CCL_ERROR_SPLINE;
        ccl_cosmology_set_status_message(
          cosmo, "ccl_background.c: ccl_cosmology_compute_growth(): Error evaluating Delta f(z) spline\n");
      }
    }

    // Generate Delta(f) spline
    if (*status == 0) {
      df_a_spline = gsl_spline_alloc(cosmo->spline_params.A_SPLINE_TYPE,na);

      if (df_a_spline == NULL) {
        *status = CCL_ERROR_MEMORY;
        ccl_cosmology_set_status_message(
          cosmo, "ccl_background.c: ccl_cosmology_compute_growth(): ran out of memory\n");
      }
    }

    if (*status == 0) {
      chistatus = gsl_spline_init(df_a_spline, a, df_arr, na);

      if (chistatus) {
        *status = CCL_ERROR_SPLINE;
        ccl_cosmology_set_status_message(
          cosmo, "ccl_background.c: ccl_cosmology_compute_growth(): Error creating Delta f(a) spline\n");
      }
    }

    if (*status == 0) {
      workspace = gsl_integration_cquad_workspace_alloc(cosmo->gsl_params.N_ITERATION);
      F.function = &df_integrand;
      F.params = df_a_spline;
    }
  }

  // allocate space for y, which will be all three
  // of E(a), chi(a), D(a) and f(a) in turn.
  if (*status == 0) {
    y = malloc(sizeof(double)*na);
    if (y == NULL) {
      *status = CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(
          cosmo, "ccl_background.c: ccl_cosmology_compute_growth(): ran out of memory\n");
    }
  }

  if (*status == 0) {
    y2 = malloc(sizeof(double)*na);
    if (y2 == NULL) {
      *status = CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(
          cosmo, "ccl_background.c: ccl_cosmology_compute_growth(): ran out of memory\n");
    }
  }

  if (*status == 0) {
    // Get the growth factor and growth rate at z=0
    chistatus |= growth_factor_and_growth_rate(1., &growth0, &fgrowth0, cosmo, status);

    // Get the growth factor and growth rate at other redshifts
    for(int i=0; i<na; i++) {
      chistatus |= growth_factor_and_growth_rate(a[i], &(y[i]), &(y2[i]), cosmo, status);

      if(cosmo->params.has_mgrowth) {
        if(a[i]>0) {
          // Add modification to f
          gslstatus = gsl_spline_eval_e(df_a_spline, a[i], NULL, &df);
          if (gslstatus != GSL_SUCCESS) {
            ccl_raise_gsl_warning(gslstatus, "ccl_background.c: ccl_cosmology_compute_growth():");
            status_mg |= gslstatus;
          }
          y2[i] += df;

          // Multiply D by exp(-int(df))
          gslstatus = gsl_integration_cquad(
            &F, a[i], 1.0, 0.0, cosmo->gsl_params.INTEGRATION_DISTANCE_EPSREL,
            workspace, &integ, NULL, NULL);
          if (gslstatus != GSL_SUCCESS) {
            ccl_raise_gsl_warning(gslstatus, "ccl_background.c: ccl_cosmology_compute_growth():");
            status_mg |= gslstatus;
          }
          y[i] *= exp(-integ);
        }
      }
      // Normalizing to the growth factor to the growth today
      y[i] /= growth0;
    }

    if (chistatus || status_mg || *status) {
      if (chistatus) {
        *status = CCL_ERROR_INTEG;
        ccl_cosmology_set_status_message(
          cosmo, "ccl_background.c: ccl_cosmology_compute_growth(): integral for linear growth factor didn't converge\n");
      }
      if(status_mg) {
        *status = CCL_ERROR_INTEG;
        ccl_cosmology_set_status_message(
          cosmo, "ccl_background.c: ccl_cosmology_compute_growth(): integral for MG growth factor didn't converge\n");
      }
    }
  }

  gsl_spline *growth = NULL;
  gsl_spline *fgrowth = NULL;

  if (*status == 0) {
    growth = gsl_spline_alloc(cosmo->spline_params.A_SPLINE_TYPE, na);
    fgrowth = gsl_spline_alloc(cosmo->spline_params.A_SPLINE_TYPE, na);

    if (growth == NULL || fgrowth == NULL) {
      *status = CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(
          cosmo, "ccl_background.c: ccl_cosmology_compute_growth(): ran out of memory\n");
    }
  }

  if (*status == 0) {
    chistatus = gsl_spline_init(growth, a, y, na);

    if (chistatus) {
      *status = CCL_ERROR_SPLINE;
      ccl_cosmology_set_status_message(
        cosmo, "ccl_background.c: ccl_cosmology_compute_growth(): Error creating D(a) spline\n");
    }
  }

  if (*status == 0) {
    chistatus = gsl_spline_init(fgrowth, a, y2, na);
    if (chistatus) {
      *status = CCL_ERROR_SPLINE;
      ccl_cosmology_set_status_message(
        cosmo, "ccl_background.c: ccl_cosmology_compute_growth(): Error creating f(a) spline\n");
    }
  }

  if (*status == 0) {
    // assign all the splines we've just made to the structure.
    cosmo->data.growth = growth;
    cosmo->data.fgrowth = fgrowth;
    cosmo->data.growth0 = growth0;
    cosmo->computed_growth = true;
  }
  else {
    gsl_spline_free(growth);
    gsl_spline_free(fgrowth);
  }

  free(a);
  free(y);
  free(y2);
  free(df_arr);
  gsl_spline_free(df_z_spline);
  gsl_spline_free(df_a_spline);
  gsl_integration_cquad_workspace_free(workspace);
}

//Expansion rate normalized to 1 today

double ccl_h_over_h0(ccl_cosmology * cosmo, double a, int* status)
{

  if(!cosmo->computed_distances) {
    *status = CCL_ERROR_DISTANCES_INIT;
    ccl_cosmology_set_status_message(
      cosmo,
      "ccl_background.c: ccl_h_over_h0(): distance splines have not been precomputed!");
    return NAN;
  }

  double h_over_h0;
  int gslstatus = gsl_spline_eval_e(cosmo->data.E, a, NULL, &h_over_h0);
  if(gslstatus != GSL_SUCCESS) {
    ccl_raise_gsl_warning(gslstatus, "ccl_background.c: ccl_h_over_h0():");
    *status = gslstatus;
    ccl_cosmology_set_status_message(cosmo, "ccl_background.c: ccl_h_over_h0(): Scale factor outside interpolation range.\n");
    return NAN;
  }

  return h_over_h0;
}

void ccl_h_over_h0s(ccl_cosmology * cosmo, int na, double a[], double output[], int * status)
{
  int _status;

  for (int i=0; i<na; i++) {
    _status = 0;
    output[i] = ccl_h_over_h0(cosmo, a[i], &_status);
    *status |= _status;
  }
}

// Distance-like function examples, all in Mpc
double ccl_comoving_radial_distance(ccl_cosmology * cosmo, double a, int * status)
{
  if((a > (1.0 - 1.e-8)) && (a<=1.0)) {
    return 0.;
  }
  else if(a>1.) {
    *status = CCL_ERROR_COMPUTECHI;
    ccl_cosmology_set_status_message(cosmo, "ccl_background.c: scale factor cannot be larger than 1.\n");
    return NAN;
  }
  else {
    if(!cosmo->computed_distances) {
      *status = CCL_ERROR_DISTANCES_INIT;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_background.c: ccl_comoving_radial_distance(): distance splines have not been precomputed!");
      return NAN;
    }

    double crd;
    int gslstatus = gsl_spline_eval_e(cosmo->data.chi, a, NULL, &crd);
    if(gslstatus != GSL_SUCCESS) {
      ccl_raise_gsl_warning(gslstatus, "ccl_background.c: ccl_comoving_radial_distance():");
      *status = gslstatus;
      ccl_cosmology_set_status_message(
        cosmo, "ccl_background.c: ccl_comoving_radial_distance(): Scale factor outside interpolation range.\n");
      return NAN;
    }
    return crd;
  }
}

void ccl_comoving_radial_distances(ccl_cosmology * cosmo, int na, double a[], double output[], int* status)
{
  int _status;

  for (int i=0; i<na; i++) {
    _status = 0;
    output[i] = ccl_comoving_radial_distance(cosmo, a[i], &_status);
    *status |= _status;
  }
}

double ccl_sinn(ccl_cosmology *cosmo, double chi, int *status)
{
  //////
  //         { sin(x)  , if k==1
  // sinn(x)={  x      , if k==0
  //         { sinh(x) , if k==-1
  switch(cosmo->params.k_sign) {
  case -1:
    return sinh(cosmo->params.sqrtk * chi) / cosmo->params.sqrtk;
  case 1:
    return sin(cosmo->params.sqrtk*chi) / cosmo->params.sqrtk;
  case 0:
    return chi;
  default:
    *status = CCL_ERROR_PARAMETERS;
    ccl_cosmology_set_status_message(cosmo, "ccl_background.c: ccl_sinn: ill-defined cosmo->params.k_sign = %d",
                                     cosmo->params.k_sign);
    return NAN;
  }
}

double ccl_comoving_angular_distance(ccl_cosmology * cosmo, double a, int* status)
{
  if((a > (1.0 - 1.e-8)) && (a<=1.0)) {
    return 0.;
  }
  else if(a>1.) {
    *status = CCL_ERROR_COMPUTECHI;
    ccl_cosmology_set_status_message(cosmo, "ccl_background.c: scale factor cannot be larger than 1.\n");
    return NAN;
  }
  else {
    if (!cosmo->computed_distances) {
      *status = CCL_ERROR_DISTANCES_INIT;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_background.c: ccl_comoving_angular_distance(): distance splines have not been precomputed!");
      return NAN;
    }

    double chi;
    int gslstatus = gsl_spline_eval_e(cosmo->data.chi, a, NULL, &chi);
    if(gslstatus != GSL_SUCCESS) {
      ccl_raise_gsl_warning(gslstatus, "ccl_background.c: ccl_comoving_angular_distance():");
      *status |= gslstatus;
      ccl_cosmology_set_status_message(cosmo, "ccl_background.c: ccl_comoving_angular_distance(): Scale factor outside interpolation range.\n");
      return NAN;
    }
    return ccl_sinn(cosmo,chi,status);
  }
}

void ccl_comoving_angular_distances(ccl_cosmology * cosmo, int na, double a[],
                                    double output[], int* status)
{
  int _status;

  for (int i=0; i < na; i++) {
    _status = 0;
    output[i] = ccl_comoving_angular_distance(cosmo, a[i], &_status);
    *status |= _status;
  }
}

double ccl_angular_diameter_distance(ccl_cosmology * cosmo, double a1, double a2, int* status)
{
  if(a1>1. || a2>1. || a1<a2) {
    *status = CCL_ERROR_COMPUTECHI;
    ccl_cosmology_set_status_message(cosmo,"ccl_background.c: ccl_angular_diameter_distance(): error on input scale factor.");
    return NAN;
  } else {
      if (!cosmo->computed_distances) {
          *status = CCL_ERROR_DISTANCES_INIT;
          ccl_cosmology_set_status_message(cosmo,"ccl_background.c: ccl_angular_diameter_distance(): distance splines have not been precomputed!");
          return NAN;
      }
      double chi1,chi2;
      int gslstatus = gsl_spline_eval_e(cosmo->data.chi, a1, NULL, &chi1);
      if(gslstatus != GSL_SUCCESS) {
          ccl_raise_gsl_warning(gslstatus, "ccl_background.c: ccl_angular_diameter_distance():");
          *status |= gslstatus;
          ccl_cosmology_set_status_message(cosmo,"ccl_background.c: ccl_angular_diameter_distance(): Scale factor outside interpolation range.\n");
          return NAN;
      }
      gslstatus = gsl_spline_eval_e(cosmo->data.chi, a2, NULL, &chi2);
      if(gslstatus != GSL_SUCCESS) {
          ccl_raise_gsl_warning(gslstatus, "ccl_background.c: ccl_angular_diameter_distance():");
          *status |= gslstatus;
          ccl_cosmology_set_status_message(cosmo,"ccl_background.c: ccl_angular_diameter_distance(): Scale factor outside interpolation range.\n");
          return NAN;
      }
      return a2*ccl_sinn(cosmo,chi2-chi1,status);
  }
}

void ccl_angular_diameter_distances(ccl_cosmology * cosmo, int na, double a1[], double a2[],
                                    double output[], int* status)
{
  int _status;

  for (int i=0; i < na; i++) {
    _status = 0;
    output[i] = ccl_angular_diameter_distance(cosmo, a1[i], a2[i], &_status);
    *status |= _status;
  }
}

double ccl_luminosity_distance(ccl_cosmology * cosmo, double a, int* status)
{
  return ccl_comoving_angular_distance(cosmo, a, status) / a;
}

void ccl_luminosity_distances(ccl_cosmology * cosmo, int na, double a[], double output[], int * status)
{
  int _status;

  for (int i=0; i<na; i++) {
    _status = 0;
    output[i] = ccl_luminosity_distance(cosmo, a[i], &_status);
    *status |= _status;
  }
}

double ccl_distance_modulus(ccl_cosmology * cosmo, double a, int* status)
{
  if((a > (1.0 - 1.e-8)) && (a<=1.0)) {
    *status = CCL_ERROR_COMPUTECHI;
    ccl_cosmology_set_status_message(cosmo, "ccl_background.c: distance_modulus undefined for a=1.\n");
    return NAN;
  } else if(a>1.) {
    *status = CCL_ERROR_COMPUTECHI;
    ccl_cosmology_set_status_message(cosmo, "ccl_background.c: scale factor cannot be larger than 1.\n");
    return NAN;
  } else {
    if (!cosmo->computed_distances) {
      *status = CCL_ERROR_DISTANCES_INIT;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_background.c: ccl_distance_modulus(): distance splines have not been precomputed!");
      return NAN;
    }
    /* distance modulus = 5 * log10(d) - 5
       Since d in CCL is in Mpc, you get
         5*log10(10^6) - 5 = 30 - 5 = 25
      for the constant.
    */
    return 5 * log10(ccl_luminosity_distance(cosmo, a, status)) + 25;
  }
}

void ccl_distance_moduli(ccl_cosmology * cosmo, int na, double a[], double output[], int * status)
{
  int _status;

  for (int i=0; i<na; i++) {
    _status = 0;
    output[i] = ccl_distance_modulus(cosmo, a[i], &_status);
    *status |= _status;
  }
}

//Scale factor for a given distance
double ccl_scale_factor_of_chi(ccl_cosmology * cosmo, double chi, int * status)
{
  if((chi < 1.e-8) && (chi>=0.)) {
    return 1.;
  }
  else if(chi<0.) {
    *status = CCL_ERROR_COMPUTECHI;
    ccl_cosmology_set_status_message(cosmo, "ccl_background.c: distance cannot be smaller than 0.\n");
    return NAN;
  }
  else {
    if (!cosmo->computed_distances) {
      *status = CCL_ERROR_DISTANCES_INIT;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_background.c: ccl_scale_factor_of_chi(): distance splines have not been precomputed!");
      return NAN;
    }
    double a;
    int gslstatus = gsl_spline_eval_e(cosmo->data.achi, chi, NULL, &a);
    if(gslstatus != GSL_SUCCESS) {
      ccl_raise_gsl_warning(gslstatus, "ccl_background.c: ccl_scale_factor_of_chi():");
      *status |= gslstatus;
    }
    return a;
  }
}

void ccl_scale_factor_of_chis(ccl_cosmology * cosmo, int nchi, double chi[], double output[], int * status)
{
  int _status;

  for (int i=0; i<nchi; i++) {
    _status = 0;
    output[i] = ccl_scale_factor_of_chi(cosmo, chi[i], &_status);
    *status |= _status;
  }
}

double ccl_growth_factor(ccl_cosmology * cosmo, double a, int * status)
{
  if(a==1.){
    return 1.;
  }
  else if(a>1.) {
    *status = CCL_ERROR_COMPUTECHI;
    ccl_cosmology_set_status_message(
      cosmo, "ccl_background.c: scale factor cannot be larger than 1.");
    return NAN;
  }
  else {
    if (!cosmo->computed_growth){
      *status = CCL_ERROR_GROWTH_INIT;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_background.c: ccl_growth_factor(): growth factor splines have not been precomputed!");
      return NAN;
    }
    if (*status != CCL_ERROR_NOT_IMPLEMENTED) {
      double D;
      int gslstatus = gsl_spline_eval_e(cosmo->data.growth, a, NULL, &D);
      if(gslstatus != GSL_SUCCESS) {
        ccl_raise_gsl_warning(gslstatus, "ccl_background.c: ccl_growth_factor():");
        *status |= gslstatus;
        ccl_cosmology_set_status_message(
          cosmo, "ccl_background.c: ccl_growth_factor(): Scale factor outside interpolation range.");
        return NAN;
      }
      return D;
    }
    else {
      return NAN;
    }
  }
}

void ccl_growth_factors(ccl_cosmology * cosmo, int na, double a[], double output[], int * status)
{
  int _status;

  for (int i=0; i<na; i++) {
    _status = 0;
    output[i] = ccl_growth_factor(cosmo, a[i], &_status);
    *status |= _status;
  }
}

double ccl_growth_factor_unnorm(ccl_cosmology * cosmo, double a, int * status)
{

  if (!cosmo->computed_growth) {
    *status = CCL_ERROR_GROWTH_INIT;
    ccl_cosmology_set_status_message(
      cosmo,
      "ccl_background.c: ccl_growth_factor_unnorm(): growth factor splines have not been precomputed!");
  }

  if (*status != CCL_ERROR_NOT_IMPLEMENTED) {
      return cosmo->data.growth0 * ccl_growth_factor(cosmo, a, status);
  } else {
    return NAN;
  }
}

void ccl_growth_factors_unnorm(ccl_cosmology * cosmo, int na, double a[], double output[], int * status)
{
  int _status;

  for (int i=0; i<na; i++) {
    _status = 0;
    output[i] = ccl_growth_factor_unnorm(cosmo, a[i], &_status);
    *status |= _status;
  }
}

double ccl_growth_rate(ccl_cosmology * cosmo, double a, int * status)
{
  if(a>1.) {
    *status = CCL_ERROR_COMPUTECHI;
    ccl_cosmology_set_status_message(cosmo, "ccl_background.c: scale factor cannot be larger than 1.\n");
    return NAN;
  } else {
    if (!cosmo->computed_growth) {
      *status = CCL_ERROR_GROWTH_INIT;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_background.c: ccl_growth_rate(): growth rate splines have not been precomputed!");
    }
    if(*status != CCL_ERROR_NOT_IMPLEMENTED) {
      double g;
      int gslstatus = gsl_spline_eval_e(cosmo->data.fgrowth, a, NULL ,&g);
      if(gslstatus != GSL_SUCCESS) {
        ccl_raise_gsl_warning(gslstatus, "ccl_background.c: ccl_growth_rate():");
        *status |= gslstatus;
        ccl_cosmology_set_status_message(cosmo, "ccl_background.c: ccl_growth_rate(): Scale factor outside interpolation range.\n");
        return NAN;
      }
      return g;
    } else {
        return NAN;
    }
  }
}

void ccl_growth_rates(ccl_cosmology * cosmo, int na, double a[], double output[], int * status)
{
  int _status;

  for (int i=0; i<na; i++) {
    _status = 0;
    output[i] = ccl_growth_rate(cosmo, a[i], &_status);
    *status |= _status;
  }
}

ext: c

max_stars_repo_name: borisbolliet/CCL

max_issues_repo_name: borisbolliet/CCL

max_forks_repo_name: borisbolliet/CCL


text:
// Copyright (c) 2013-2017 Anton Kozhevnikov, Thomas Schulthess
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without modification, are permitted provided that
// the following conditions are met:
//
// 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the
//    following disclaimer.
// 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions
//    and the following disclaimer in the documentation and/or other materials provided with the distribution.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED
// WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
// PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
// ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
// CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

/** \file sht.h
 *
 *  \brief Contains declaration and particular implementation of sirius::SHT class.
 */

#ifndef __SHT_H__
#define __SHT_H__

#include <math.h>
#include <stddef.h>
#include <gsl/gsl_sf_coupling.h>
#include <gsl/gsl_sf_legendre.h>
#include <string.h>
#include <vector>
#include "typedefs.h"
#include "utils.h"
#include "linalg.hpp"
#include "LebedevLaikov.h"

namespace sirius
{

/// Spherical harmonics transformation.
class SHT // TODO: better name
{
    private:

        int lmax_;

        int lmmax_;

        int num_points_;

        mdarray<double, 2> coord_;

        mdarray<double, 2> tp_;

        std::vector<double> w_;

        /// backward transformation from Ylm to spherical coordinates
        mdarray<double_complex, 2> ylm_backward_;

        /// forward transformation from spherical coordinates to Ylm
        mdarray<double_complex, 2> ylm_forward_;

        /// backward transformation from Rlm to spherical coordinates
        mdarray<double, 2> rlm_backward_;

        /// forward transformation from spherical coordinates to Rlm
        mdarray<double, 2> rlm_forward_;

        int mesh_type_;

    public:

        /// Default constructor.
        SHT(int lmax__)
            : lmax_(lmax__), mesh_type_(0)
        {
            lmmax_ = (lmax_ + 1) * (lmax_ + 1);

            if (mesh_type_ == 0) num_points_ = Lebedev_Laikov_npoint(2 * lmax_);
            if (mesh_type_ == 1) num_points_ = lmmax_;

            std::vector<double> x(num_points_);
            std::vector<double> y(num_points_);
            std::vector<double> z(num_points_);

            coord_ = mdarray<double, 2>(3, num_points_);

            tp_ = mdarray<double, 2>(2, num_points_);

            w_.resize(num_points_);

            if (mesh_type_ == 0) Lebedev_Laikov_sphere(num_points_, &x[0], &y[0], &z[0], &w_[0]);
            if (mesh_type_ == 1) uniform_coverage();

            ylm_backward_ = mdarray<double_complex, 2>(lmmax_, num_points_);

            ylm_forward_ = mdarray<double_complex, 2>(num_points_, lmmax_);

            rlm_backward_ = mdarray<double, 2>(lmmax_, num_points_);

            rlm_forward_ = mdarray<double, 2>(num_points_, lmmax_);

            for (int itp = 0; itp < num_points_; itp++)
            {
                if (mesh_type_ == 0)
                {
                    coord_(0, itp) = x[itp];
                    coord_(1, itp) = y[itp];
                    coord_(2, itp) = z[itp];

                    auto vs = spherical_coordinates(vector3d<double>(x[itp], y[itp], z[itp]));
                    spherical_harmonics(lmax_, vs[1], vs[2], &ylm_backward_(0, itp));
                    spherical_harmonics(lmax_, vs[1], vs[2], &rlm_backward_(0, itp));
                    for (int lm = 0; lm < lmmax_; lm++)
                    {
                        ylm_forward_(itp, lm) = conj(ylm_backward_(lm, itp)) * w_[itp] * fourpi;
                        rlm_forward_(itp, lm) = rlm_backward_(lm, itp) * w_[itp] * fourpi;
                    }
                }
                if (mesh_type_ == 1)
                {
                    double t = tp_(0, itp);
                    double p = tp_(1, itp);

                    coord_(0, itp) = sin(t) * cos(p);
                    coord_(1, itp) = sin(t) * sin(p);
                    coord_(2, itp) = cos(t);

                    spherical_harmonics(lmax_, t, p, &ylm_backward_(0, itp));
                    spherical_harmonics(lmax_, t, p, &rlm_backward_(0, itp));

                    for (int lm = 0; lm < lmmax_; lm++)
                    {
                        ylm_forward_(lm, itp) = ylm_backward_(lm, itp);
                        rlm_forward_(lm, itp) = rlm_backward_(lm, itp);
                    }
                }
            }

            if (mesh_type_ == 1)
            {
                linalg<CPU>::geinv(lmmax_, ylm_forward_);
                linalg<CPU>::geinv(lmmax_, rlm_forward_);
            }

            #if (__VERIFICATION > 0)
            {
                double dr = 0;
                double dy = 0;

                for (int lm = 0; lm < lmmax_; lm++)
                {
                    for (int lm1 = 0; lm1 < lmmax_; lm1++)
                    {
                        double t = 0;
                        double_complex zt(0, 0);
                        for (int itp = 0; itp < num_points_; itp++)
                        {
                            zt += ylm_forward_(itp, lm) * ylm_backward_(lm1, itp);
                            t += rlm_forward_(itp, lm) * rlm_backward_(lm1, itp);
                        }

                        if (lm == lm1)
                        {
                            zt -= 1.0;
                            t -= 1.0;
                        }
                        dr += std::abs(t);
                        dy += std::abs(zt);
                    }
                }
                dr = dr / lmmax_ / lmmax_;
                dy = dy / lmmax_ / lmmax_;

                if (dr > 1e-15 || dy > 1e-15)
                {
                    std::stringstream s;
                    s << "spherical mesh error is too big" << std::endl
                      << "  real spherical integration error " << dr << std::endl
                      << "  complex spherical integration error " << dy;
                    WARNING(s.str())
                }

                std::vector<double> flm(lmmax_);
                std::vector<double> ftp(num_points_);
                for (int lm = 0; lm < lmmax_; lm++)
                {
                    std::memset(&flm[0], 0, lmmax_ * sizeof(double));
                    flm[lm] = 1.0;
                    backward_transform(lmmax_, &flm[0], 1, lmmax_, &ftp[0]);
                    forward_transform(&ftp[0], 1, lmmax_, lmmax_, &flm[0]);
                    flm[lm] -= 1.0;

                    double t = 0.0;
                    for (int lm1 = 0; lm1 < lmmax_; lm1++) t += std::abs(flm[lm1]);

                    t /= lmmax_;

                    if (t > 1e-15)
                    {
                        std::stringstream s;
                        s << "test of backward / forward real SHT failed" << std::endl
                          << "  total error " << t;
                        WARNING(s.str());
                    }
                }
            }
            #endif
        }

        /// Perform a backward transformation from spherical harmonics to spherical coordinates.
        /** \f[
         *      f(\theta, \phi, r) = \sum_{\ell m} f_{\ell m}(r) Y_{\ell m}(\theta, \phi)
         *  \f]
         *
         *  \param [in] ld Size of leading dimension of flm.
         *  \param [in] flm Raw pointer to \f$ f_{\ell m}(r) \f$.
         *  \param [in] nr Number of radial points.
         *  \param [in] lmmax Maximum number of lm- harmonics to take into sum.
         *  \param [out] ftp Raw pointer to \f$ f(\theta, \phi, r) \f$.
         */
        template <typename T>
        void backward_transform(int ld, T const* flm, int nr, int lmmax, T* ftp);

        /// Perform a forward transformation from spherical coordinates to spherical harmonics.
        /** \f[
         *      f_{\ell m}(r) = \iint  f(\theta, \phi, r) Y_{\ell m}^{*}(\theta, \phi) \sin \theta d\phi d\theta =
         *        \sum_{i} f(\theta_i, \phi_i, r) Y_{\ell m}^{*}(\theta_i, \phi_i) w_i
         *  \f]
         *
         *  \param [in] ftp Raw pointer to \f$ f(\theta, \phi, r) \f$.
         *  \param [in] nr Number of radial points.
         *  \param [in] lmmax Maximum number of lm- coefficients to generate.
         *  \param [in] ld Size of leading dimension of flm.
         *  \param [out] flm Raw pointer to \f$ f_{\ell m}(r) \f$.
         */
        template <typename T>
        void forward_transform(T const* ftp, int nr, int lmmax, int ld, T* flm);

        /// Convert form Rlm to Ylm representation.
        static void convert(int lmax__, double const* f_rlm__, double_complex* f_ylm__);

        /// Convert from Ylm to Rlm representation.
        static void convert(int lmax__, double_complex const* f_ylm__, double* f_rlm__);

        //void rlm_forward_iterative_transform(double *ftp__, int lmmax, int ncol, double* flm)
        //{
        //    Timer t("sirius::SHT::rlm_forward_iterative_transform");
        //
        //    assert(lmmax <= lmmax_);

        //    mdarray<double, 2> ftp(ftp__, num_points_, ncol);
        //    mdarray<double, 2> ftp1(num_points_, ncol);
        //
        //    blas<cpu>::gemm(1, 0, lmmax, ncol, num_points_, 1.0, &rlm_forward_(0, 0), num_points_, &ftp(0, 0), num_points_, 0.0,
        //                    flm, lmmax);
        //
        //    for (int i = 0; i < 2; i++)
        //    {
        //        rlm_backward_transform(flm, lmmax, ncol, &ftp1(0, 0));
        //        double tdiff = 0.0;
        //        for (int ir = 0; ir < ncol; ir++)
        //        {
        //            for (int itp = 0; itp < num_points_; itp++)
        //            {
        //                ftp1(itp, ir) = ftp(itp, ir) - ftp1(itp, ir);
        //                //tdiff += fabs(ftp1(itp, ir));
        //            }
        //        }
        //
        //        for (int itp = 0; itp < num_points_; itp++)
        //        {
        //            tdiff += fabs(ftp1(itp, ncol - 1));
        //        }
        //        std::cout << "iter : " << i << " avg. MT diff = " << tdiff / num_points_ << std::endl;
        //        blas<cpu>::gemm(1, 0, lmmax, ncol, num_points_, 1.0, &rlm_forward_(0, 0), num_points_, &ftp1(0, 0), num_points_, 1.0,
        //                        flm, lmmax);
        //    }
        //}

        /// Transform Cartesian coordinates [x,y,z] to spherical coordinates [r,theta,phi]
        static vector3d<double> spherical_coordinates(vector3d<double> vc);

        /// Generate complex spherical harmonics Ylm
        static void spherical_harmonics(int lmax, double theta, double phi, double_complex* ylm);

        /// Generate real spherical harmonics Rlm
        /** Mathematica code:
         *  \verbatim
         *  R[l_, m_, th_, ph_] :=
         *   If[m > 0, std::sqrt[2]*ComplexExpand[Re[SphericalHarmonicY[l, m, th, ph]]],
         *   If[m < 0, std::sqrt[2]*ComplexExpand[Im[SphericalHarmonicY[l, m, th, ph]]],
         *   If[m == 0, ComplexExpand[Re[SphericalHarmonicY[l, 0, th, ph]]]]]]
         *  \endverbatim
         */
        static void spherical_harmonics(int lmax, double theta, double phi, double* rlm);

        /// Compute element of the transformation matrix from complex to real spherical harmonics.
        /** Real spherical harmonic can be written as a linear combination of complex harmonics:
         *
         *    \f[
         *        R_{\ell m}(\theta, \phi) = \sum_{m'} a^{\ell}_{m' m}Y_{\ell m'}(\theta, \phi)
         *    \f]
         *    where
         *    \f[
         *        a^{\ell}_{m' m} = \langle Y_{\ell m'} | R_{\ell m} \rangle
         *    \f]
         *
         *    Transformation from real to complex spherical harmonics is conjugate transpose:
         *
         *    \f[
         *        Y_{\ell m}(\theta, \phi) = \sum_{m'} a^{\ell*}_{m m'}R_{\ell m'}(\theta, \phi)
         *    \f]
         *
         *    Mathematica code:
         *    \verbatim
         *    b[m1_, m2_] :=
         *     If[m1 == 0, 1,
         *     If[m1 < 0 && m2 < 0, -I/std::sqrt[2],
         *     If[m1 > 0 && m2 < 0, (-1)^m1*I/std::sqrt[2],
         *     If[m1 < 0 && m2 > 0, (-1)^m2/std::sqrt[2],
         *     If[m1 > 0 && m2 > 0, 1/std::sqrt[2]]]]]]
         *
         *    a[m1_, m2_] := If[Abs[m1] == Abs[m2], b[m1, m2], 0]
         *
         *    R[l_, m_, t_, p_] := Sum[a[m1, m]*SphericalHarmonicY[l, m1, t, p], {m1, -l, l}]
         *    \endverbatim
         */
        static inline double_complex ylm_dot_rlm(int l, int m1, int m2)
        {
            const double isqrt2 = 0.70710678118654752440;

            assert(l >= 0 && std::abs(m1) <= l && std::abs(m2) <= l);

            if (!((m1 == m2) || (m1 == -m2))) {
                return double_complex(0, 0);
            }

            if (m1 == 0) {
                return double_complex(1, 0);
            }

            if (m1 < 0) {
                if (m2 < 0) {
                    return -double_complex(0, isqrt2);
                } else {
                    return std::pow(-1.0, m2) * double_complex(isqrt2, 0);
                }
            } else {
                if (m2 < 0) {
                    return pow(-1.0, m1) * double_complex(0, isqrt2);
                } else {
                    return double_complex(isqrt2, 0);
                }
            }
        }

        static inline double_complex rlm_dot_ylm(int l, int m1, int m2)
        {
            return std::conj(ylm_dot_rlm(l, m2, m1));
        }

        /// Gaunt coefficent of three complex spherical harmonics.
        /**
         *  \f[
         *    \langle Y_{\ell_1 m_1} | Y_{\ell_2 m_2} | Y_{\ell_3 m_3} \rangle
         *  \f]
         */
        static double gaunt_ylm(int l1, int l2, int l3, int m1, int m2, int m3);

        /// Gaunt coefficent of three real spherical harmonics.
        /**
         *  \f[
         *    \langle R_{\ell_1 m_1} | R_{\ell_2 m_2} | R_{\ell_3 m_3} \rangle
         *  \f]
         */
        static double gaunt_rlm(int l1, int l2, int l3, int m1, int m2, int m3);

        /// Gaunt coefficent of two complex and one real spherical harmonics.
        /**
         *  \f[
         *    \langle Y_{\ell_1 m_1} | R_{\ell_2 m_2} | Y_{\ell_3 m_3} \rangle
         *  \f]
         */
        static double_complex gaunt_hybrid(int l1, int l2, int l3, int m1, int m2, int m3);

        void uniform_coverage();

        /// Return Clebsch-Gordan coefficient.
        /** Clebsch-Gordan coefficients arise when two angular momenta are combined into a
         *  total angular momentum.
         */
        static inline double clebsch_gordan(int l1, int l2, int l3, int m1, int m2, int m3)
        {
            assert(l1 >= 0);
            assert(l2 >= 0);
            assert(l3 >= 0);
            assert(m1 >= -l1 && m1 <= l1);
            assert(m2 >= -l2 && m2 <= l2);
            assert(m3 >= -l3 && m3 <= l3);

            return pow(-1, l1 - l2 + m3) * sqrt(double(2 * l3 + 1)) *
                   gsl_sf_coupling_3j(2 * l1, 2 * l2, 2 * l3, 2 * m1, 2 * m2, -2 * m3);
        }

        inline double_complex ylm_backward(int lm,  int itp)
        {
            return ylm_backward_(lm, itp);
        }

        inline double rlm_backward(int lm,  int itp)
        {
            return rlm_backward_(lm, itp);
        }

        inline double coord(int x, int itp)
        {
            return coord_(x, itp);
        }

        inline int num_points()
        {
            return num_points_;
        }

        inline int lmax()
        {
            return lmax_;
        }

        inline int lmmax()
        {
            return lmmax_;
        }

        static void wigner_d_matrix(int l, double beta, mdarray<double, 2>& d_mtrx__)
        {
            long double cos_b2 = std::cos((long double)beta / 2.0L);
            long double sin_b2 = std::sin((long double)beta / 2.0L);

            for (int m1 = -l; m1 <= l; m1++)
            {
                for (int m2 = -l; m2 <= l; m2++)
                {
                    long double d = 0;
                    for (int j = 0; j <= std::min(l + m1, l - m2); j++)
                    {
                        if ((l - m2 - j) >= 0 && (l + m1 - j) >= 0 && (j + m2 - m1) >= 0)
                        {
                            long double g = (std::sqrt(Utils::factorial(l + m1)) / Utils::factorial(l - m2 - j)) *
                                            (std::sqrt(Utils::factorial(l - m1)) / Utils::factorial(l + m1 - j)) *
                                            (std::sqrt(Utils::factorial(l - m2)) / Utils::factorial(j + m2 - m1)) *
                                            (std::sqrt(Utils::factorial(l + m2)) / Utils::factorial(j));
                            d += g * std::pow(-1, j) * std::pow(cos_b2, 2 * l + m1 - m2 - 2 * j) * std::pow(sin_b2, 2 * j + m2 - m1);
                        }
                    }
                    d_mtrx__(m1 + l, m2 + l) = (double)d;
                }
            }
        }

        static void rotation_matrix_l(int l, vector3d<double> euler_angles, int proper_rotation,
                                      double_complex* rot_mtrx__, int ld)
        {
            mdarray<double_complex, 2> rot_mtrx(rot_mtrx__, ld, 2 * l + 1);

            mdarray<double, 2> d_mtrx(2 * l + 1, 2 * l + 1);
            wigner_d_matrix(l, euler_angles[1], d_mtrx);

            for (int m1 = -l; m1 <= l; m1++)
            {
                for (int m2 = -l; m2 <= l; m2++)
                {
                    rot_mtrx(m1 + l, m2 + l) = std::exp(double_complex(0, -euler_angles[0] * m1 - euler_angles[2] * m2)) *
                                               d_mtrx(m1 + l, m2 + l) * std::pow(proper_rotation, l);
                }
            }
        }

        static void rotation_matrix_l(int l, vector3d<double> euler_angles, int proper_rotation,
                                      double* rot_mtrx__, int ld)
        {
            mdarray<double, 2> rot_mtrx_rlm(rot_mtrx__, ld, 2 * l + 1);
            mdarray<double_complex, 2> rot_mtrx_ylm(2 * l + 1, 2 * l + 1);

            mdarray<double, 2> d_mtrx(2 * l + 1, 2 * l + 1);
            wigner_d_matrix(l, euler_angles[1], d_mtrx);

            for (int m1 = -l; m1 <= l; m1++)
            {
                for (int m2 = -l; m2 <= l; m2++)
                {
                    rot_mtrx_ylm(m1 + l, m2 + l) = std::exp(double_complex(0, -euler_angles[0] * m1 - euler_angles[2] * m2)) *
                                                   d_mtrx(m1 + l, m2 + l) * std::pow(proper_rotation, l);
                }
            }
            for (int m1 = -l; m1 <= l; m1++)
            {
                auto i13 = (m1 == 0) ? std::vector<int>({0}) : std::vector<int>({-m1, m1});

                for (int m2 = -l; m2 <= l; m2++)
                {
                    auto i24 = (m2 == 0) ? std::vector<int>({0}) : std::vector<int>({-m2, m2});

                    for (int m3: i13)
                    {
                        for (int m4: i24)
                        {
                            rot_mtrx_rlm(m1 + l, m2 + l) += std::real(rlm_dot_ylm(l, m1, m3) *
                                                                      rot_mtrx_ylm(m3 + l, m4 + l) *
                                                                      ylm_dot_rlm(l, m4, m2));
                        }
                    }
                }
            }
        }

        static void rotation_matrix(int lmax, vector3d<double> euler_angles, int proper_rotation,
                                    mdarray<double_complex, 2>& rotm)
        {
            rotm.zero();

            for (int l = 0; l <= lmax; l++)
            {
                rotation_matrix_l(l, euler_angles, proper_rotation, &rotm(l * l, l * l), rotm.ld());
            }
        }

        static void rotation_matrix(int lmax, vector3d<double> euler_angles, int proper_rotation,
                                    mdarray<double, 2>& rotm)
        {
            rotm.zero();

            for (int l = 0; l <= lmax; l++)
            {
                rotation_matrix_l(l, euler_angles, proper_rotation, &rotm(l * l, l * l), rotm.ld());
            }
        }

        /// Compute derivative of real-spherical harmonic with respect to theta angle.
        static double dRlm_dtheta(int lm, double theta, double phi)
        {
            switch (lm) {
                case 0: return 0;

                case 1: return -(std::sqrt(3/pi)*std::cos(theta)*std::sin(phi))/2.;

                case 2: return -(std::sqrt(3/pi)*std::sin(theta))/2.;

                case 3: return -(std::sqrt(3/pi)*std::cos(phi)*std::cos(theta))/2.;

                case 4: return -(std::sqrt(15/pi)*std::cos(phi)*std::cos(theta)*std::sin(phi)*std::sin(theta));

                case 5: return -(std::sqrt(15/pi)*std::cos(2*theta)*std::sin(phi))/2.;

                case 6: return (-3*std::sqrt(5/pi)*std::cos(theta)*std::sin(theta))/2.;

                case 7: return -(std::sqrt(15/pi)*std::cos(phi)*std::cos(2*theta))/2.;

                case 8: return (std::sqrt(15/pi)*std::cos(2*phi)*std::sin(2*theta))/4.;

                case 9: return (-3*std::sqrt(35/(2.*pi))*std::cos(theta)*std::sin(3*phi)*std::pow(std::sin(theta),2))/4.;

                case 10: return (std::sqrt(105/pi)*std::sin(2*phi)*(std::sin(theta) - 3*std::sin(3*theta)))/16.;

                case 11: return (std::sqrt(21/(2.*pi))*std::cos(theta)*(7 - 15*std::cos(2*theta))*std::sin(phi))/8.;

                case 12: return (-3*std::sqrt(7/pi)*(3 + 5*std::cos(2*theta))*std::sin(theta))/8.;

                case 13: return (std::sqrt(21/(2.*pi))*std::cos(phi)*std::cos(theta)*(7 - 15*std::cos(2*theta)))/8.;

                case 14: return (std::sqrt(105/pi)*std::cos(2*phi)*(1 + 3*std::cos(2*theta))*std::sin(theta))/8.;

                case 15: return (-3*std::sqrt(35/(2.*pi))*std::cos(3*phi)*std::cos(theta)*std::pow(std::sin(theta),2))/4.;

                case 16: return (-3*std::sqrt(35/pi)*std::cos(theta)*std::sin(4*phi)*std::pow(std::sin(theta),3))/4.;

                case 17: return (-3*std::sqrt(35/(2.*pi))*(1 + 2*std::cos(2*theta))*std::sin(3*phi)*std::pow(std::sin(theta),2))/4.;

                case 18: return (3*std::sqrt(5/pi)*(1 - 7*std::cos(2*theta))*std::sin(2*phi)*std::sin(2*theta))/8.;

                case 19: return (-3*std::sqrt(5/(2.*pi))*(std::cos(2*theta) + 7*std::cos(4*theta))*std::sin(phi))/8.;

                case 20: return (15*std::cos(theta)*(3 - 7*std::pow(std::cos(theta),2))*std::sin(theta))/(4.*std::sqrt(pi));

                case 21: return (-3*std::sqrt(5/(2.*pi))*std::cos(phi)*(std::cos(2*theta) + 7*std::cos(4*theta)))/8.;

                case 22: return (3*std::sqrt(5/pi)*std::cos(2*phi)*(-2*std::sin(2*theta) + 7*std::sin(4*theta)))/16.;

                case 23: return (-3*std::sqrt(35/(2.*pi))*std::cos(3*phi)*(1 + 2*std::cos(2*theta))*std::pow(std::sin(theta),2))/4.;

                case 24: return (3*std::sqrt(35/pi)*std::cos(4*phi)*std::cos(theta)*std::pow(std::sin(theta),3))/4.;

                default: {
                    TERMINATE_NOT_IMPLEMENTED
                }
            }
            return 0; // make compiler happy
        }

        ///  Compute derivative of real-spherical harmonic with respect to phi angle and divide by sin(theta).
        static double dRlm_dphi_sin_theta(int lm, double theta, double phi)
        {
            switch (lm) {
                case 0: return 0;

                case 1: return -(std::sqrt(3/pi)*std::cos(phi))/2.;

                case 2: return 0;

                case 3: return (std::sqrt(3/pi)*std::sin(phi))/2.;

                case 4: return -(std::sqrt(15/pi)*std::cos(2*phi)*std::sin(theta))/2.;

                case 5: return -(std::sqrt(15/pi)*std::cos(phi)*std::cos(theta))/2.;

                case 6: return 0;

                case 7: return (std::sqrt(15/pi)*std::cos(theta)*std::sin(phi))/2.;

                case 8: return -(std::sqrt(15/pi)*std::cos(phi)*std::sin(phi)*std::sin(theta));

                case 9: return (-3*std::sqrt(35/(2.*pi))*std::cos(3*phi)*std::pow(std::sin(theta),2))/4.;

                case 10: return -(std::sqrt(105/pi)*std::cos(2*phi)*std::sin(2*theta))/4.;

                case 11: return -(std::sqrt(21/(2.*pi))*std::cos(phi)*(3 + 5*std::cos(2*theta)))/8.;

                case 12: return 0;

                case 13: return (std::sqrt(21/(2.*pi))*(3 + 5*std::cos(2*theta))*std::sin(phi))/8.;

                case 14: return -(std::sqrt(105/pi)*std::cos(phi)*std::cos(theta)*std::sin(phi)*std::sin(theta));

                case 15: return (3*std::sqrt(35/(2.*pi))*std::sin(3*phi)*std::pow(std::sin(theta),2))/4.;

                case 16: return (-3*std::sqrt(35/pi)*std::cos(4*phi)*std::pow(std::sin(theta),3))/4.;

                case 17: return (-9*std::sqrt(35/(2.*pi))*std::cos(3*phi)*std::cos(theta)*std::pow(std::sin(theta),2))/4.;

                case 18: return (-3*std::sqrt(5/pi)*std::cos(2*phi)*(3*std::sin(theta) + 7*std::sin(3*theta)))/16.;

                case 19: return (-3*std::sqrt(5/(2.*pi))*std::cos(phi)*(9*std::cos(theta) + 7*std::cos(3*theta)))/16.;

                case 20: return 0;

                case 21: return (3*std::sqrt(5/(2.*pi))*std::cos(theta)*(1 + 7*std::cos(2*theta))*std::sin(phi))/8.;

                case 22: return (-3*std::sqrt(5/pi)*std::sin(2*phi)*(3*std::sin(theta) + 7*std::sin(3*theta)))/16.;

                case 23: return (9*std::sqrt(35/(2.*pi))*std::cos(theta)*std::sin(3*phi)*std::pow(std::sin(theta),2))/4.;

                case 24: return (-3*std::sqrt(35/pi)*std::sin(4*phi)*std::pow(std::sin(theta),3))/4.;

                default: {
                    TERMINATE_NOT_IMPLEMENTED
                }
            }
            return 0; // make compiler happy
        }
};

template <>
inline void SHT::backward_transform<double>(int ld, double const* flm, int nr, int lmmax, double* ftp)
{
    assert(lmmax <= lmmax_);
    assert(ld >= lmmax);
    linalg<CPU>::gemm(1, 0, num_points_, nr, lmmax, &rlm_backward_(0, 0), lmmax_, flm, ld, ftp, num_points_);
}

template <>
inline void SHT::backward_transform<double_complex>(int ld, double_complex const* flm, int nr, int lmmax, double_complex* ftp)
{
    assert(lmmax <= lmmax_);
    assert(ld >= lmmax);
    linalg<CPU>::gemm(1, 0, num_points_, nr, lmmax, &ylm_backward_(0, 0), lmmax_, flm, ld, ftp, num_points_);
}

template <>
inline void SHT::forward_transform<double>(double const* ftp, int nr, int lmmax, int ld, double* flm)
{
    assert(lmmax <= lmmax_);
    assert(ld >= lmmax);
    linalg<CPU>::gemm(1, 0, lmmax, nr, num_points_, &rlm_forward_(0, 0), num_points_, ftp, num_points_, flm, ld);
}

template <>
inline void SHT::forward_transform<double_complex>(double_complex const* ftp, int nr, int lmmax, int ld, double_complex* flm)
{
    assert(lmmax <= lmmax_);
    assert(ld >= lmmax);
    linalg<CPU>::gemm(1, 0, lmmax, nr, num_points_, &ylm_forward_(0, 0), num_points_, ftp, num_points_, flm, ld);
}

inline double SHT::gaunt_ylm(int l1, int l2, int l3, int m1, int m2, int m3)
{
    assert(l1 >= 0);
    assert(l2 >= 0);
    assert(l3 >= 0);
    assert(m1 >= -l1 && m1 <= l1);
    assert(m2 >= -l2 && m2 <= l2);
    assert(m3 >= -l3 && m3 <= l3);

    return std::pow(-1.0, std::abs(m1)) * std::sqrt(double(2 * l1 + 1) * double(2 * l2 + 1) * double(2 * l3 + 1) / fourpi) *
           gsl_sf_coupling_3j(2 * l1, 2 * l2, 2 * l3, 0, 0, 0) *
           gsl_sf_coupling_3j(2 * l1, 2 * l2, 2 * l3, -2 * m1, 2 * m2, 2 * m3);
}

inline double SHT::gaunt_rlm(int l1, int l2, int l3, int m1, int m2, int m3)
{
    assert(l1 >= 0);
    assert(l2 >= 0);
    assert(l3 >= 0);
    assert(m1 >= -l1 && m1 <= l1);
    assert(m2 >= -l2 && m2 <= l2);
    assert(m3 >= -l3 && m3 <= l3);

    double d = 0;
    for (int k1 = -l1; k1 <= l1; k1++)
    {
        for (int k2 = -l2; k2 <= l2; k2++)
        {
            for (int k3 = -l3; k3 <= l3; k3++)
            {
                d += std::real(std::conj(SHT::ylm_dot_rlm(l1, k1, m1)) *
                                         SHT::ylm_dot_rlm(l2, k2, m2) *
                                         SHT::ylm_dot_rlm(l3, k3, m3)) * SHT::gaunt_ylm(l1, l2, l3, k1, k2, k3);
            }
        }
    }
    return d;
}

inline double_complex SHT::gaunt_hybrid(int l1, int l2, int l3, int m1, int m2, int m3)
{
    assert(l1 >= 0);
    assert(l2 >= 0);
    assert(l3 >= 0);
    assert(m1 >= -l1 && m1 <= l1);
    assert(m2 >= -l2 && m2 <= l2);
    assert(m3 >= -l3 && m3 <= l3);

    if (m2 == 0)
    {
        return double_complex(gaunt_ylm(l1, l2, l3, m1, m2, m3), 0.0);
    }
    else
    {
        return (ylm_dot_rlm(l2, m2, m2) * gaunt_ylm(l1, l2, l3, m1, m2, m3) +
                ylm_dot_rlm(l2, -m2, m2) * gaunt_ylm(l1, l2, l3, m1, -m2, m3));
    }
}

inline vector3d<double> SHT::spherical_coordinates(vector3d<double> vc)
{
    vector3d<double> vs;

    double eps = 1e-12;

    vs[0] = vc.length();

    if (vs[0] <= eps)
    {
        vs[1] = 0.0;
        vs[2] = 0.0;
    }
    else
    {
        vs[1] = std::acos(vc[2] / vs[0]); // theta = cos^{-1}(z/r)

        if (std::abs(vc[0]) > eps || std::abs(vc[1]) > eps)
        {
            vs[2] = std::atan2(vc[1], vc[0]); // phi = tan^{-1}(y/x)
            if (vs[2] < 0.0) vs[2] += twopi;
        }
        else
        {
            vs[2] = 0.0;
        }
    }

    return vs;
}

inline void SHT::spherical_harmonics(int lmax, double theta, double phi, double_complex* ylm)
{
    double x = std::cos(theta);
    std::vector<double> result_array(lmax + 1);

    for (int l = 0; l <= lmax; l++)
    {
        for (int m = 0; m <= l; m++)
        {
            double_complex z = std::exp(double_complex(0.0, m * phi));
            ylm[Utils::lm_by_l_m(l, m)] = gsl_sf_legendre_sphPlm(l, m, x) * z;
            if (m % 2)
            {
                ylm[Utils::lm_by_l_m(l, -m)] = -std::conj(ylm[Utils::lm_by_l_m(l, m)]);
            }
            else
            {
                ylm[Utils::lm_by_l_m(l, -m)] = std::conj(ylm[Utils::lm_by_l_m(l, m)]);
            }
        }
    }
}

inline void SHT::spherical_harmonics(int lmax, double theta, double phi, double* rlm)
{
    int lmmax = (lmax + 1) * (lmax + 1);
    std::vector<double_complex> ylm(lmmax);
    spherical_harmonics(lmax, theta, phi, &ylm[0]);

    double t = std::sqrt(2.0);

    rlm[0] = y00;

    for (int l = 1; l <= lmax; l++)
    {
        for (int m = -l; m < 0; m++)
            rlm[Utils::lm_by_l_m(l, m)] = t * ylm[Utils::lm_by_l_m(l, m)].imag();

        rlm[Utils::lm_by_l_m(l, 0)] = ylm[Utils::lm_by_l_m(l, 0)].real();

        for (int m = 1; m <= l; m++)
            rlm[Utils::lm_by_l_m(l, m)] = t * ylm[Utils::lm_by_l_m(l, m)].real();
    }
}

inline void SHT::uniform_coverage()
{
    tp_(0, 0) = pi;
    tp_(1, 0) = 0;

    for (int k = 1; k < num_points_ - 1; k++)
    {
        double hk = -1.0 + double(2 * k) / double(num_points_ - 1);
        tp_(0, k) = acos(hk);
        double t = tp_(1, k - 1) + 3.80925122745582 / sqrt(double(num_points_)) / sqrt(1 - hk * hk);
        tp_(1, k) = fmod(t, twopi);
    }

    tp_(0, num_points_ - 1) = 0;
    tp_(1, num_points_ - 1) = 0;
}

inline void SHT::convert(int lmax__, double const* f_rlm__, double_complex* f_ylm__)
{
    int lm = 0;
    for (int l = 0; l <= lmax__; l++)
    {
        for (int m = -l; m <= l; m++)
        {
            if (m == 0)
            {
                f_ylm__[lm] = f_rlm__[lm];
            }
            else
            {
                int lm1 = Utils::lm_by_l_m(l, -m);
                f_ylm__[lm] = ylm_dot_rlm(l, m, m) * f_rlm__[lm] + ylm_dot_rlm(l, m, -m) * f_rlm__[lm1];
            }
            lm++;
        }
    }
}

inline void SHT::convert(int lmax__, double_complex const* f_ylm__, double* f_rlm__)
{
    int lm = 0;
    for (int l = 0; l <= lmax__; l++)
    {
        for (int m = -l; m <= l; m++)
        {
            if (m == 0)
            {
                f_rlm__[lm] = std::real(f_ylm__[lm]);
            }
            else
            {
                int lm1 = Utils::lm_by_l_m(l, -m);
                f_rlm__[lm] = std::real(rlm_dot_ylm(l, m, m) * f_ylm__[lm] + rlm_dot_ylm(l, m, -m) * f_ylm__[lm1]);
            }
            lm++;
        }
    }
}

};

#endif // __SHT_H__

max_stars_repo_name: cocteautwins/SIRIUS-develop

max_issues_repo_name: cocteautwins/SIRIUS-develop

max_forks_repo_name: cocteautwins/SIRIUS-develop