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

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

#include "ccl.h"

/*
 * We use the remapping procedure in https://arxiv.org/pdf/1601.07230.pdf
 * to map w0-wa models onto wa == 0 models.
 * The functions below implement this procedure.
*/

// I snuck a private function from background.c into here. :P
// The default routines create a full spline but we just need one value.
// This may seem like over optimizing, but in testing initializing halofit
// is an order of magnitude or more slower if we don't do this.
void compute_chi(double a, ccl_cosmology *cosmo, double * chi, int * stat);

static double zdrag_eh(ccl_parameters *params) {
  // eqn 4 of Eisenstein & Hu 1998
  double OMh2 = (params->Omega_c + params->Omega_b) * params->h * params->h;
  double OBh2 = params->Omega_b * params->h * params->h;
  double b1 = 0.313 * pow(OMh2, -0.419) * (1 + 0.607*pow(OMh2, 0.674));
  double b2 = 0.238 * pow(OMh2, 0.223);
  return 1291 * pow(OMh2, 0.251) * (1 + b1*pow(OBh2, b2)) / (1 + 0.659*pow(OMh2, 0.828));
}

struct hf_model_match_data {
  double chi_drag;
  double a;
  ccl_cosmology *cosmo;
  int *status;
};

static ccl_cosmology *create_w0eff_cosmo(double w0eff, ccl_cosmology *cosmo, int *status) {
  // create a cosmology with the same parameters as the input except w0-wa. Instead
  // the cosmology is created with w0 = w0eff.
  ccl_parameters params_w0eff;
  double norm_pk;
  double mnu[3];
  int i;

  for(i=0; i<3; ++i)
    mnu[i] = 0;
  for(i=0; i<cosmo->params.N_nu_mass; ++i)
    mnu[i] = cosmo->params.m_nu[i];

  if (isnan(cosmo->params.A_s))
    norm_pk = cosmo->params.sigma8;
  else
    norm_pk = cosmo->params.A_s;

  params_w0eff = ccl_parameters_create(
    cosmo->params.Omega_c, cosmo->params.Omega_b, cosmo->params.Omega_k,
    cosmo->params.Neff, mnu, cosmo->params.N_nu_mass,
    w0eff, 0, cosmo->params.h, norm_pk,
    cosmo->params.n_s, cosmo->params.bcm_log10Mc, cosmo->params.bcm_etab,
    cosmo->params.bcm_ks, cosmo->params.mu_0, cosmo->params.sigma_0,
    cosmo->params.c1_mg, cosmo->params.c2_mg, cosmo->params.lambda_mg,
    cosmo->params.nz_mgrowth,
    cosmo->params.z_mgrowth, cosmo->params.df_mgrowth, status);

    if(*status != 0)
      return NULL;

    return ccl_cosmology_create(params_w0eff, cosmo->config);
}

static double w0eff_func(double w0eff, void *p) {
  // function used to compare the distance to the CMB in a test cosmology to
  // to the value in the original cosmology
  // returns chi_eff - chi
  struct hf_model_match_data *hfd = (struct hf_model_match_data*)p;
  ccl_cosmology *cosmo_w0eff = NULL;
  double chi_drag_w0eff, tmp, zdrag_w0eff;

  // make the equivalent cosmology
  cosmo_w0eff = create_w0eff_cosmo(w0eff, hfd->cosmo, hfd->status);
  if (cosmo_w0eff == NULL) {
    *(hfd->status) = CCL_ERROR_MEMORY;
    return NAN;
  }

  // get the comoving distance to zdrag
  zdrag_w0eff = zdrag_eh(&(cosmo_w0eff->params));
  compute_chi(1.0 / (1.0 + zdrag_w0eff), cosmo_w0eff, &tmp, hfd->status);
  chi_drag_w0eff = tmp;
  compute_chi(hfd->a, cosmo_w0eff, &tmp, hfd->status);
  chi_drag_w0eff -= tmp;
  if (*(hfd->status) != 0) {
    return NAN;
  }

  ccl_parameters_free(&(cosmo_w0eff->params));
  ccl_cosmology_free(cosmo_w0eff);
  return chi_drag_w0eff - hfd->chi_drag;
}

static double get_w0eff(double a, struct hf_model_match_data data) {
  // For a given input w0-wa cosmology, this function solves for the value of
  // w0eff such that the comoving distance from a to the CMB in a cosmology
  // with the same parameters, but with w0, wa = w0eff, 0, is the same as the
  // original cosmology.
  double w0eff, w0eff_low = -2.0, w0eff_high = -0.35;
  double flow, fhigh;
  int itr, max_itr = 1000, gsl_status;
  const gsl_root_fsolver_type *T;
  gsl_root_fsolver *s;
  gsl_function F;

  data.a = a;
  data.chi_drag = ccl_comoving_radial_distance(data.cosmo, 1.0 / (1.0 + zdrag_eh(&(data.cosmo->params))), data.status);
  data.chi_drag -= ccl_comoving_radial_distance(data.cosmo, a, data.status);
  if(*(data.status) != 0) {
    ccl_cosmology_set_status_message(
      data.cosmo,
      "ccl_halofit.c: get_w0eff(): "
      "could not compute chi_drag for cosmology\n");
    return NAN;
  }

  F.function = &w0eff_func;
  F.params = &data;

  // we have to bound the root, otherwise return -1
  // we will fiil in any -1's in the calling routine
  flow = w0eff_func(w0eff_low, &data);
  fhigh = w0eff_func(w0eff_high, &data);
  if (flow * fhigh > 0) {
    return -1;
  }

  T = gsl_root_fsolver_brent;
  s = gsl_root_fsolver_alloc(T);
  if (s == NULL) {
    *(data.status) = CCL_ERROR_MEMORY;
  }
  else {
    gsl_root_fsolver_set(s, &F, w0eff_low, w0eff_high);

    itr = 0;
    do {
      itr++;
      gsl_status = gsl_root_fsolver_iterate(s);
      if (gsl_status == GSL_EBADFUNC)
        break;

      w0eff = gsl_root_fsolver_root(s);
      w0eff_low = gsl_root_fsolver_x_lower(s);
      w0eff_high = gsl_root_fsolver_x_upper(s);

      gsl_status = gsl_root_test_interval(
        w0eff_low, w0eff_high,
        1e-6,
        1e-6);
    } while (gsl_status == GSL_CONTINUE && itr < max_itr);

    gsl_root_fsolver_free(s);

    if (gsl_status != GSL_SUCCESS || itr >= max_itr) {
      ccl_raise_gsl_warning(
        gsl_status, "ccl_halofit.c: get_w0eff(): error in root finding for the halofit matching cosmology\n");
      *(data.status) |= gsl_status;
    }
  }

  return w0eff;
}

/* helper data and functions for integrals

 the integral is

 \int dlnk \Delta^{2}(k) \exp(-k^{2}R^{2})

 for halofit, we also need the first and second derivatives of this
 integral
*/
struct hf_int_data {
  double r;
  double r2;
  double a;
  ccl_cosmology *cosmo;
  ccl_f2d_t *plin;
  int *status;
  gsl_integration_cquad_workspace *workspace;
};

static double gauss_norm_int_func(double lnk, void *p) {
  struct hf_int_data *hfd = (struct hf_int_data*)p;
  double k = exp(lnk);
  double k2 = k*k;

  return (
    ccl_f2d_t_eval(hfd->plin, lnk, hfd->a, hfd->cosmo, hfd->status) *
    k*k2/2.0/M_PI/M_PI *
    exp(-k2 * (hfd->r2)));
}

static double onederiv_gauss_norm_int_func(double lnk, void *p) {
  struct hf_int_data *hfd = (struct hf_int_data*)p;
  double k = exp(lnk);
  double k2 = k*k;

  return (
    ccl_f2d_t_eval(hfd->plin, lnk, hfd->a, hfd->cosmo, hfd->status) *
    k*k2/2.0/M_PI/M_PI *
    exp(-k2 * (hfd->r2)) *
    (-k2 * 2.0 * (hfd->r)));
}

static double twoderiv_gauss_norm_int_func(double lnk, void *p) {
  struct hf_int_data *hfd = (struct hf_int_data*)p;
  double k = exp(lnk);
  double k2 = k*k;

  return (
    ccl_f2d_t_eval(hfd->plin, lnk, hfd->a, hfd->cosmo, hfd->status) *
    k*k2/2.0/M_PI/M_PI *
    exp(-k2 * (hfd->r2)) *
    (-2.0*k2 + 4.0*k2*k2 * (hfd->r2)));
}

// function whose root is \sigma^2{rsigma, a} = 1
static double rsigma_func(double rsigma, void *p) {
  struct hf_int_data *hfd = (struct hf_int_data*)p;
  double result, lnkmin, lnkmax;
  gsl_function F;
  int gsl_status;

  lnkmin = hfd->plin->lkmin;
  lnkmax = hfd->plin->lkmax;
  hfd->r = rsigma;
  hfd->r2 = rsigma * rsigma;
  F.function = &gauss_norm_int_func;
  F.params = (void *)hfd;
  gsl_status = gsl_integration_cquad(
    &F, lnkmin, lnkmax,
    0.0, hfd->cosmo->gsl_params.INTEGRATION_SIGMAR_EPSREL,
    hfd->workspace, &result, NULL, NULL);

  if (gsl_status != GSL_SUCCESS) {
    ccl_raise_gsl_warning(
      gsl_status,
      "ccl_halofit.c: rsigma_func(): error in integration "
      "for finding the halofit non-linear scale\n");
    *(hfd->status) |= gsl_status;
  }

  return result - 1.0;
}

static double get_rsigma(double a, struct hf_int_data data) {
  double rsigma, rlow = 1e-2, rhigh = 1e2;
  double flow, fhigh;
  int itr, max_itr = 1000, gsl_status;
  const gsl_root_fsolver_type *T;
  gsl_root_fsolver *s;
  gsl_function F;

  data.a = a;
  F.function = &rsigma_func;
  F.params = &data;

  // we have to bound the root, otherwise return -1
  // we will fiil in any -1's in the calling routine
  flow = rsigma_func(rlow, &data);
  fhigh = rsigma_func(rhigh, &data);
  if (flow * fhigh > 0) {
    return -1;
  }

  T = gsl_root_fsolver_brent;
  s = gsl_root_fsolver_alloc(T);
  if (s == NULL) {
    *(data.status) = CCL_ERROR_MEMORY;
  }
  else {
    gsl_root_fsolver_set(s, &F, rlow, rhigh);

    itr = 0;
    do {
      itr++;
      gsl_status = gsl_root_fsolver_iterate(s);
      if (gsl_status == GSL_EBADFUNC)
        break;

      rsigma = gsl_root_fsolver_root(s);
      rlow = gsl_root_fsolver_x_lower(s);
      rhigh = gsl_root_fsolver_x_upper(s);

      gsl_status = gsl_root_test_interval(
        rlow, rhigh,
        data.cosmo->gsl_params.INTEGRATION_SIGMAR_EPSREL,
        data.cosmo->gsl_params.INTEGRATION_SIGMAR_EPSREL);
    } while (gsl_status == GSL_CONTINUE && itr < max_itr);

    gsl_root_fsolver_free(s);

    if (gsl_status != GSL_SUCCESS || itr >= max_itr) {
      ccl_raise_gsl_warning(
        gsl_status, "ccl_halofit.c: get_rsigma(): error in root finding for the halofit non-linear scale\n");
      *(data.status) |= gsl_status;
    }
  }

  return rsigma;
}

/*
 * Allocate a new struct for storing halofit data
 * @param cosmo Cosmological data
 * @param int, status of computations
 */
halofit_struct* ccl_halofit_struct_new(ccl_cosmology *cosmo,
                                       ccl_f2d_t *plin, int *status) {
  size_t n_a;
  int i, gsl_status;
  double lnkmin, lnkmax;
  double *a_vec = NULL;
  double *vals = NULL;
  double *vals_om = NULL;
  double *vals_de = NULL;
  halofit_struct *hf = NULL;
  struct hf_int_data data;
  gsl_function F;
  gsl_integration_cquad_workspace *workspace = NULL;
  double result;
  double sigma2, rsigma, dsigma2drsigma;
  struct hf_model_match_data data_w0eff;
  ccl_cosmology *cosmo_w0eff = NULL;

  // compute spline point locations and integral bounds
  // note that the spline point locations in `a` determine a radius by
  // solving sigma2(R, a) = 1
  // it is this radius that is needed for the subsequent splines of
  // the derivatives.
  lnkmin = plin->lkmin;
  lnkmax = plin->lkmax;
  if(plin->fa != NULL) {
    n_a = plin->fa->size;
    a_vec = plin->fa->x;
  }
  else if(plin->fka != NULL) {
    n_a = plin->fka->interp_object.ysize;
    a_vec = plin->fka->yarr;
  }
  else {
    *status = CCL_ERROR_SPLINE;
    ccl_cosmology_set_status_message(cosmo,
           "ccl_halofit.c: ccl_halofit_struct_new(): "
           "input pk2d has no splines.\n");
  }

  if(*status == 0) {
    ///////////////////////////////////////////////////////
    // memory allocation
    hf = (halofit_struct*)malloc(sizeof(halofit_struct));
    if (hf == NULL) {
      *status = CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(
                                       cosmo,
                                       "ccl_halofit.c: ccl_halofit_struct_new(): "
                                       "memory could not be allocated for halofit_struct\n");
    } else {
      hf->rsigma = NULL;
      hf->sigma2 = NULL;
      hf->n_eff = NULL;
      hf->C = NULL;
      hf->weff = NULL;
      hf->omeff = NULL;
      hf->deeff = NULL;
    }
  }

  if (*status == 0) {
    workspace = gsl_integration_cquad_workspace_alloc(cosmo->gsl_params.N_ITERATION);
    if (workspace == NULL) {
      *status = CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_halofit.c: ccl_halofit_struct_new(): "
        "memory could not be allocated for cquad workspace\n");
    }
  }

  if (*status == 0) {
    vals = (double*)malloc(sizeof(double) * n_a);
    if (vals == NULL) {
      *status = CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_halofit.c: ccl_halofit_struct_new(): "
        "memory could not be allocated for results vector\n");
    }
  }

  if (*status == 0) {
    vals_om = (double*)malloc(sizeof(double) * n_a);
    if (vals_om == NULL) {
      *status = CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_halofit.c: ccl_halofit_struct_new(): "
        "memory could not be allocated for OmegaM results vector\n");
    }
  }

  if (*status == 0) {
    vals_de = (double*)malloc(sizeof(double) * n_a);
    if (vals_de == NULL) {
      *status = CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_halofit.c: ccl_halofit_struct_new(): "
        "memory could not be allocated for OmegaDE results vector\n");
    }
  }

  ////////////////////////////////////////////////////////
  // if wa != 0, then we need to find an equivalent
  // cosmology with wa = 0
  data_w0eff.cosmo = cosmo;
  data_w0eff.status = status;

  if (*status == 0) {
    if (cosmo->params.wa != 0) {
      for(i=0; i<n_a; ++i) {
        vals[i] = get_w0eff(a_vec[i], data_w0eff);
        if (*status != 0) {
          *status = CCL_ERROR_ROOT;
          ccl_cosmology_set_status_message(
            cosmo,
            "ccl_halofit.c: ccl_halofit_struct_new(): "
            "could not solve for effective value of w0 for w0-wa cosmology\n");
          break;
        }

        // now get omeff and deff
        cosmo_w0eff = create_w0eff_cosmo(vals[i], cosmo, status);
        if (cosmo_w0eff == NULL) {
          *status = CCL_ERROR_MEMORY;
          ccl_cosmology_set_status_message(
            cosmo,
            "ccl_halofit.c: ccl_halofit_struct_new(): "
            "could not allocat memory for effective w0 for w0-wa cosmology\n");
          break;
        }

        vals_om[i] = ccl_omega_x(cosmo_w0eff, a_vec[i], ccl_species_m_label, status) +
          ccl_omega_x(cosmo_w0eff, a_vec[i], ccl_species_nu_label, status);
        vals_de[i] = ccl_omega_x(cosmo_w0eff, a_vec[i], ccl_species_l_label, status);

        ccl_parameters_free(&(cosmo_w0eff->params));
        ccl_cosmology_free(cosmo_w0eff);
        cosmo_w0eff = NULL;

        if (*status != 0) {
          ccl_cosmology_set_status_message(
            cosmo,
            "ccl_halofit.c: ccl_halofit_struct_new(): "
            "could not compute OmegaM and OmegaDE for cosmology\n");
          break;
        }
      }
    } else {
      for(i=0; i<n_a; ++i) {
        vals[i] = cosmo->params.w0;
        vals_om[i] = ccl_omega_x(cosmo, a_vec[i], ccl_species_m_label, status) +
          ccl_omega_x(cosmo, a_vec[i], ccl_species_nu_label, status);
        vals_de[i] = ccl_omega_x(cosmo, a_vec[i], ccl_species_l_label, status);
        if (*status != 0) {
          ccl_cosmology_set_status_message(
            cosmo,
            "ccl_halofit.c: ccl_halofit_struct_new(): "
            "could not compute OmegaM and OmegaDE for cosmology\n");
          break;
        }
      }
    }
  }

  // spline the weff values
  if (*status == 0) {
    hf->weff = gsl_spline_alloc(gsl_interp_akima, n_a);
    if (hf->weff == NULL) {
      *status = CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_halofit.c: ccl_halofit_struct_new(): "
        "memory could not be allocated for weff spline\n");
    }
  }

  if (*status == 0) {
    gsl_status = gsl_spline_init(hf->weff, a_vec, vals, n_a);
    if (gsl_status != GSL_SUCCESS) {
      *status = CCL_ERROR_SPLINE;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_halofit.c: ccl_halofit_struct_new(): could not build weff spline\n");
    }
  }

  // spline the omeff values
  if (*status == 0) {
    hf->omeff = gsl_spline_alloc(gsl_interp_akima, n_a);
    if (hf->omeff == NULL) {
      *status = CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_halofit.c: ccl_halofit_struct_new(): "
        "memory could not be allocated for omeff spline\n");
    }
  }

  if (*status == 0) {
    gsl_status = gsl_spline_init(hf->omeff, a_vec, vals_om, n_a);
    if (gsl_status != GSL_SUCCESS) {
      *status = CCL_ERROR_SPLINE;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_halofit.c: ccl_halofit_struct_new(): could not build omeff spline\n");
    }
  }

  // spline the deeff values
  if (*status == 0) {
    hf->deeff = gsl_spline_alloc(gsl_interp_akima, n_a);
    if (hf->deeff == NULL) {
      *status = CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_halofit.c: ccl_halofit_struct_new(): "
        "memory could not be allocated for deeff spline\n");
    }
  }

  if (*status == 0) {
    gsl_status = gsl_spline_init(hf->deeff, a_vec, vals_de, n_a);
    if (gsl_status != GSL_SUCCESS) {
      *status = CCL_ERROR_SPLINE;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_halofit.c: ccl_halofit_struct_new(): could not build deeff spline\n");
    }
  }

  ///////////////////////////////////////////////////////
  // find the nonlinear scale at each scale factor
  if (*status == 0) {
    // setup for integrations
    data.status = status;
    data.cosmo = cosmo;
    data.plin = plin;
    data.workspace = workspace;

    for (i=0; i<n_a; ++i) {
      vals[i] = get_rsigma(a_vec[i], data);
      if (*status != 0) {
        *status = CCL_ERROR_ROOT;
        ccl_cosmology_set_status_message(
          cosmo,
          "ccl_halofit.c: ccl_halofit_struct_new(): "
          "could not solve for non-linear scale for halofit at scale factor %f\n", a_vec[i]);
        break;
      }
    }

    // now go backwards and fill any -1's
    // one must work, so set an error if not
    if (vals[n_a-1] == -1) {
      *status = CCL_ERROR_ROOT;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_halofit.c: ccl_halofit_struct_new(): "
        "could not solve for non-linear scale for halofit at least once\n");
    }
    if (*status == 0) {
      // linearly set rsigma to a very small number at high-z
      double min_a = -1, max_a = -1;
      double max_val = -1;
      double w;

      // first non -1 value and scale factor of last -1 value
      for (i=1; i<n_a; ++i) {
        if (vals[i] != -1) {
          max_a = a_vec[i];
          max_val = vals[i];
          break;
        }
      }
      // scale factor of first -1 value
      for (i=0; i<n_a; ++i) {
        if (vals[i] == -1) {
          min_a = a_vec[i];
          break;
        }
      }

      if (min_a != -1) {
        // at least one value is -1 so set the zeroth value
        vals[0] = 1e-6;

        // interp any values that remain -1
        for(i=1; i<n_a-1; ++i) {
          if (vals[i] == -1) {
            w = (a_vec[i] - min_a) / (max_a - min_a);
            vals[i] = w * max_val + (1.0 - w) * vals[0];
          }
        }
      }
    }
  }

  // spline the non-linear scales
  if (*status == 0) {
    hf->rsigma = gsl_spline_alloc(gsl_interp_akima, n_a);
    if (hf->rsigma == NULL) {
      *status = CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_halofit.c: ccl_halofit_struct_new(): "
        "memory could not be allocated for Rsigma spline\n");
    }
  }

  if (*status == 0) {
    gsl_status = gsl_spline_init(hf->rsigma, a_vec, vals, n_a);
    if (gsl_status != GSL_SUCCESS) {
      *status = CCL_ERROR_SPLINE;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_halofit.c: ccl_halofit_struct_new(): could not build Rsigma spline\n");
    }
  }

  ///////////////////////////////////////////////////////
  // now compute sigma2(R) at each spline point
  // and spline that
  // this should be close to 1 OFC, but better to use the exact value
  if (*status == 0) {
    for (i=0; i<n_a; ++i) {
      data.a = a_vec[i];
      vals[i] = rsigma_func(gsl_spline_eval(hf->rsigma, a_vec[i], NULL), &data) + 1;
      if (*status != 0) {
        ccl_cosmology_set_status_message(
          cosmo,
          "ccl_halofit.c: ccl_halofit_struct_new(): could not eval "
          "points for sigma2(R) spline\n");
        break;
      }
    }
  }

  if (*status == 0) {
    hf->sigma2 = gsl_spline_alloc(gsl_interp_akima, n_a);
    if (hf->sigma2 == NULL) {
      *status = CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_halofit.c: ccl_halofit_struct_new(): "
        "memory could not be allocated for sigma2(R) spline\n");
    }
  }

  if (*status == 0) {
    gsl_status = gsl_spline_init(hf->sigma2, a_vec, vals, n_a);
    if (gsl_status != GSL_SUCCESS) {
      *status = CCL_ERROR_SPLINE;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_halofit.c: ccl_halofit_struct_new(): could not build sigma2(R) spline\n");
    }
  }

  ///////////////////////////////////////////////////////
  // now compute the effective spectral index
  if (*status == 0) {
    F.function = &onederiv_gauss_norm_int_func;
    F.params = &data;

    for (i=0; i<n_a; ++i) {
      rsigma = gsl_spline_eval(hf->rsigma, a_vec[i], NULL);
      sigma2 = gsl_spline_eval(hf->sigma2, a_vec[i], NULL);

      data.a = a_vec[i];
      data.r = rsigma;
      data.r2 = rsigma * rsigma;

      gsl_status = gsl_integration_cquad(
        &F,
        lnkmin, lnkmax,
        0.0, cosmo->gsl_params.INTEGRATION_SIGMAR_EPSREL,
        workspace, &result, NULL, NULL);

      if (gsl_status != GSL_SUCCESS) {
        *status = CCL_ERROR_INTEG;
        ccl_raise_gsl_warning(
          gsl_status,
          "ccl_power.c: ccl_halofit_struct_new(): could not eval "
          "points for n_eff spline\n");
        ccl_cosmology_set_status_message(
          cosmo,
          "ccl_halofit.c: ccl_halofit_struct_new(): could not eval "
          "points for n_eff spline\n");
        break;
      }

      // this is n_eff but expressed in terms of linear derivs
      // see eqn A5 of Takahashi et al.
      vals[i] = -rsigma / sigma2 * result - 3.0;
    }
  }

  if (*status == 0) {
    hf->n_eff = gsl_spline_alloc(gsl_interp_akima, n_a);
    if (hf->n_eff == NULL) {
      *status = CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_halofit.c: ccl_halofit_struct_new(): "
        "memory could not be allocated for n_eff spline\n");
    }
  }

  if (*status == 0) {
    gsl_status = gsl_spline_init(hf->n_eff, a_vec, vals, n_a);
    if (gsl_status != GSL_SUCCESS) {
      *status = CCL_ERROR_SPLINE;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_halofit.c: ccl_halofit_struct_new(): could not build n_eff spline\n");
    }
  }

  ///////////////////////////////////////////////////////
  // now compute the curvature C
  if (*status == 0) {
    F.function = &twoderiv_gauss_norm_int_func;
    F.params = &data;

    for (i=0; i<n_a; ++i) {
      rsigma = gsl_spline_eval(hf->rsigma, a_vec[i], NULL);
      sigma2 = gsl_spline_eval(hf->sigma2, a_vec[i], NULL);

      // we need to solve for the deriv we need from n_eff here
      dsigma2drsigma = gsl_spline_eval(hf->n_eff, a_vec[i], NULL);
      dsigma2drsigma = (dsigma2drsigma + 3.0) / (-rsigma / sigma2);

      data.a = a_vec[i];
      data.r = rsigma;
      data.r2 = rsigma * rsigma;

      gsl_status = gsl_integration_cquad(
        &F,
        lnkmin, lnkmax,
        0.0, cosmo->gsl_params.INTEGRATION_SIGMAR_EPSREL,
        workspace, &result, NULL, NULL);

      if (gsl_status != GSL_SUCCESS) {
        *status = CCL_ERROR_INTEG;
        ccl_raise_gsl_warning(
          gsl_status,
          "ccl_power.c: ccl_halofit_struct_new(): could not eval "
          "points for C spline\n");
        ccl_cosmology_set_status_message(
          cosmo,
          "ccl_halofit.c: ccl_halofit_struct_new(): could not eval "
          "points for C spline\n");
        break;
      }

      // this is C but expressed in terms of linear derivs
      // see eqn A5 of Takahashi et al.
      vals[i] = (
        -1.0 * (
          result * rsigma * rsigma / sigma2 +
          dsigma2drsigma * rsigma / sigma2 -
          dsigma2drsigma * dsigma2drsigma * rsigma * rsigma / sigma2 / sigma2));
    }
  }

  if (*status == 0) {
    hf->C = gsl_spline_alloc(gsl_interp_akima, n_a);
    if (hf->C == NULL) {
      *status = CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_halofit.c: ccl_halofit_struct_new(): "
        "memory could not be allocated for C spline\n");
    }
  }

  if (*status == 0) {
    gsl_status = gsl_spline_init(hf->C, a_vec, vals, n_a);
    if (gsl_status != GSL_SUCCESS) {
      *status = CCL_ERROR_SPLINE;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_halofit.c: ccl_halofit_struct_new(): could not build C spline\n");
    }
  }

  // free stuff on the way out
  // if the status is non-zero, then we should free any data already
  // accumulated
  if (*status != 0) {
    ccl_halofit_struct_free(hf);
    hf = NULL;
  }
  gsl_integration_cquad_workspace_free(workspace);
  free(vals);
  free(vals_om);
  free(vals_de);
  if (cosmo_w0eff != NULL) {
    ccl_parameters_free(&(cosmo_w0eff->params));
    ccl_cosmology_free(cosmo_w0eff);
  }

  return hf;
}

/*
 * Free a halofit struct
 * @param hf, pointer to halofit struct to free
 */
void ccl_halofit_struct_free(halofit_struct *hf) {
  if (hf != NULL) {
    if (hf->rsigma != NULL)
      gsl_spline_free(hf->rsigma);

    if (hf->sigma2 != NULL)
      gsl_spline_free(hf->sigma2);

    if (hf->n_eff != NULL)
      gsl_spline_free(hf->n_eff);

    if (hf->C != NULL)
      gsl_spline_free(hf->C);

    if (hf->weff != NULL)
      gsl_spline_free(hf->weff);

    if (hf->omeff != NULL)
      gsl_spline_free(hf->omeff);

    if (hf->deeff != NULL)
      gsl_spline_free(hf->deeff);

    free(hf);
  }
}

/**
 * Computes the halofit non-linear power spectrum
 * @param cosmo: cosmology object containing parameters
 * @param lk: natural logarithm of wavenumber in units of Mpc^{-1}
 * @param a: scale factor normalised to a=1 today
 * @param status: Status flag: 0 if there are no errors, non-zero otherwise
 * @param hf: halofit splines for evaluating the power spectrum
 * @return halofit_matter_power: halofit power spectrum, P(k), units of Mpc^{3}
 */
double ccl_halofit_power(ccl_cosmology *cosmo, ccl_f2d_t *plin,
                         double lk, double a, halofit_struct *hf, int *status) {
  double rsigma, neff, C;
  double ksigma, weffa, omegaMz, omegaDEwz, kh;
  double PkL, PkNL, f1, f2, f3, an, bn, cn, gamman, alphan, betan, nun, mun, y, fy;
  double DeltakL, DeltakL_tilde, DeltakQ, DeltakH, DeltakHprime, DeltakNL;
  double Qnu, fnu;
  double f1a, f2a, f3a, f1b, f2b, f3b, fb_frac;
  double neff2, neff3, neff4;
  double kh2, y2;
  double delta2_norm, om_nu;
  double k=exp(lk);

  // all eqns are from Takahashi et al. unless stated otherwise
  // eqns A4 - A5
  rsigma = gsl_spline_eval(hf->rsigma, a, NULL);
  neff = gsl_spline_eval(hf->n_eff, a, NULL);
  C = gsl_spline_eval(hf->C, a, NULL);

  weffa = cosmo->params.w0;
  omegaMz = ccl_omega_x(cosmo, a, ccl_species_m_label, status);
  omegaDEwz = ccl_omega_x(cosmo, a, ccl_species_l_label, status);

  // not using these to match CLASS better - might be a bug in CLASS
  // weffa = gsl_spline_eval(hf->weff, a, NULL);
  // omegaMz = gsl_spline_eval(hf->omeff, a, NULL);
  // omegaDEwz = gsl_spline_eval(hf->deeff, a, NULL);

  ksigma = 1.0 / rsigma;
  neff2 = neff * neff;
  neff3 = neff2 * neff;
  neff4 = neff3 * neff;

  delta2_norm = k*k*k/2.0/M_PI/M_PI;

  // compute the present day neutrino massive neutrino fraction
  // uses all neutrinos even if they are moving fast
  om_nu = cosmo->params.sum_nu_masses / 93.14 / cosmo->params.h / cosmo->params.h;
  fnu = om_nu / (cosmo->params.Omega_m);

  // eqns A6 - A13 of Takahashi et al.
  an = pow(
    10.0,
    1.5222 + 2.8553*neff + 2.3706*neff2 + 0.9903*neff3 +
    0.2250*neff4 - 0.6038*C + 0.1749*omegaDEwz*(1.0 + weffa));
  bn = pow(10.0, -0.5642 + 0.5864*neff + 0.5716*neff2 - 1.5474*C + 0.2279*omegaDEwz*(1.0 + weffa));
  cn = pow(10.0, 0.3698 + 2.0404*neff + 0.8161*neff2 + 0.5869*C);
  gamman = 0.1971 - 0.0843*neff + 0.8460*C;
  alphan = fabs(6.0835 + 1.3373*neff - 0.1959*neff2 - 5.5274*C);
  betan = 2.0379 - 0.7354*neff + 0.3157*neff2 + 1.2490*neff3 + 0.3980*neff4 - 0.1682*C;
  mun = 0.0;
  nun = pow(10.0, 5.2105 + 3.6902*neff);

  // eqns C17 and C18 for Smith et al.
  if (fabs(1.0 - omegaMz) > 0.01) {
    f1a = pow(omegaMz, -0.0732);
    f2a = pow(omegaMz, -0.1423);
    f3a = pow(omegaMz, 0.0725);
    f1b = pow(omegaMz, -0.0307);
    f2b = pow(omegaMz, -0.0585);
    f3b = pow(omegaMz, 0.0743);
    fb_frac = omegaDEwz / (1.0 - omegaMz);
    f1 = fb_frac * f1b + (1.0 - fb_frac) * f1a;
    f2 = fb_frac * f2b + (1.0 - fb_frac) * f2a;
    f3 = fb_frac * f3b + (1.0 - fb_frac) * f3a;
  } else {
    f1 = 1.0;
    f2 = 1.0;
    f3 = 1.0;
  }

  // correction to betan from Bird et al., eqn A10
  betan += (fnu * (1.081 + 0.395*neff2));

  // eqns A1 - A3
  PkL = ccl_f2d_t_eval(plin, lk, a, cosmo, status);
  y = k/ksigma;
  y2 = y * y;
  fy = y/4.0 + y2/8.0;
  DeltakL = PkL * delta2_norm;

  // correction to DeltakL from Bird et al., eqn A9
  kh = k / cosmo->params.h;
  kh2 = kh * kh;
  DeltakL_tilde = DeltakL * (1.0 + fnu * (47.48 * kh2) / (1.0 + 1.5 * kh2));
  DeltakQ = DeltakL * pow(1.0 + DeltakL_tilde, betan) / (1.0 + alphan*DeltakL_tilde) * exp(-fy);

  DeltakHprime = an * pow(y, 3.0*f1) / (1.0 + bn*pow(y, f2) + pow(cn*f3*y, 3.0 - gamman));
  DeltakH = DeltakHprime / (1.0 + mun/y + nun/y2);

  // correction to DeltakH from Bird et al., eqn A6-A7
  Qnu = fnu * (0.977 - 18.015 * (cosmo->params.Omega_m - 0.3));
  DeltakH *= (1.0 + Qnu);

  DeltakNL = DeltakQ + DeltakH;
  PkNL = DeltakNL / delta2_norm;

  // we check the status once
  if(*status != 0)
    return NAN;
  else
    return PkNL;
}

max_stars_repo_name: Jappenn/CCL

max_issues_repo_name: Jappenn/CCL

max_forks_repo_name: Jappenn/CCL


text:
#ifndef SurrogateHeatFluidDriver_h
#define SurrogateHeatFluidDriver_h

#include <memory>
#include <vector>

#include "Teuchos_ParameterList.hpp"
#include "Teuchos_RCP.hpp"

#include <gsl/gsl>

#include "Assembly_Model.h"
#include "Single_Pin_Subchannel.h"
#include "Single_Pin_Conduction.h"

namespace enrico {

//! Driver for coupled subchannel/heat conduction simulations. This driver
//! decouples the fluid and solid phase in an explicit manner (i.e. no
//! iteration between the convective heat transfer linking the phases within
//! a single solve), which is sufficient given that this project targets
//! steady-state simulations via fixed point iteration.
class SurrogateHeatFluidDriver {
public:
  //@{
  //! Typedefs
  using SP_Assembly = std::shared_ptr<Assembly_Model>;
  using RCP_PL = Teuchos::RCP<Teuchos::ParameterList>;
  //@}

  //! solve the thermal-hydraulic problem for the fluid and solid phases
  //! for a given power distribution
  void solve(const std::vector<double>& powers);

  //! fluid temperature
  virtual std::vector<double> fluid_temperature() const { return d_pin_temps; }

  //! fluid density
  virtual std::vector<double> density() const { return d_pin_densities; }

  //! solid temperature
  virtual std::vector<double> solid_temperature() const { return d_solid_temps; }

  double pressure_bc_; //! System pressure in [MPa]

private:
  // >>> DATA
  SP_Assembly d_assembly;

  // Note that for channel-centered equations, Nx and Ny are one greater
  //  than number of pins in each direction.
  int d_Nx, d_Ny, d_Nz;

  // Cross sectional area of each channel
  std::vector<double> d_areas;

  // Mass flow rate (kg/s) in each channel
  std::vector<double> d_mdots;

  //! subchannel solver for a single channel
  std::unique_ptr<Single_Pin_Subchannel> d_pin_subchannel;

  //! heat conduction solver for a single rod
  std::unique_ptr<Single_Pin_Conduction> d_pin_conduction;

  //! coolant temperature in [K] for each channel, of total length given by the
  //! product of the number of pins by the number of axial cells
  std::vector<double> d_pin_temps;

  //! coolant density in [g/cm^3] for each channel, of total length given by the
  //! product of the number of pins by the number of axial cells
  std::vector<double> d_pin_densities;

  //! solid temperature in [K]
  std::vector<double> d_solid_temps;

public:
  // Constructor
  SurrogateHeatFluidDriver(SP_Assembly assembly,
                       RCP_PL subchannel_params,
                       RCP_PL conduction_params,
                       const std::vector<double>& dz);

  //! heat source
  std::vector<double> d_pin_powers;

private:
  //! Set up the sizes of solution arrays
  void generate_arrays();

  int channel_index(int ix, int iy) const
  {
    Expects(ix < d_Nx);
    Expects(iy < d_Ny);
    return ix + d_Nx * iy;
  }

  // Solve the fluid flow equations
  void solve_fluid(const std::vector<double>& power);

  // Solve the solid equations
  void solve_heat(const std::vector<double>& power,
    const std::vector<double>& channel_temp, std::vector<double>& fuel_temp);
};

//---------------------------------------------------------------------------//
} // end namespace enrico

//---------------------------------------------------------------------------//
#endif // SurrogateHeatFluidDriver_h

//---------------------------------------------------------------------------//
// end of SurrogateHeatFluidDriver.h
//---------------------------------------------------------------------------//

ext: h

max_stars_repo_name: sphamil/enrico

max_issues_repo_name: sphamil/enrico

max_forks_repo_name: sphamil/enrico


text:
#include <stdio.h>
#include <gsl/gsl_spline.h>

void spline_interp(long data_size, long out_size, \
        double *data_x, double *data_y, \
        double *out_x, double *out_y) {

    // initialize
    gsl_interp_accel *acc = gsl_interp_accel_alloc();
    gsl_spline *spline    = gsl_spline_alloc(gsl_interp_cspline, data_size);
    gsl_spline_init(spline, data_x, data_y, data_size);

    // evaluate
    int ii;
    for (ii=0; ii < out_size; ii++) {
        out_y[ii] = gsl_spline_eval (spline, out_x[ii], acc);
    }

    // free memory
    gsl_spline_free(spline);
    gsl_interp_accel_free(acc);
}

ext: c

max_stars_repo_name: thompsonphys/gwsurrogate


text:
#ifndef __Q_INCS
#define __Q_INCS
#include <alloca.h>
#include <assert.h>
#include <ctype.h>
#include <errno.h>
#include <fcntl.h>
#include <float.h>
// #include <gsl/gsl_vector.h>
// #include <gsl/gsl_matrix.h>
// #include <gsl/gsl_blas.h>
// #include <gsl/gsl_linalg.h>
#include <inttypes.h>
#include <limits.h>
// TODO P4 do not think this is needed:#include <malloc.h>
#include <math.h>
#include <memory.h>
#include <omp.h>
#include <stdbool.h>
#include <stddef.h>
#include <stdlib.h>
#include <stdint.h>
#include <stdio.h>
#include <string.h>
#include <sys/mman.h>
#include <sys/stat.h>
#ifndef MAC_OSX
// TODO P1 Is this needed?#include <sys/sysinfo.h>
#endif
#include <sys/time.h>
#include <time.h>
#include <unistd.h>
#include "q_macros.h"
#endif

max_stars_repo_name: subramon/qlu

max_issues_repo_name: subramon/qlu

max_forks_repo_name: subramon/qlu


text:
#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_math.h>

#include "allvars.h"
#include "proto.h"

/*! \file snIa_heating.c
 *  \brief routines for heating by type Ia supernovae
 */

#if defined(SNIA_HEATING)

static struct snheatingdata_in
{
  MyDouble Pos[3];
  MyFloat Density;
  MyFloat Energy;
  MyFloat Hsml;
  int NodeList[NODELISTLENGTH];
}
 *SnheatingdataIn, *SnheatingdataGet;

static double hubble_a, ascale, a3inv;

void snIa_heating(void)
{
  int i, j, k;
  int ndone_flag, ndone;
  int ngrp, sendTask, recvTask, place, nexport, nimport, dummy;
  double dt;
  MPI_Status status;

  if(ThisTask == 0)
    {
      printf("Beginning snIa heating\n");
      fflush(stdout);
    }

  CPU_Step[CPU_MISC] += measure_time();

  if(All.ComovingIntegrationOn)
    {
      ascale = All.Time;
      a3inv = 1.0 / (ascale * ascale * ascale);
      hubble_a = hubble_function(All.Time);
    }
  else
    hubble_a = ascale = a3inv = 1;

  /* allocate buffers to arrange communication */

  Ngblist = (int *) mymalloc(NumPart * sizeof(int));

  All.BunchSize =
    (int) ((All.BufferSize * 1024 * 1024) / (sizeof(struct data_index) + sizeof(struct data_nodelist) +
					     2 * sizeof(struct snheatingdata_in)));
  DataIndexTable = (struct data_index *) mymalloc(All.BunchSize * sizeof(struct data_index));
  DataNodeList = (struct data_nodelist *) mymalloc(All.BunchSize * sizeof(struct data_nodelist));

  /** Let's first spread the feedback energy */

  i = FirstActiveParticle;	/* first particle for this task */

  do
    {
      for(j = 0; j < NTask; j++)
	{
	  Send_count[j] = 0;
	  Exportflag[j] = -1;
	}

      /* do local particles and prepare export list */

      for(nexport = 0; i >= 0; i = NextActiveParticle[i])
	if(P[i].Type == 4)
	  if(snIaheating_evaluate(i, 0, &nexport, Send_count) < 0)
	    break;

#ifdef MYSORT
      mysort_dataindex(DataIndexTable, nexport, sizeof(struct data_index), data_index_compare);
#else
      qsort(DataIndexTable, nexport, sizeof(struct data_index), data_index_compare);
#endif

      MPI_Allgather(Send_count, NTask, MPI_INT, Sendcount_matrix, NTask, MPI_INT, MPI_COMM_WORLD);

      for(j = 0, nimport = 0, Recv_offset[0] = 0, Send_offset[0] = 0; j < NTask; j++)
	{
	  Recv_count[j] = Sendcount_matrix[j * NTask + ThisTask];
	  nimport += Recv_count[j];

	  if(j > 0)
	    {
	      Send_offset[j] = Send_offset[j - 1] + Send_count[j - 1];
	      Recv_offset[j] = Recv_offset[j - 1] + Recv_count[j - 1];
	    }
	}

      SnheatingdataGet = (struct snheatingdata_in *) mymalloc(nimport * sizeof(struct snheatingdata_in));
      SnheatingdataIn = (struct snheatingdata_in *) mymalloc(nexport * sizeof(struct snheatingdata_in));

      for(j = 0; j < nexport; j++)
	{
	  place = DataIndexTable[j].Index;

	  for(k = 0; k < 3; k++)
	    SnheatingdataIn[j].Pos[k] = P[place].Pos[k];

	  SnheatingdataIn[j].Hsml = PPP[place].Hsml;
	  SnheatingdataIn[j].Density = P[place].DensAroundStar;

	  dt = (P[place].TimeBin ? (1 << P[place].TimeBin) : 0) * All.Timebase_interval / hubble_a;

	  SnheatingdataIn[j].Energy = All.SnIaHeatingRate * dt * P[place].Mass;

	  memcpy(SnheatingdataIn[j].NodeList,
		 DataNodeList[DataIndexTable[j].IndexGet].NodeList, NODELISTLENGTH * sizeof(int));
	}

      /* exchange particle data */
      for(ngrp = 1; ngrp < (1 << PTask); ngrp++)
	{
	  sendTask = ThisTask;
	  recvTask = ThisTask ^ ngrp;

	  if(recvTask < NTask)
	    {
	      if(Send_count[recvTask] > 0 || Recv_count[recvTask] > 0)
		{
		  /* get the particles */
		  MPI_Sendrecv(&SnheatingdataIn[Send_offset[recvTask]],
			       Send_count[recvTask] * sizeof(struct snheatingdata_in), MPI_BYTE,
			       recvTask, TAG_DENS_A,
			       &SnheatingdataGet[Recv_offset[recvTask]],
			       Recv_count[recvTask] * sizeof(struct snheatingdata_in), MPI_BYTE,
			       recvTask, TAG_DENS_A, MPI_COMM_WORLD, &status);
		}
	    }
	}

      myfree(SnheatingdataIn);

      /* now do the particles that were sent to us */

      for(j = 0; j < nimport; j++)
	snIaheating_evaluate(j, 1, &dummy, &dummy);

      if(i < 0)
	ndone_flag = 1;
      else
	ndone_flag = 0;

      MPI_Allreduce(&ndone_flag, &ndone, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);

      myfree(SnheatingdataGet);
    }
  while(ndone < NTask);

  myfree(DataNodeList);
  myfree(DataIndexTable);
  myfree(Ngblist);

  CPU_Step[CPU_COOLINGSFR] += measure_time();
}

int snIaheating_evaluate(int target, int mode, int *nexport, int *nsend_local)
{
  int startnode, numngb, j, n, listindex = 0;
  MyFloat *pos, h_i, dt, rho;
  double dx, dy, dz, h_i2, r2, r, u, hinv, hinv3, wk, energy, egy;

  if(mode == 0)
    {
      pos = P[target].Pos;
      rho = P[target].DensAroundStar;

      dt = (P[target].TimeBin ? (1 << P[target].TimeBin) : 0) * All.Timebase_interval / hubble_a;

      energy = All.SnIaHeatingRate * dt * P[target].Mass;

      h_i = PPP[target].Hsml;
    }
  else
    {
      pos = SnheatingdataGet[target].Pos;
      rho = SnheatingdataGet[target].Density;
      energy = SnheatingdataGet[target].Energy;
      h_i = SnheatingdataGet[target].Hsml;
    }

  /* initialize variables before SPH loop is started */
  h_i2 = h_i * h_i;

  /* Now start the actual SPH computation for this particle */
  if(mode == 0)
    {
      startnode = All.MaxPart;	/* root node */
    }
  else
    {
      startnode = SnheatingdataGet[target].NodeList[0];
      startnode = Nodes[startnode].u.d.nextnode;	/* open it */
    }

  while(startnode >= 0)
    {
      while(startnode >= 0)
	{
	  numngb = ngb_treefind_variable(pos, h_i, target, &startnode, mode, nexport, nsend_local);

	  if(numngb < 0)
	    return -1;

	  for(n = 0; n < numngb; n++)
	    {
	      j = Ngblist[n];

	      if(P[j].Type == 0 && P[j].Mass > 0)
		{
		  dx = pos[0] - P[j].Pos[0];
		  dy = pos[1] - P[j].Pos[1];
		  dz = pos[2] - P[j].Pos[2];
#ifdef PERIODIC			/*  now find the closest image in the given box size  */
		  if(dx > boxHalf_X)
		    dx -= boxSize_X;
		  if(dx < -boxHalf_X)
		    dx += boxSize_X;
		  if(dy > boxHalf_Y)
		    dy -= boxSize_Y;
		  if(dy < -boxHalf_Y)
		    dy += boxSize_Y;
		  if(dz > boxHalf_Z)
		    dz -= boxSize_Z;
		  if(dz < -boxHalf_Z)
		    dz += boxSize_Z;
#endif
		  r2 = dx * dx + dy * dy + dz * dz;

		  if(r2 < h_i2)
		    {
		      r = sqrt(r2);
		      hinv = 1 / h_i;
#ifndef  TWODIMS
		      hinv3 = hinv * hinv * hinv;
#else
		      hinv3 = hinv * hinv / boxSize_Z;
#endif

		      u = r * hinv;

		      if(u < 0.5)
			wk = hinv3 * (KERNEL_COEFF_1 + KERNEL_COEFF_2 * (u - 1) * u * u);
		      else
			wk = hinv3 * KERNEL_COEFF_5 * (1.0 - u) * (1.0 - u) * (1.0 - u);

		      egy = P[j].Mass * wk / rho * energy;

		      SphP[j].Entropy += egy / P[j].Mass * GAMMA_MINUS1 /
			pow(SphP[j].d.Density * a3inv, GAMMA_MINUS1);

		    }
		}
	    }
	}

      if(mode == 1)
	{
	  listindex++;
	  if(listindex < NODELISTLENGTH)
	    {
	      startnode = SnheatingdataGet[target].NodeList[listindex];
	      if(startnode >= 0)
		startnode = Nodes[startnode].u.d.nextnode;	/* open it */
	    }
	}
    }

  return 0;
}

#endif

max_stars_repo_name: egpbos/egp

max_issues_repo_name: egpbos/egp

max_forks_repo_name: egpbos/egp


text:
#include "q_incs.h"
#include <lapacke.h>
#include "eigenvectors.h"

int
eigenvectors(
             uint64_t n,
             double *W,
             double *A,
             double **X
            )
{
  int status = 0;

  char jobz = 'V'; /*want eigenvectors, use 'N' for eigenvalues only*/
  char uplo = 'U'; /*'U' for upper triangle of X, L for lower triangle of X*/
  int N = n;
  int LDA = N; /* dimensions of X = LDA by N*/

  for (uint64_t i = 0; i < n; i++ ) {
    for (uint64_t j = 0; j < n; j++ ) {
      A[i * n + j] = X[i][j];
    }
  }

  status = LAPACKE_dsyev(LAPACK_COL_MAJOR, jobz, uplo, N, A, LDA, W);
  return status;
}


text:
/* statistics/gsl_statistics_double.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Jim Davies, 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 2 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.
 */

#ifndef __GSL_STATISTICS_DOUBLE_H__
#define __GSL_STATISTICS_DOUBLE_H__

#include <stddef.h>
#include <gsl/gsl_types.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

GSL_EXPORT double gsl_stats_mean (const double data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_variance (const double data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_sd (const double data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_variance_with_fixed_mean (const double data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_sd_with_fixed_mean (const double data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_absdev (const double data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_skew (const double data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_kurtosis (const double data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_lag1_autocorrelation (const double data[], const size_t stride, const size_t n);

GSL_EXPORT double gsl_stats_covariance (const double data1[], const size_t stride1,const double data2[], const size_t stride2, const size_t n);

GSL_EXPORT double gsl_stats_variance_m (const double data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_sd_m (const double data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_absdev_m (const double data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_skew_m_sd (const double data[], const size_t stride, const size_t n, const double mean, const double sd);
GSL_EXPORT double gsl_stats_kurtosis_m_sd (const double data[], const size_t stride, const size_t n, const double mean, const double sd);
GSL_EXPORT double gsl_stats_lag1_autocorrelation_m (const double data[], const size_t stride, const size_t n, const double mean);

GSL_EXPORT double gsl_stats_covariance_m (const double data1[], const size_t stride1,const double data2[], const size_t stride2, const size_t n, const double mean1, const double mean2);

/* DEFINED FOR FLOATING POINT TYPES ONLY */

GSL_EXPORT double gsl_stats_wmean (const double w[], const size_t wstride, const double data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_wvariance (const double w[], const size_t wstride, const double data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_wsd (const double w[], const size_t wstride, const double data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_wvariance_with_fixed_mean (const double w[], const size_t wstride, const double data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_wsd_with_fixed_mean (const double w[], const size_t wstride, const double data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_wabsdev (const double w[], const size_t wstride, const double data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_wskew (const double w[], const size_t wstride, const double data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_wkurtosis (const double w[], const size_t wstride, const double data[], const size_t stride, const size_t n);

GSL_EXPORT double gsl_stats_wvariance_m (const double w[], const size_t wstride, const double data[], const size_t stride, const size_t n, const double wmean);
GSL_EXPORT double gsl_stats_wsd_m (const double w[], const size_t wstride, const double data[], const size_t stride, const size_t n, const double wmean);
GSL_EXPORT double gsl_stats_wabsdev_m (const double w[], const size_t wstride, const double data[], const size_t stride, const size_t n, const double wmean);
GSL_EXPORT double gsl_stats_wskew_m_sd (const double w[], const size_t wstride, const double data[], const size_t stride, const size_t n, const double wmean, const double wsd);
GSL_EXPORT double gsl_stats_wkurtosis_m_sd (const double w[], const size_t wstride, const double data[], const size_t stride, const size_t n, const double wmean, const double wsd);

/* END OF FLOATING POINT TYPES */

GSL_EXPORT double gsl_stats_pvariance (const double data1[], const size_t stride1, const size_t n1, const double data2[], const size_t stride2, const size_t n2);
GSL_EXPORT double gsl_stats_ttest (const double data1[], const size_t stride1, const size_t n1, const double data2[], const size_t stride2, const size_t n2);

GSL_EXPORT double gsl_stats_max (const double data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_min (const double data[], const size_t stride, const size_t n);
GSL_EXPORT void gsl_stats_minmax (double * min, double * max, const double data[], const size_t stride, const size_t n);

GSL_EXPORT size_t gsl_stats_max_index (const double data[], const size_t stride, const size_t n);
GSL_EXPORT size_t gsl_stats_min_index (const double data[], const size_t stride, const size_t n);
GSL_EXPORT void gsl_stats_minmax_index (size_t * min_index, size_t * max_index, const double data[], const size_t stride, const size_t n);

GSL_EXPORT double gsl_stats_median_from_sorted_data (const double sorted_data[], const size_t stride, const size_t n) ;
GSL_EXPORT double gsl_stats_quantile_from_sorted_data (const double sorted_data[], const size_t stride, const size_t n, const double f) ;

__END_DECLS

#endif /* __GSL_STATISTICS_DOUBLE_H__ */

max_stars_repo_name: dynaryu/vaws

max_issues_repo_name: dynaryu/vaws

max_forks_repo_name: dynaryu/vaws


text:
#ifndef __INTEGRAL_H__
#define __INTEGRAL_H__
#include "integral.h"
#endif

#ifndef __GSLEXTRA_H__
#define __GSLEXTRA_H__
#include "../include/gslextra.h"
#endif

#ifndef __GSLPRINT_H__
#define __GSLPRINT_H__
#include "../include/gslprint.h"
#endif

#include <gsl/gsl_eigen.h>

// calculate the attraction energy matrix element (Z integral)
double nuclear_attraction_energy_matrix_element(orbital * a, orbital * b, atomic_orbital * atom_HEAD);

// calculate the single electron hamiltonian matrix (core hamiltonian matrix) element
double single_electron_hamiltonian_matrix_element(orbital * a, orbital * b, atomic_orbital * atom_HEAD);

// calculate the fock matrix element
double fock_matrix_element(gsl_quad_tensor * v, gsl_matrix * density_matrix, gsl_matrix * h_matrix, int i, int j, int length);

// calculate the Hartree-Fock energy of the system
double HF_energy(gsl_quad_tensor * v, gsl_matrix * density_matrix, gsl_matrix * h_matrix, int length);

// calculate the kinetic energy matrix
void kinetic_energy_matrix(gsl_matrix * dest, orbital * HEAD, int length);

// calculate the attraction energy matrix (Z integrals)
void nuclear_attraction_energy_matrix(gsl_matrix * dest, orbital * HEAD, atomic_orbital * atom_HEAD, int length);

// calculate the single electron hamiltonian matrix (core hamiltonian matrix)
void core_hamiltonian_matrix(gsl_matrix * dest, orbital * HEAD, atomic_orbital * atom_HEAD, int length);

//obtain the quad tensor of the two-electron Coulomb integrals
void two_electron_quad_tensor(gsl_quad_tensor * dest, orbital * HEAD, int length);

//obtain the fock matrix
void fock_matrix(gsl_matrix * dest, gsl_quad_tensor * v, gsl_matrix * density_matrix, gsl_matrix * h_matrix, int length);

// perform initial guess of the coefficient matrix by performing diagonalization of core hamiltonian matrix
void initial_guess(gsl_matrix * dest, gsl_matrix * core_hamiltonian, gsl_matrix * S, int length);

// calculate density matrix from coefficient matrix and the number of electrons
void density_matrix(gsl_matrix * dest, gsl_matrix * coef, int el_num, int length);

//perform RHF as well as printing information
int RHF_SCF_print(double * tot_energy, gsl_vector * energy, gsl_matrix * coef, orbital * HEAD, atomic_orbital * atom_HEAD, int length, int el_num, int iteration_max, double errmax, int countmax, double alpha, int mixing_type, int SCF_INITIAL_FLAG, int SCF_FOCK_FLAG, int SCF_COEF_FLAG, int FOCK_FLAG);

//calculate the nuclei repulsion energy
double nuclei_repulsion(atomic_orbital * atomlist_HEAD);

max_stars_repo_name: Walter-Feng/Hartree-Fock

max_issues_repo_name: Walter-Feng/Hartree-Fock

max_forks_repo_name: Walter-Feng/Hartree-Fock


text:
#pragma once

#include "exp_curve_fitter.h"
#include "optimizer_job.h"
#include "random_utils.h"
#include "schedule_params.h"
#include <gsl/gsl-lite.hpp>
#include <vector>

namespace angonoka::stun {
/**
    How many iterations without improvement before
    considering optimization complete.
*/
enum class MaxIdleIters : std::int_fast32_t;

/**
    Optimization algorithm based on stochastic tunneling.

    This is the primary facade for doing stochastic tunneling
    optimization.
*/
class Optimizer {
public:
    /**
        Constructor.

        @param params           Scheduling parameters
        @param batch_size       Number of iterations per update
        @param max_idle_iters   Stopping condition
    */
    Optimizer(
        const ScheduleParams& params,
        BatchSize batch_size,
        MaxIdleIters max_idle_iters);

    /**
        Run stochastic tunneling optimization batch.

        Does batch_size number of iterations and adjusts the estimated
        progress accordingly.
    */
    void update() noexcept;

    /**
        Checks if the stopping condition has been met.

        @return True when further improvements are unlikely
    */
    [[nodiscard]] bool has_converged() const noexcept;

    /**
        Estimated optimization progress from 0.0 to 1.0.

        @return Progress from 0.0 to 1.0
    */
    [[nodiscard]] float estimated_progress() const noexcept;

    /**
        The best schedule so far.

        @return A schedule.
    */
    [[nodiscard]] Schedule schedule() const;

    /**
        The best makespan so far.

        @return Makespan.
    */
    [[nodiscard]] float normalized_makespan() const;

    /**
        Reset the optimization to initial state.
    */
    void reset();

    /**
        Get the current ScheduleParams object.

        @return Schedule parameters.
    */
    [[nodiscard]] const ScheduleParams& params() const;

    /**
        Set the ScheduleParams object.

        @param params ScheduleParams object
    */
    void params(const ScheduleParams& params);

private:
    struct Impl;

    int16 batch_size;
    int16 max_idle_iters;
    int16 idle_iters{0};
    int16 epochs{0};
    float last_progress{0.F};
    float last_makespan{0.F};
    ExpCurveFitter exp_curve;

    /**
        An optimization job and a PRNG.

        @var random_utils   Random number generator utilities
        @var job            Optimization job
    */
    struct Job {
        /**
            Constructor.

            @param params       Scheduling parameters
            @param batch_size   Number of iterations per update
        */
        Job(const ScheduleParams& params, BatchSize batch_size);
        RandomUtils random_utils;
        OptimizerJob job;
    };
    std::vector<Job> jobs;
};
} // namespace angonoka::stun

max_stars_repo_name: coffee-lord/angonoka

max_issues_repo_name: coffee-lord/angonoka

max_forks_repo_name: coffee-lord/angonoka


text:
/* multimin/diff.c
 *
 * Copyright (C) 2000 David Morrison
 *
 * 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_multimin.h>

int
gsl_multimin_diff (const gsl_multimin_function * f,
                   const gsl_vector * x, gsl_vector * g)
{
  size_t i, n = f->n;

  double h = GSL_SQRT_DBL_EPSILON;

  gsl_vector * x1 = gsl_vector_alloc (n);  /* FIXME: pass as argument */

  gsl_vector_memcpy (x1, x);

  for (i = 0; i < n; i++)
    {
      double fl, fh;

      double xi = gsl_vector_get (x, i);
      double dx = fabs(xi) * h;

      if (dx == 0.0) dx = h;

      gsl_vector_set (x1, i, xi + dx);
      fh = GSL_MULTIMIN_FN_EVAL(f, x1);

      gsl_vector_set (x1, i, xi - dx);
      fl = GSL_MULTIMIN_FN_EVAL(f, x1);

      gsl_vector_set (x1, i, xi);
      gsl_vector_set (g, i, (fh - fl) / (2.0 * dx));
    }

  gsl_vector_free (x1);

  return GSL_SUCCESS;
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#include <stdio.h>

//#ifndef DARWIN
//#include <malloc.h>
//#endif

#include <stddef.h>
#include <Python.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>
#include <string.h>
#include <math.h>

/*
void inplace(double* npyArray3D, int npyLength1D, int npylength2D,
              int npylength2D, double* invec, int n)
{
  int i;
  for (i=0; i<n; i++) {
    invec[i] = npyArray3D[i];
  }
}
*/

int print_matrix(FILE *f, const gsl_matrix *m)
{
  int status, n = 0;

  for (size_t i = 0; i < m->size1; i++) {
    for (size_t j = 0; j < m->size2; j++) {
      if ((status = fprintf(f,"%g ", gsl_matrix_get(m, i, j))) < 0)
        return -1;
      n += status;
    }

    if ((status = fprintf(f, "\n")) < 0)
      return -1;
    n += status;
  }
  return n;
}

int print_vector(FILE *f, const gsl_vector *m)
{
  int status, n = 0;

  for (size_t i = 0; i < m->size; i++) {
    if ((status = fprintf(f, "%g ", gsl_vector_get(m, i))) < 0)
      return -1;
    n += status;
  }
  return n;
}

void print_mat(double* mat, int dim1, int dim2)
{
  int i, j;
  printf("\n");
  for (i=0; i<dim1; i++) {
    for (j=0; j<dim2; j++) {
      printf("%6.2f,", mat[i*dim1 + j]);
    }
    printf("\n");
  }
}
void print_vec(double* vec, int dim1)
{
  int i;
  printf("\n");
  for (i=0; i<dim1; i++) {
    printf("%6.2f,", vec[i]);
  }
  printf("\n");
}

int sum(int* npyArray3D, int npyLength1D, int npyLength2D, int npyLength3D)
{
  int i, j, k;
  int sum = 0;

  for (i=0; i<npyLength1D; i++)
    for (j=0; j<npyLength2D; j++)
      for (k=0; k<npyLength3D; k++)
        sum += npyArray3D[i*npyLength3D*npyLength2D + k*npyLength2D + j];

  return sum;
}

double get_det(PyObject *A)
{
  int MAT_DIM = 6;
  int i, signum;
  double det;
  int nInts = PyList_Size(A);

  gsl_matrix *m = gsl_matrix_alloc(MAT_DIM, MAT_DIM);
  gsl_permutation *p;

  p = gsl_permutation_alloc(m->size1);

  for (i=0; i<nInts; i++)
  {
    PyObject *oo = PyList_GetItem(A, i);
    //printf("%6.2f\n",PyFloat_AS_DOUBLE(oo));
    gsl_matrix_set (m, i%MAT_DIM, i/MAT_DIM, PyFloat_AS_DOUBLE(oo));
  }

  gsl_linalg_LU_decomp(m, p, &signum);
  det = gsl_linalg_LU_det(m, signum);
  //printf("%6.2f\n",det);
  return det;
}

double get_overlap(double* gr_icov, int gr_dim1, int gr_dim2,
                   double* gr_mn, int gr_mn_dim, double gr_icov_det,
                   double* st_icov, int st_dim1, int st_dim2,
                   double* st_mn, int st_mn_dim, double st_icov_det)
{
  int MAT_DIM = gr_dim1;
  int i, j, signum;
  double ApB_det, d_temp, result;
  gsl_permutation *p;

  gsl_matrix *A        = gsl_matrix_alloc(MAT_DIM, MAT_DIM);
  gsl_matrix *B        = gsl_matrix_alloc(MAT_DIM, MAT_DIM);
  gsl_matrix *ApB      = gsl_matrix_alloc(MAT_DIM, MAT_DIM);
  gsl_vector *a        = gsl_vector_alloc(MAT_DIM);
  gsl_vector *b        = gsl_vector_alloc(MAT_DIM);
  gsl_vector *AapBb    = gsl_vector_alloc(MAT_DIM);
  gsl_vector *c        = gsl_vector_alloc(MAT_DIM);
  gsl_vector *v_temp   = gsl_vector_alloc(MAT_DIM);
  gsl_vector *v_temp2  = gsl_vector_alloc(MAT_DIM);
  gsl_vector *amc      = gsl_vector_alloc(MAT_DIM); //will hold a - c
  gsl_vector *bmc      = gsl_vector_alloc(MAT_DIM); //will hold b - c

  p = gsl_permutation_alloc(A->size1);

//for (l=0; l<1; l++){ //DELETE

  //Inserting values into matricies and vectors
  for (i=0; i<MAT_DIM; i++)
  {
    for (j=0; j<MAT_DIM; j++)
    {
      gsl_matrix_set (A, i, j, gr_icov[i*MAT_DIM + j]);
      gsl_matrix_set (B, i, j, st_icov[i*MAT_DIM + j]);
    }
  }

  for (i=0; i<MAT_DIM; i++)
  {
    gsl_vector_set (a, i, gr_mn[i]);
    gsl_vector_set (b, i, st_mn[i]);
  }

  // Adding A and B together and storing in ApB
  gsl_matrix_set_zero(ApB);
  gsl_matrix_add(ApB, A);
  gsl_matrix_add(ApB, B);

  // Storing the result A*a + B*b in AapBb
  gsl_vector_set_zero(AapBb);
  gsl_blas_dsymv(CblasUpper, 1.0, A, a, 1.0, AapBb);
  gsl_blas_dsymv(CblasUpper, 1.0, B, b, 1.0, AapBb);

  // Getting determinant of ApB
  gsl_linalg_LU_decomp(ApB, p, &signum);
  //ApB_det = gsl_linalg_LU_det(ApB, signum);
  ApB_det = fabs(gsl_linalg_LU_det(ApB, signum)); //temp doctoring determinant

  // Solve for c
  gsl_linalg_LU_solve(ApB, p, AapBb, c);

  // Compute the overlap formula
  gsl_vector_set_zero(v_temp);
  gsl_blas_dcopy(a, v_temp);       //v_temp holds a
  gsl_blas_daxpy(-1.0, c, v_temp); //v_temp holds a - c
  gsl_blas_dcopy(v_temp, amc);     //amc holds a - c

  // CAN'T HAVE v_temp and v_temp2 be the same vector.
  // Results in 0's being stored in v_temp2.
  gsl_blas_dsymv(CblasUpper, 1.0, A, v_temp, 0.0, v_temp2);
  //v_temp2 holds A (a-c)

  result = 0.0;
  gsl_blas_ddot(v_temp2, amc, &d_temp); //d_temp holds (a-c)^T A (a-c)

  result += d_temp;

  gsl_vector_set_zero(v_temp);
  gsl_blas_dcopy(b, v_temp);       //v_temp holds b
  gsl_blas_daxpy(-1.0, c, v_temp); //v_temp holds b - c
  gsl_blas_dcopy(v_temp, bmc);     //bmc holds b - c

  // CAN'T HAVE v_temp and v_temp2 be the same vector.
  // Results in 0's being stored in v_temp2.
  gsl_blas_dsymv(CblasUpper, 1.0, B, v_temp, 0.0, v_temp2);
  //v_temp2 holds A (b-c)

  gsl_blas_ddot(v_temp2, bmc, &d_temp); //d_temp holds (b-c)^T B (b-c)
  result += d_temp;

  result = -0.5 * result;
  result = exp(result);

  result *= sqrt((gr_icov_det * st_icov_det/ApB_det)
                / pow(2*M_PI, MAT_DIM));
//} //DELETE THIS

  // Freeing memory
  gsl_matrix_free(A);
  gsl_matrix_free(B);
  gsl_matrix_free(ApB);
  gsl_vector_free(a);
  gsl_vector_free(b);
  gsl_vector_free(AapBb);
  gsl_vector_free(c);
  gsl_vector_free(v_temp);
  gsl_vector_free(v_temp2);
  gsl_vector_free(amc);
  gsl_vector_free(bmc);

  gsl_permutation_free(p);

  return result;
}

double get_overlap2(PyObject *gr_icov, PyObject *gr_mn, double gr_icov_det,
                    PyObject *st_icov, PyObject *st_mn, double st_icov_det)
{
  int MAT_DIM = 6;
  int i, j, signum;
  double ApB_det, d_temp, result;
  PyObject *o1, *o2;
  gsl_permutation *p;

  gsl_matrix *A        = gsl_matrix_alloc(MAT_DIM, MAT_DIM);
  gsl_matrix *B        = gsl_matrix_alloc(MAT_DIM, MAT_DIM);
  gsl_matrix *ApB      = gsl_matrix_alloc(MAT_DIM, MAT_DIM);
  gsl_vector *a        = gsl_vector_alloc(MAT_DIM);
  gsl_vector *b        = gsl_vector_alloc(MAT_DIM);
  gsl_vector *AapBb    = gsl_vector_alloc(MAT_DIM);
  gsl_vector *c        = gsl_vector_alloc(MAT_DIM);
  gsl_vector *v_temp   = gsl_vector_alloc(MAT_DIM);
  gsl_vector *v_temp2  = gsl_vector_alloc(MAT_DIM);
  gsl_vector *amc      = gsl_vector_alloc(MAT_DIM); //will hold a - c
  gsl_vector *bmc      = gsl_vector_alloc(MAT_DIM); //will hold b - c

  p = gsl_permutation_alloc(A->size1);

  //Inserting values into matricies and vectors
  for (i=0; i<MAT_DIM; i++)
  {
    for (j=0; j<MAT_DIM; j++)
    {
      o1 = PyList_GetItem(gr_icov, i*MAT_DIM + j);
      gsl_matrix_set (A, i, j, PyFloat_AS_DOUBLE(o1));
      o2 = PyList_GetItem(st_icov, i*MAT_DIM + j);
      gsl_matrix_set (B, i, j, PyFloat_AS_DOUBLE(o2));
    }
  }

  for (i=0; i<MAT_DIM; i++)
  {
    o1 = PyList_GetItem(gr_mn, i);
    gsl_vector_set (a, i, PyFloat_AS_DOUBLE(o1));
    o2 = PyList_GetItem(st_mn, i);
    gsl_vector_set (b, i, PyFloat_AS_DOUBLE(o2));
  }

  // Adding A and B together and storing in ApB
  gsl_matrix_set_zero(ApB);
  gsl_matrix_add(ApB, A);
  gsl_matrix_add(ApB, B);

  // Storing the result A*a + B*b in AapBb
  gsl_vector_set_zero(AapBb);
  gsl_blas_dgemv(CblasNoTrans,
                 1.0, A, a,
                 1.0, AapBb);
  gsl_blas_dgemv(CblasNoTrans,
                 1.0, B, b,
                 1.0, AapBb);

  // Getting determinant of ApB
  gsl_linalg_LU_decomp(ApB, p, &signum);
  //ApB_det = gsl_linalg_LU_det(ApB, signum);
  ApB_det = fabs(gsl_linalg_LU_det(ApB, signum)); //temp doctoring determinant

  // Solve for c
  gsl_linalg_LU_solve(ApB, p, AapBb, c);

  // Compute the overlap formula
  gsl_vector_set_zero(v_temp);
  gsl_blas_dcopy(a, v_temp);       //v_temp holds a
  gsl_blas_daxpy(-1.0, c, v_temp); //v_temp holds a - c
  gsl_blas_dcopy(v_temp, amc);     //amc holds a - c

  // CAN'T HAVE v_temp and v_temp2 be the same vector.
  // Results in 0's being stored in v_temp2.
  gsl_blas_dgemv(CblasNoTrans, 1.0, A, v_temp, 0.0, v_temp2);
  //v_temp2 holds A (a-c)

  result = 0.0;
  gsl_blas_ddot(v_temp2, amc, &d_temp); //d_temp holds (a-c)^T A (a-c)

  result += d_temp;

  gsl_vector_set_zero(v_temp);
  gsl_blas_dcopy(b, v_temp);       //v_temp holds b
  gsl_blas_daxpy(-1.0, c, v_temp); //v_temp holds b - c
  gsl_blas_dcopy(v_temp, bmc);     //bmc holds b - c

  // CAN'T HAVE v_temp and v_temp2 be the same vector.
  // Results in 0's being stored in v_temp2.
  gsl_blas_dgemv(CblasNoTrans, 1.0, B, v_temp, 0.0, v_temp2);
  //v_temp2 holds A (b-c)

  gsl_blas_ddot(v_temp2, bmc, &d_temp); //d_temp holds (b-c)^T B (b-c)
  result += d_temp;

  result = -0.5 * result;
  result = exp(result);

  result *= sqrt((gr_icov_det * st_icov_det/ApB_det)
                / pow(2*M_PI, MAT_DIM));

  // Freeing memory
  gsl_matrix_free(A);
  gsl_matrix_free(B);
  gsl_matrix_free(ApB);
  gsl_vector_free(a);
  gsl_vector_free(b);
  gsl_vector_free(AapBb);
  gsl_vector_free(c);
  gsl_vector_free(v_temp);
  gsl_vector_free(v_temp2);
  gsl_vector_free(amc);
  gsl_vector_free(bmc);

  gsl_permutation_free(p);

  return result;
}

/* Main function which performs fastest so far:
 * --parameters--
 *  group_icov     (6*6 npyArray) the group's inverse covariance matrix
 *  group_mn       (1*6 npyArray) which is the group's mean kinematic info
 *  group_icov_det (flt)          the determinent of the group_icov
 *  Bs             (nstars*6*6)   an array of each star's icov matrix
 *  bs:            (nstars*6)     an array of each star's mean kinematic info
 *  B_dets:        (nstars)       an array of the determinent of each icov
 *  nstars:        (int)          number of stars, used to determine the size
 *                          of npyArray which will return calculated overlaps)
 *
 * returns: (nstars) array of calculated overlaps of every star with 1 group
 *
 * todo: instead of calling internal function actually use cblas functions
 *          this will save time on the reallocation and deallocation
 */
void get_overlaps(double* gr_icov, int gr_dim1, int gr_dim2,
                  double* gr_mn, int gr_mn_dim,
                  double gr_icov_det,
                  double* st_icovs, int st_dim1, int st_dim2, int st_dim3,
                  double* st_mns, int st_mn_dim1, int st_mn_dim2,
                  double* st_icov_dets, int st_icov_dets_dim,
                  double* rangevec, int n)
{
  // ALLOCATE MEMORY
  int star_count = 0;
  int MAT_DIM = gr_dim1; //Typically set to 6
  int i, j, signum;
  double ApB_det, d_temp, result;
  gsl_permutation *p;

  gsl_matrix *A        = gsl_matrix_alloc(MAT_DIM, MAT_DIM);
  gsl_matrix *B        = gsl_matrix_alloc(MAT_DIM, MAT_DIM);
  gsl_matrix *ApB      = gsl_matrix_alloc(MAT_DIM, MAT_DIM);
  gsl_vector *a        = gsl_vector_alloc(MAT_DIM);
  gsl_vector *b        = gsl_vector_alloc(MAT_DIM);
  gsl_vector *AapBb    = gsl_vector_alloc(MAT_DIM);
  gsl_vector *c        = gsl_vector_alloc(MAT_DIM);
  gsl_vector *v_temp   = gsl_vector_alloc(MAT_DIM);
  gsl_vector *v_temp2  = gsl_vector_alloc(MAT_DIM);
  gsl_vector *amc      = gsl_vector_alloc(MAT_DIM); //will hold a - c
  gsl_vector *bmc      = gsl_vector_alloc(MAT_DIM); //will hold b - c

  p = gsl_permutation_alloc(A->size1);

  // INITIALISE GROUP MATRICES
  for (i=0; i<MAT_DIM; i++)
    for (j=0; j<MAT_DIM; j++)
      gsl_matrix_set (A, i, j, gr_icov[i*MAT_DIM + j]);

  for (i=0; i<MAT_DIM; i++)
    gsl_vector_set (a, i, gr_mn[i]);

  for (star_count=0; star_count<n; star_count++) {
    // INITIALISE STAR MATRICES
    for (i=0; i<MAT_DIM; i++)
      for (j=0; j<MAT_DIM; j++)
        gsl_matrix_set(B,i,j, st_icovs[star_count*MAT_DIM*MAT_DIM+i*MAT_DIM+j]);

    for (i=0; i<MAT_DIM; i++)
      gsl_vector_set (b, i, st_mns[star_count*MAT_DIM + i]);

    // FIND OVERLAP
    // Adding A and B together and storing in ApB
    gsl_matrix_set_zero(ApB);
    gsl_matrix_add(ApB, A);
    gsl_matrix_add(ApB, B);

    // Storing the result A*a + B*b in AapBb
    gsl_vector_set_zero(AapBb);
    gsl_blas_dsymv(CblasUpper, 1.0, A, a, 1.0, AapBb);
    gsl_blas_dsymv(CblasUpper, 1.0, B, b, 1.0, AapBb);

    // Getting determinant of ApB
    gsl_linalg_LU_decomp(ApB, p, &signum);
    //ApB_det = gsl_linalg_LU_det(ApB, signum);
    ApB_det = fabs(gsl_linalg_LU_det(ApB, signum)); //temp doctoring determinant

    // Solve for c
    gsl_linalg_LU_solve(ApB, p, AapBb, c);

    // Compute the overlap formula
    gsl_vector_set_zero(v_temp);
    gsl_blas_dcopy(a, v_temp);       //v_temp holds a
    gsl_blas_daxpy(-1.0, c, v_temp); //v_temp holds a - c
    gsl_blas_dcopy(v_temp, amc);     //amc holds a - c

    // CAN'T HAVE v_temp and v_temp2 be the same vector.
    // Results in 0's being stored in v_temp2.
    gsl_blas_dsymv(CblasUpper, 1.0, A, v_temp, 0.0, v_temp2);
    //v_temp2 holds A (a-c)

    result = 0.0;
    gsl_blas_ddot(v_temp2, amc, &d_temp); //d_temp holds (a-c)^T A (a-c)

    result += d_temp;

    gsl_vector_set_zero(v_temp);
    gsl_blas_dcopy(b, v_temp);       //v_temp holds b
    gsl_blas_daxpy(-1.0, c, v_temp); //v_temp holds b - c
    gsl_blas_dcopy(v_temp, bmc);     //bmc holds b - c

    // CAN'T HAVE v_temp and v_temp2 be the same vector.
    // Results in 0's being stored in v_temp2.
    gsl_blas_dsymv(CblasUpper, 1.0, B, v_temp, 0.0, v_temp2);
    //v_temp2 holds A (b-c)

    gsl_blas_ddot(v_temp2, bmc, &d_temp); //d_temp holds (b-c)^T B (b-c)
    result += d_temp;

    result = -0.5 * result;
    result = exp(result);

    result *= sqrt((gr_icov_det * st_icov_dets[star_count]/ApB_det)
                  / pow(2*M_PI, MAT_DIM));

    // STORE IN 'rangevec'
    rangevec[star_count] = result;
  }

  // DEALLOCATE THE MEMORY
  gsl_matrix_free(A);
  gsl_matrix_free(B);
  gsl_matrix_free(ApB);
  gsl_vector_free(a);
  gsl_vector_free(b);
  gsl_vector_free(AapBb);
  gsl_vector_free(c);
  gsl_vector_free(v_temp);
  gsl_vector_free(v_temp2);
  gsl_vector_free(amc);
  gsl_vector_free(bmc);

  gsl_permutation_free(p);
}

/* New main function, speed not yet tested
 * --parameters--
 *  group_icov     (6*6 npyArray) the group's inverse covariance matrix
 *  group_mn       (1*6 npyArray) which is the group's mean kinematic info
 *  group_icov_det (flt)          the determinent of the group_icov
 *  Bs             (nstars*6*6)   an array of each star's icov matrix
 *  bs:            (nstars*6)     an array of each star's mean kinematic info
 *  B_dets:        (nstars)       an array of the determinent of each icov
 *  nstars:        (int)          number of stars, used to determine the size
 *                          of npyArray which will return calculated overlaps)
 *
 * returns: (nstars) array of calculated overlaps of every star with 1 group
 *
 * todo: instead of calling internal function actually use cblas functions
 *          this will save time on the reallocation and deallocation
 *
 *      look up how to find inverse
 *      look up how to access math.pi
 */
void new_get_lnoverlaps(
  double* gr_cov, int gr_dim1, int gr_dim2,
  double* gr_mn, int gr_mn_dim,
  double* st_covs, int st_dim1, int st_dim2, int st_dim3,
  double* st_mns, int st_mn_dim1, int st_mn_dim2,
  double* rangevec, int n
  )
{
  //printf("Inside new_get_lnoverlaps function\n");
  // ALLOCATE MEMORY
  int star_count = 0;
  int MAT_DIM = gr_dim1; //Typically set to 6
  int i, j, signum;
  double d_temp, result, ln_det_BpA;
  FILE* fout = stdout;
  gsl_permutation *p1;

  gsl_matrix *BpA      = gsl_matrix_alloc(MAT_DIM, MAT_DIM); //(B+A)
  //gsl_matrix *BpAi     = gsl_matrix_alloc(MAT_DIM, MAT_DIM); //(B+A)^-1
  gsl_vector *bma      = gsl_vector_alloc(MAT_DIM); //will hold b - a
  gsl_vector *v_temp   = gsl_vector_alloc(MAT_DIM);

  p1 = gsl_permutation_alloc(BpA->size1);

  //printf("Memory allocated\n");
  for (star_count=0; star_count<n; star_count++) {
    // INITIALISE STAR MATRICES
    for (i=0; i<MAT_DIM; i++)
      for (j=0; j<MAT_DIM; j++)
        //perform B+A as part of the initialisation
        gsl_matrix_set(
          BpA,i,j,
          st_covs[star_count*MAT_DIM*MAT_DIM+i*MAT_DIM+j] +
          gr_cov[i*MAT_DIM+j]
        );
//    printf("Printing BpA\n");
//    print_matrix(fout, BpA);

    for (i=0; i<MAT_DIM; i++) {
      gsl_vector_set(
        bma, i,
        st_mns[star_count*MAT_DIM + i] -
        gr_mn[i]
      );
    }
    //printf("Printing bma\n");
    //print_vec(bma->data, 6);
    //printf("Matrices initialised\n\n");

    result = 6*log(2*M_PI);
    // To Do! put 6ln(2 pi) in here ^^

    // Get inverse of BpA, this line is wrong, fix when have internet
    gsl_linalg_LU_decomp(BpA, p1, &signum);
    ln_det_BpA = log(fabs(gsl_linalg_LU_det(BpA, signum)));
    result += ln_det_BpA;

    //printf("ln(det(BpA)) added\n");
    //printf("%6.2f\n\n",ln_det_BpA);
    //printf("result so far:\n%6.2f\n",result);
    //// Solve for c
    //gsl_linalg_LU_solve(ApB, p, AapBb, c);

    // Don't use invert, use solve like example above ^^
    //gsl_linalg_LU_invert(BpA, p1, BpAi);
    //
    gsl_vector_set_zero(v_temp);
    gsl_linalg_LU_solve(BpA, p1, bma, v_temp); //v_temp holds (B+A)^-1 (b-a)
    gsl_blas_ddot(v_temp, bma, &d_temp); //d_temp holds (b-a)^T (B+A)-1 (b-a)
    //printf("Printing bma_BpAi_bma\n");
    //printf("%6.2f\n\n", d_temp);

    result += d_temp;
    //printf("result after bma_BpAi_bma:\n%6.2f\n",result);

    result *= -0.5;
    //printf("Everything calculated\n");
    //printf("Final result:\n%6.2f\n", result);
    //
    // STORE IN 'rangevec'
    rangevec[star_count] = result;
  }

  // DEALLOCATE THE MEMORY
  gsl_matrix_free(BpA);
  //gsl_matrix_free(BpAi);
  gsl_vector_free(bma);
  gsl_vector_free(v_temp);

  gsl_permutation_free(p1);

  //printf("At end of new_get_lnoverlaps function\n");
}

max_stars_repo_name: tcrundall/chronostar

max_issues_repo_name: tcrundall/chronostar

max_forks_repo_name: tcrundall/chronostar


text:
#ifdef __cplusplus
#include <string>
#include <algorithm>
#include "halley/data_structures/vector.h"
#include <map>
#include <set>
#include <array>
#include <charconv>
#include <gsl/gsl>
#include <cstddef>
#include <limits>
#include <memory>
#include <chrono>
#include <thread>
#include <mutex>
#include <iomanip>
#include <optional>
#include <atomic>
#include <deque>
#include <iostream>
#include <sstream>
#include <istream>
#include <ostream>
#endif

max_stars_repo_name: amrezzd/halley

max_issues_repo_name: amrezzd/halley

max_forks_repo_name: amrezzd/halley


text:
#ifndef SIR_H
#define SIR_H

#include <gsl/gsl_odeiv.h>
#include <gsl/gsl_errno.h>
#include <iostream>
#include "DiffEq_Sim.h"

using namespace std;

class SIR : public DiffEq_Sim {

    private:
        const double beta;
        const double gamma;

    public:
        SIR() : beta(0.0), gamma(0.0) { nbins=3;}
        SIR(double b, double g): beta(b), gamma(g) { nbins=3; }
        ~SIR() {};

        void initialize( double S, double I, double R) {
            y = new double[nbins];
            y[0] = S; y[1] = I; y[2] = R;
        }

        void derivative(double const y[], double dydt[]) {
            dydt[0] = -beta*y[0]*y[1];
            dydt[1] = +beta*y[0]*y[1] - gamma*y[1];
            dydt[2] = +gamma*y[1];
        }

};

#endif


text:
#include <math.h>
#include <galpy_potentials.h>
#include <stdio.h>
#include <gsl/gsl_sf_gegenbauer.h>
#include <gsl/gsl_sf_legendre.h>

#ifndef GSL_MAJOR_VERSION
#define GSL_MAJOR_VERSION 1
#endif
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
//SCF Disk potential
//4 arguments: amp, Acos, Asin, a

const int FORCE =1;
const int DERIV =2;

//Useful Functions

//Converts from cylindrical coordinates to spherical
static inline void cyl_to_spher(double R, double Z,double *r, double *theta)
{
    *r = sqrt(R*R + Z*Z);
    *theta = atan2(R, Z);
}

//Integer power
double power(double x, int i)
{
    if (i==0)
        return 1;
    return x*power(x,i - 1);
}

//Calculates xi
static inline void calculateXi(double r, double a, double *xi)
{
    *xi = (r - a)/(r + a);
}

//Potentials, forces, and derivative functions
//LCOV_EXCL_START
double computePhi(double Acos_val, double Asin_val, double mCos, double mSin, double P, double phiTilde, int m)
{
    return (Acos_val*mCos + Asin_val*mSin)*P*phiTilde;
}
//LCOV_EXCL_STOP
double computeAxiPhi(double Acos_val, double P, double phiTilde)
{
    return Acos_val*P*phiTilde;
}

double computeF_r(double Acos_val, double Asin_val, double mCos, double mSin, double P, double dphiTilde, int m)
{
    return -(Acos_val*mCos + Asin_val*mSin)*P*dphiTilde;
}
double computeAxiF_r(double Acos_val, double P, double dphiTilde)
{
    return -Acos_val*P*dphiTilde;
}

double computeF_theta(double Acos_val, double Asin_val, double mCos, double mSin, double dP, double phiTilde, int m)
{
    return -(Acos_val*mCos + Asin_val*mSin)*dP*phiTilde;
}
double computeAxiF_theta(double Acos_val, double dP, double phiTilde)
{
    return -Acos_val*dP*phiTilde;
}

double computeF_phi(double Acos_val, double Asin_val, double mCos, double mSin, double P, double phiTilde, int m)
{
    return m*(Acos_val*mSin - Asin_val*mCos)*P*phiTilde;
}
double computeAxiF_phi(double Acos_val, double P, double phiTilde)
{
    return 0.;
}

double computeF_rr(double Acos_val, double Asin_val, double mCos, double mSin, double P, double d2phiTilde, int m)
{
    return -(Acos_val*mCos + Asin_val*mSin)*P*d2phiTilde;
}
double computeAxiF_rr(double Acos_val, double P, double d2phiTilde)
{
    return -Acos_val*P*d2phiTilde;
}

double computeF_rphi(double Acos_val, double Asin_val, double mCos, double mSin, double P, double dphiTilde, int m)
{
    return m*(Acos_val*mSin - Asin_val*mCos)*P*dphiTilde;
}
double computeAxiF_rphi(double Acos_val, double P, double d2phiTilde)
{
    return 0.;
}

double computeF_phiphi(double Acos_val, double Asin_val, double mCos, double mSin, double P, double phiTilde, int m)
{
    return m*m*(Acos_val*mCos + Asin_val*mSin)*P*phiTilde;
}
double computeAxiF_phiphi(double Acos_val, double P, double d2phiTilde)
{
    return 0.;
}

//Calculates the Gegenbauer polynomials
void compute_C(double xi, int N, int L, double * C_array)
{
    int l;
    for (l = 0; l < L; l++)
    {
        gsl_sf_gegenpoly_array(N - 1, 3./2 + 2*l, xi, C_array + l*N);
    }
}

//Calculates the derivative of the Gegenbauer polynomials
void compute_dC(double xi, int N, int L, double * dC_array)
{
    int n,l;
    for (l = 0; l < L; l++)
    {
        *(dC_array +l*N) = 0;
        if (N != 1)
            gsl_sf_gegenpoly_array(N - 2, 5./2 + 2*l, xi, dC_array + l*N + 1);
        for (n = 0; n<N; n++)
        {
            *(dC_array +l*N + n) *= 2*(2*l + 3./2);
        }
    }

}

//Calculates the second derivative of the Gegenbauer polynomials
void compute_d2C(double xi, int N, int L, double * d2C_array)
{
    int n,l;
    for (l = 0; l < L; l++)
    {
        *(d2C_array +l*N) = 0;
        if (N >1)
            *(d2C_array +l*N + 1) = 0;
        if (N > 2)
            gsl_sf_gegenpoly_array(N - 3, 7./2 + 2*l, xi, d2C_array + l*N + 2);
        for (n = 0; n<N; n++)
        {
            *(d2C_array +l*N + n) *= 4*(2*l + 3./2)*(2*l + 5./2);
        }
    }

}

//Compute phi_Tilde
void compute_phiTilde(double r, double a, int N, int L, double* C, double * phiTilde)
{
    double xi;
    calculateXi(r, a, &xi);
    double rterms = -1./(r + a);
    int n,l;
    for (l = 0; l < L; l++)
    {

        if (l != 0)
            rterms *= r*a/((a + r)*(a + r));

        for (n = 0; n < N; n++)
        {
            *(phiTilde + l*N + n)  = rterms*(*(C + n + l*N));
        }
    }

}

//Computes the derivative of phiTilde with respect to r
void compute_dphiTilde(double r, double a, int N, int L, double * C, double * dC, double * dphiTilde)
{
    double xi;
    calculateXi(r, a, &xi);
    double rterm = 1./(r*power(a + r, 3));
    int n,l;
    for (l = 0; l < L; l++)
    {
        if (l != 0)
        {
            rterm *= (a*r)/power(a + r, 2) ;
        }
        for (n = 0; n < N; n++)
        {
            *( dphiTilde + l*N + n) = rterm *(((2*l + 1)*r*(a + r) - l*power(a + r,2))*(*(C + l*N + n)) -
                                              2*a*r*(*(dC + l*N + n)));

        }
    }
}

//Computes the second derivative of phiTilde with respect to r
void compute_d2phiTilde(double r, double a, int N, int L, double * C, double * dC,double * d2C, double * d2phiTilde)
{
    double xi;
    calculateXi(r, a, &xi);
    double rterm = 1./(r*r) / power(a + r,5);
    int n,l;
    for (l = 0; l < L; l++)
    {

        if (l != 0)
        {
            rterm *= (a*r)/power(a + r, 2);
        }
        for (n = 0; n < N; n++)
        {

            double C_val = *(C + l*N + n);
            double dC_val = *(dC + l*N + n);
            double d2C_val = *(d2C + l*N + n);
            *( d2phiTilde + l*N + n) = rterm*(C_val*(l*(1 - l)*power(a + r, 4) - (4*l*l + 6*l + 2.)*r*r*power(a + r,2) +
                                              l*(4*l + 2)*r*power(a + r,3))
                                              + a*r*((4*r*r + 4*a*r + (8*l + 4)*r*(a + r)-
                                                      4*l*power(a + r,2))*dC_val - 4*a*r*d2C_val));

        }
    }
}

//Computes the associated Legendre polynomials
void compute_P(double x, int L, int M, double * P_array)
{
    if (M == 1){
        gsl_sf_legendre_Pl_array (L - 1, x, P_array);
    } else {
        #if GSL_MAJOR_VERSION == 2
            gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,L - 1, x, -1, P_array);
        #else
            int m;
            for (m = 0; m < M; m++)
            {
                gsl_sf_legendre_Plm_array(L - 1, m, x, P_array);
                P_array += L - m;
            }
        #endif
    }

}

//Computes the associated Legendre polynomials and its derivative
void compute_P_dP(double x, int L, int M, double * P_array, double *dP_array)
{
    if (M == 1){
        gsl_sf_legendre_Pl_deriv_array (L - 1, x, P_array, dP_array);
    } else {
        #if GSL_MAJOR_VERSION == 2
            gsl_sf_legendre_deriv_array_e(GSL_SF_LEGENDRE_NONE, L - 1, x, -1,P_array, dP_array);

        #else
            int m;
            for (m = 0; m < M; m++)
            {
                gsl_sf_legendre_Plm_deriv_array(L - 1, m, x, P_array, dP_array);
                P_array += L - m;
                dP_array += L - m;
            }
        #endif
    }
}

typedef struct equations equations;
struct equations
{
    double ((**Eq)(double, double, double, double, double, double, int));
    double *(*phiTilde);
    double *(*P);
    double *Constant;
};

typedef struct axi_equations axi_equations;
struct axi_equations
{
    double ((**Eq)(double, double, double));
    double *(*phiTilde);
    double *(*P);
    double *Constant;
};

//Compute axi symmetric potentials.
void compute(double a, int N, int L, int M,
	     double r, double theta, double phi,
	     double *Acos, int eq_size,
	     axi_equations e,
	     double *F)
{
    int i,n,l;
    for (i = 0; i < eq_size; i++)
    {
        *(F + i) =0; //Initialize each F
    }

    for (l = 0; l < L; l++)
    {

        for (n = 0; n < N; n++)
        {

            double Acos_val = *(Acos + M*l + M*L*n);
            for (i = 0; i < eq_size; i++)
            {
                double (*Eq)(double, double, double) = *(e.Eq + i);
                double *P = *(e.P + i);
                double *phiTilde = *(e.phiTilde + i);
                *(F + i) += (*Eq)(Acos_val, P[l], phiTilde[l*N + n]);
            }

        }
    }

    //Multiply F by constants
    for (i = 0; i < eq_size; i++)
    {
        double constant = *(e.Constant + i);
        *(F + i) *= constant*sqrt(4*M_PI);
    }

}

//Compute Non Axi symmetric Potentials
void computeNonAxi(double a, int N, int L, int M,
		   double r, double theta, double phi,
		   double *Acos, double *Asin, int eq_size,
		   equations e,
		   double *F)
{
    int i,n,l,m;
    for (i = 0; i < eq_size; i++)
    {
        *(F + i) =0; //Initialize each F
    }

    int v1counter=0;
    int v2counter = 0;
    int Pindex = 0;
    for (l = 0; l < L; l++)
    {
        v1counter = 0;
        for (m = 0; m<=l; m++)
        {
            #if GSL_MAJOR_VERSION == 2
                Pindex = v2counter;
            #else
                Pindex = v1counter + l;
            #endif

            double mCos = cos(m*phi);
            double mSin = sin(m*phi);
            for (n = 0; n < N; n++)
            {
                double Acos_val = *(Acos +m + M*l + M*L*n);
                double Asin_val = *(Asin +m + M*l + M*L*n);
                for (i = 0; i < eq_size; i++)
                {
                    double (*Eq)(double, double, double, double, double, double, int) = *(e.Eq + i);
                    double *P = *(e.P + i);
                    double *phiTilde = *(e.phiTilde + i);
                    *(F + i) += (*Eq)(Acos_val, Asin_val, mCos, mSin, P[Pindex], phiTilde[l*N + n], m);
                }

            }

                v2counter++;
                v1counter+=M-m - 1;

        }

    }
    //Multiply F by constants
    for (i = 0; i < eq_size; i++)
    {
        double constant = *(e.Constant + i);
        *(F + i) *= constant*sqrt(4*M_PI);
    }

}

//Compute the Forces
void computeForce(double R,double Z, double phi,
		  double t,
		  struct potentialArg * potentialArgs, double * F)
{
    double * args= potentialArgs->args;
    //Get args
    double a = *args++;
    int isNonAxi = (int)*args++;
    int N = (int)*args++;
    int L = (int)*args++;
    int M = (int)*args++;

    double* Acos = args;

    double* caching_i = (args + (isNonAxi + 1)*N*L*M);
    double *Asin;
    if (isNonAxi == 1)
    {
        Asin = args + N*L*M;
    }
    double *cached_type = caching_i;
    double * cached_coords = (caching_i+ 1);
    double * cached_values = (caching_i + 4);
    if ((int)*cached_type==FORCE)
    {
        if (*cached_coords == R && *(cached_coords + 1) == Z && *(cached_coords + 2) == phi)
        {
            *F = *cached_values;
            *(F + 1) = *(cached_values + 1);
            *(F + 2) = *(cached_values + 2);
            return;
        }
    }
    double r;
    double theta;
    cyl_to_spher(R, Z, &r, &theta);

    double xi;
    calculateXi(r, a, &xi);

//Compute the gegenbauer polynomials and its derivative.
    double *C= (double *) malloc ( N*L * sizeof(double) );
    double *dC= (double *) malloc ( N*L * sizeof(double) );
    double *phiTilde= (double *) malloc ( N*L * sizeof(double) );
    double *dphiTilde= (double *) malloc ( N*L * sizeof(double) );

    compute_C(xi, N, L, C);
    compute_dC(xi, N, L, dC);

//Compute phiTilde and its derivative
    compute_phiTilde(r, a, N, L, C, phiTilde);

    compute_dphiTilde(r, a, N, L, C, dC, dphiTilde);

//Compute Associated Legendre Polynomials
    int M_eff = M;
    int size = 0;

    if (isNonAxi==0)
    {
    M_eff = 1;
    size = L;
    } else{
    size = L*L - L*(L-1)/2;
    }

    double *P= (double *) malloc ( size * sizeof(double) );
    double *dP= (double *) malloc ( size * sizeof(double) );
    compute_P_dP(cos(theta), L, M_eff, P, dP);

    double (*PhiTilde_Pointer[3]) = {dphiTilde,phiTilde,phiTilde};
    double (*P_Pointer[3]) = {P, dP, P};

    double Constant[3] = {1., -sin(theta), 1.};

    if (isNonAxi == 1)
    {
        double (*Eq[3])(double, double, double, double, double, double, int) = {&computeF_r, &computeF_theta, &computeF_phi};
        equations e = {Eq,&PhiTilde_Pointer[0], &P_Pointer[0], &Constant[0]};
        computeNonAxi(a, N, L, M,r, theta, phi, Acos, Asin, 3, e, F);
    }
    else
    {
        double (*Eq[3])(double, double, double) = {&computeAxiF_r, &computeAxiF_theta, &computeAxiF_phi};
        axi_equations e = {Eq,&PhiTilde_Pointer[0], &P_Pointer[0], &Constant[0]};
        compute(a, N, L, M,r, theta, phi, Acos, 3, e, F);
    }

    //Caching

    *cached_type = (double)FORCE;

    * cached_coords = R;
    * (cached_coords + 1) = Z;
    * (cached_coords + 2) = phi;
    * (cached_values) = *F;
    * (cached_values + 1) = *(F + 1);
    * (cached_values + 2) = *(F + 2);

    // Free memory
    free(C);
    free(dC);
    free(phiTilde);
    free(dphiTilde);
    free(P);
    free(dP);

}

//Compute the Derivatives
void computeDeriv(double R,double Z, double phi,
		  double t,
		  struct potentialArg * potentialArgs, double * F)
{
    double * args= potentialArgs->args;
    //Get args
    double a = *args++;
    int isNonAxi = (int)*args++;
    int N = (int) *args++;
    int L = (int) *args++;
    int M = (int) *args++;
    double* Acos = args;

    double * caching_i = (args + (isNonAxi + 1)*N*L*M);
    double *Asin;
    if (isNonAxi == 1)
    {
        Asin = args + N*L*M;
    }

    double *cached_type = caching_i;
    double * cached_coords = (caching_i+ 1);
    double * cached_values = (caching_i + 4);
    if ((int)*cached_type==DERIV)
    {
        if (*cached_coords == R && *(cached_coords + 1) == Z && *(cached_coords + 2) == phi)
        {
            *F = *cached_values;
            *(F + 1) = *(cached_values + 1);
            *(F + 2) = *(cached_values + 2);
            return;
        }
    }

    double r;
    double theta;
    cyl_to_spher(R, Z, &r, &theta);

    double xi;
    calculateXi(r, a, &xi);

//Compute the gegenbauer polynomials and its derivative.
    double *C= (double *) malloc ( N*L * sizeof(double) );
    double *dC= (double *) malloc ( N*L * sizeof(double) );
    double *d2C= (double *) malloc ( N*L * sizeof(double) );
    double *phiTilde= (double *) malloc ( N*L * sizeof(double) );
    double *dphiTilde= (double *) malloc ( N*L * sizeof(double) );
    double *d2phiTilde= (double *) malloc ( N*L * sizeof(double) );

    compute_C(xi, N, L, C);
    compute_dC(xi, N, L, dC);
    compute_d2C(xi, N, L, d2C);

//Compute phiTilde and its derivative
    compute_phiTilde(r, a, N, L, C, phiTilde);
    compute_dphiTilde(r, a, N, L, C, dC, dphiTilde);
    compute_d2phiTilde(r, a, N, L, C, dC, d2C, d2phiTilde);

//Compute Associated Legendre Polynomials
    int M_eff = M;
    int size = 0;

    if (isNonAxi==0)
    {
    M_eff = 1;
    size = L;
    } else{
    size = L*L - L*(L-1)/2;
    }
    double *P= (double *) malloc ( size * sizeof(double) );

    compute_P(cos(theta), L,M_eff, P);

    double (*PhiTilde_Pointer[3])= {d2phiTilde,phiTilde,dphiTilde};
    double (*P_Pointer[3]) = {P, P, P};

    double Constant[3] = {1., 1., 1.};

    if (isNonAxi==1)
    {
        double (*Eq[3])(double, double, double, double, double, double, int) = {&computeF_rr, &computeF_phiphi, &computeF_rphi};
        equations e = {Eq,&PhiTilde_Pointer[0],&P_Pointer[0],&Constant[0]};
        computeNonAxi(a, N, L, M,r, theta, phi, Acos, Asin, 3, e, F);
    }
    else
    {
        double (*Eq[3])(double, double, double) = {&computeAxiF_rr, &computeAxiF_phiphi, &computeAxiF_rphi};
        axi_equations e = {Eq,&PhiTilde_Pointer[0],&P_Pointer[0],&Constant[0]};
        compute(a, N, L, M,r, theta, phi, Acos, 3, e, F);
    }

    //Caching

    *cached_type = (double)DERIV;

    * cached_coords = R;
    * (cached_coords + 1) = Z;
    * (cached_coords + 2) = phi;
    * (cached_values) = *F;
    * (cached_values + 1) = *(F + 1);
    * (cached_values + 2) = *(F + 2);

    //Free memory
    free(C);
    free(dC);
    free(d2C);
    free(phiTilde);
    free(dphiTilde);
    free(d2phiTilde);
    free(P);
}

//Compute the Potential
double SCFPotentialEval(double R,double Z, double phi,
                        double t,
                        struct potentialArg * potentialArgs)
{
    double * args= potentialArgs->args;
    //Get args
    double a = *args++;
    int isNonAxi = (int)*args++;
    int N = (int) *args++;
    int L = (int) *args++;
    int M = (int) *args++;
    double* Acos = args;
    double* Asin;
    if (isNonAxi==1) //LCOV_EXCL_START
    {
        Asin = args + N*L*M;
    } //LCOV_EXCL_STOP
    //convert R,Z to r, theta
    double r;
    double theta;
    cyl_to_spher(R, Z,&r, &theta);
    double xi;
    calculateXi(r, a, &xi);

    //Compute the gegenbauer polynomials and its derivative.
    double *C= (double *) malloc ( N*L * sizeof(double) );
    double *phiTilde= (double *) malloc ( N*L * sizeof(double) );

    compute_C(xi, N, L, C);

    //Compute phiTilde and its derivative
    compute_phiTilde(r, a, N, L, C, phiTilde);
    //Compute Associated Legendre Polynomials

    int M_eff = M;
    int size = 0;

    if (isNonAxi==0)
    {
    M_eff = 1;
    size = L;
    } else{ //LCOV_EXCL_START
    size = L*L - L*(L-1)/2;
    } //LCOV_EXCL_STOP

    double *P= (double *) malloc ( size * sizeof(double) );

    compute_P(cos(theta), L,M_eff, P);

    double potential;

    double (*PhiTilde_Pointer[1]) = {phiTilde};
    double (*P_Pointer[1]) = {P};

    double Constant[1] = {1.};

    if (isNonAxi==1) //LCOV_EXCL_START
    {
        double (*Eq[1])(double, double, double, double, double, double, int) = {&computePhi};
        equations e = {Eq,&PhiTilde_Pointer[0],&P_Pointer[0],&Constant[0]};
        computeNonAxi(a, N, L, M,r, theta, phi, Acos, Asin, 1, e, &potential);
    } //LCOV_EXCL_STOP
    else
    {
        double (*Eq[1])(double, double, double) = {&computeAxiPhi};
        axi_equations e = {Eq,&PhiTilde_Pointer[0],&P_Pointer[0],&Constant[0]};
        compute(a, N, L, M,r, theta, phi, Acos, 1, e, &potential);
    }

    //Free memory
    free(C);
    free(phiTilde);
    free(P);

    return potential;

}

//Compute the force in the R direction
double SCFPotentialRforce(double R,double Z, double phi,
                          double t,
                          struct potentialArg * potentialArgs)
{
    double r;
    double theta;
    cyl_to_spher(R, Z, &r, &theta);
    // The derivatives
    double dr_dR = R/r;
    double dtheta_dR = Z/(r*r);
    double dphi_dR = 0;

    double F[3];
    computeForce(R, Z, phi, t,potentialArgs, &F[0]) ;

    return *(F + 0)*dr_dR + *(F + 1)*dtheta_dR + *(F + 2)*dphi_dR;

}

//Compute the force in the z direction
double SCFPotentialzforce(double R,double Z, double phi,
                          double t,
                          struct potentialArg * potentialArgs)
{
    double r;
    double theta;
    cyl_to_spher(R, Z,&r, &theta);

    double dr_dz = Z/r;
    double dtheta_dz = -R/(r*r);
    double dphi_dz = 0;

    double F[3];
    computeForce(R, Z, phi, t,potentialArgs, &F[0]) ;
    return *(F + 0)*dr_dz + *(F + 1)*dtheta_dz + *(F + 2)*dphi_dz;
}

//Compute the force in the phi direction
double SCFPotentialphiforce(double R,double Z, double phi,
                            double t,
                            struct potentialArg * potentialArgs)
{

    double r;
    double theta;
    cyl_to_spher(R, Z, &r, &theta);

    double dr_dphi = 0;
    double dtheta_dphi = 0;
    double dphi_dphi = 1;

    double F[3];
    computeForce(R, Z, phi, t,potentialArgs, &F[0]) ;

    return *(F + 0)*dr_dphi + *(F + 1)*dtheta_dphi + *(F + 2)*dphi_dphi;
}

//Compute the planar force in the R direction
double SCFPotentialPlanarRforce(double R,double phi,
                                double t,
                                struct potentialArg * potentialArgs)
{
    return SCFPotentialRforce(R,0., phi,t,potentialArgs);

}

//Compute the planar force in the phi direction
double SCFPotentialPlanarphiforce(double R,double phi,
                                  double t,
                                  struct potentialArg * potentialArgs)
{
    return SCFPotentialphiforce(R,0., phi,t,potentialArgs);
}

//Compute the planar double derivative of the potential with respect to R
double SCFPotentialPlanarR2deriv(double R, double phi,
                                 double t,
                                 struct potentialArg * potentialArgs)
{
    double Farray[3];
    computeDeriv(R, 0, phi, t,potentialArgs, &Farray[0]) ;
    return *Farray;
}

//Compute the planar double derivative of the potential with respect to phi
double SCFPotentialPlanarphi2deriv(double R, double phi,
                                   double t,
                                   struct potentialArg * potentialArgs)
{
    double Farray[3];
    computeDeriv(R, 0, phi, t,potentialArgs, &Farray[0]) ;
    return *(Farray + 1);
}

//Compute the planar double derivative of the potential with respect to R, Phi
double SCFPotentialPlanarRphideriv(double R, double phi,
                                   double t,
                                   struct potentialArg * potentialArgs)
{
    double Farray[3];
    computeDeriv(R, 0, phi, t,potentialArgs, &Farray[0]) ;
    return *(Farray + 2);
}

max_stars_repo_name: turnergarrow/galpy

max_issues_repo_name: BurcuAkbulut/galpy

max_forks_repo_name: BurcuAkbulut/galpy