text:
/* ode-initval/test_odeiv.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * 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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Author:  G. Jungman
 */
#include <config.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_ieee_utils.h>
#include "gsl_odeiv.h"

int rhs_linear (double t, const double y[], double f[], void *params);
int jac_linear (double t, const double y[], double *dfdy, double dfdt[],
                void *params);
int rhs_sin (double t, const double y[], double f[], void *params);
int jac_sin (double t, const double y[], double *dfdy, double dfdt[],
             void *params);
int rhs_exp (double t, const double y[], double f[], void *params);
int jac_exp (double t, const double y[], double *dfdy, double dfdt[],
             void *params);
int rhs_stiff (double t, const double y[], double f[], void *params);
int jac_stiff (double t, const double y[], double *dfdy, double dfdt[],
               void *params);
void test_stepper_linear (const gsl_odeiv_step_type * T, double h,
                         double base_prec);
void test_stepper_sin (const gsl_odeiv_step_type * T, double h, double base_prec);
void test_stepper_exp (const gsl_odeiv_step_type * T, double h, double base_prec);
void test_stepper_stiff (const gsl_odeiv_step_type * T, double h, double base_prec);
void test_evolve_system_flat (gsl_odeiv_step * step,
                              const gsl_odeiv_system * sys,
                              double t0, double t1, double hstart,
                              double y[], double yfin[],
                              double err_target, const char *desc);

void test_evolve_system (const gsl_odeiv_step_type * T,
                         const gsl_odeiv_system * sys,
                         double t0, double t1, double hstart,
                         double y[], double yfin[],
                         double err_target, const char *desc);

void test_evolve_exp (const gsl_odeiv_step_type * T, double h, double err);
void test_evolve_sin (const gsl_odeiv_step_type * T, double h, double err);
void test_evolve_stiff1 (const gsl_odeiv_step_type * T, double h, double err);
void test_evolve_stiff5 (const gsl_odeiv_step_type * T, double h, double err);

/* RHS for a + b t */

int
rhs_linear (double t, const double y[], double f[], void *params)
{
  f[0] = 0.0;
  f[1] = y[0];
  return GSL_SUCCESS;
}

int
jac_linear (double t, const double y[], double *dfdy, double dfdt[],
	    void *params)
{
  gsl_matrix dfdy_mat;
  dfdy_mat.size1 = 2;
  dfdy_mat.size2 = 2;
  dfdy_mat.tda = 2;
  dfdy_mat.data = dfdy;
  dfdy_mat.block = 0;
  gsl_matrix_set (&dfdy_mat, 0, 0, 0.0);
  gsl_matrix_set (&dfdy_mat, 0, 1, 0.0);
  gsl_matrix_set (&dfdy_mat, 1, 0, 1.0);
  gsl_matrix_set (&dfdy_mat, 1, 1, 0.0);
  dfdt[0] = 0.0;
  dfdt[1] = 0.0;
  return GSL_SUCCESS;
}

gsl_odeiv_system rhs_func_lin = {
  rhs_linear,
  jac_linear,
  2,
  0
};

/* RHS for sin(t),cos(t) */

int
rhs_sin (double t, const double y[], double f[], void *params)
{
  f[0] = -y[1];
  f[1] = y[0];
  return GSL_SUCCESS;
}

int
jac_sin (double t, const double y[], double *dfdy, double dfdt[],
	 void *params)
{
  gsl_matrix dfdy_mat;
  dfdy_mat.data = dfdy;
  dfdy_mat.size1 = 2;
  dfdy_mat.size2 = 2;
  dfdy_mat.tda = 2;
  dfdy_mat.block = 0;
  gsl_matrix_set (&dfdy_mat, 0, 0, 0.0);
  gsl_matrix_set (&dfdy_mat, 0, 1, -1.0);
  gsl_matrix_set (&dfdy_mat, 1, 0, 1.0);
  gsl_matrix_set (&dfdy_mat, 1, 1, 0.0);
  dfdt[0] = 0.0;
  dfdt[1] = 0.0;
  return GSL_SUCCESS;
}

gsl_odeiv_system rhs_func_sin = {
  rhs_sin,
  jac_sin,
  2,
  0
};

/* RHS for a exp(t)+ b exp(-t) */

int
rhs_exp (double t, const double y[], double f[], void *params)
{
  f[0] = y[1];
  f[1] = y[0];
  return GSL_SUCCESS;
}

int
jac_exp (double t, const double y[], double *dfdy, double dfdt[],
	 void *params)
{
  gsl_matrix dfdy_mat;
  dfdy_mat.data = dfdy;
  dfdy_mat.size1 = 2;
  dfdy_mat.size2 = 2;
  dfdy_mat.tda = 2;
  dfdy_mat.block = 0;
  gsl_matrix_set (&dfdy_mat, 0, 0, 0.0);
  gsl_matrix_set (&dfdy_mat, 0, 1, 1.0);
  gsl_matrix_set (&dfdy_mat, 1, 0, 1.0);
  gsl_matrix_set (&dfdy_mat, 1, 1, 0.0);
  dfdt[0] = 0.0;
  dfdt[1] = 0.0;
  return GSL_SUCCESS;
}

gsl_odeiv_system rhs_func_exp = {
  rhs_exp,
  jac_exp,
  2,
  0
};

/* RHS for stiff example */

int
rhs_stiff (double t, const double y[], double f[], void *params)
{
  f[0] = 998.0 * y[0] + 1998.0 * y[1];
  f[1] = -999.0 * y[0] - 1999.0 * y[1];
  return GSL_SUCCESS;
}

int
jac_stiff (double t, const double y[], double *dfdy, double dfdt[],
	   void *params)
{
  gsl_matrix dfdy_mat;
  dfdy_mat.data = dfdy;
  dfdy_mat.size1 = 2;
  dfdy_mat.size2 = 2;
  dfdy_mat.tda = 2;
  dfdy_mat.block = 0;
  gsl_matrix_set (&dfdy_mat, 0, 0, 998.0);
  gsl_matrix_set (&dfdy_mat, 0, 1, 1998.0);
  gsl_matrix_set (&dfdy_mat, 1, 0, -999.0);
  gsl_matrix_set (&dfdy_mat, 1, 1, -1999.0);
  dfdt[0] = 0.0;
  dfdt[1] = 0.0;
  return GSL_SUCCESS;
}

gsl_odeiv_system rhs_func_stiff = {
  rhs_stiff,
  jac_stiff,
  2,
  0
};

void
test_stepper_linear (const gsl_odeiv_step_type * T, double h,
		     double base_prec)
{
  int s = 0;
  double y[2];
  double yerr[2];
  double t;
  double del;
  double delmax = 0.0;
  int count = 0;

  gsl_odeiv_step *stepper = gsl_odeiv_step_alloc (T, 2);

  y[0] = 1.0;
  y[1] = 0.0;

  for (t = 0.0; t < 4.0; t += h)
    {
      gsl_odeiv_step_apply (stepper, t, h, y, yerr, 0, 0, &rhs_func_lin);
      del = fabs ((y[1] - (t + h)) / y[1]);
      delmax = GSL_MAX_DBL (del, delmax);
      if (del > (count + 1.0) * base_prec)
	{
	  printf ("  LINEAR(%20.17g)  %20.17g  %20.17g  %8.4g\n", t + h, y[1],
		  t + h, del);
	  s++;
	}
      count++;
    }

  gsl_test (s, "%s, linear [0,4], max relative error = %g",
	    gsl_odeiv_step_name (stepper), delmax);

  gsl_odeiv_step_free (stepper);
}

void
test_stepper_sin (const gsl_odeiv_step_type * T, double h, double base_prec)
{
  int s = 0;
  double y[2];
  double yerr[2];
  double t;
  double del;
  double delmax = 0.0;
  int count = 0;

  gsl_odeiv_step *stepper = gsl_odeiv_step_alloc (T, 2);

  y[0] = 1.0;
  y[1] = 0.0;

  for (t = 0.0; t < M_PI; t += h)
    {
      int stat;
      double sin_th = sin (t + h);
      gsl_odeiv_step_apply (stepper, t, h, y, yerr, 0, 0, &rhs_func_sin);
      del = fabs ((y[1] - sin_th) / sin_th);
      delmax = GSL_MAX_DBL (del, delmax);
      {
	if (t < 0.5 * M_PI)
	  {
	    stat = (del > (count + 1.0) * base_prec);
	  }
	else if (t < 0.7 * M_PI)
	  {
	    stat = (del > 1.0e+04 * base_prec);
	  }
	else if (t < 0.9 * M_PI)
	  {
	    stat = (del > 1.0e+06 * base_prec);
	  }
	else
	  {
	    stat = (del > 1.0e+09 * base_prec);
	  }
	if (stat != 0)
	  {
	    printf ("  SIN(%22.18g)  %22.18g  %22.18g  %10.6g\n", t + h, y[1],
		    sin_th, del);
	  }
	s += stat;
      }
      count++;
    }
  if (delmax > 1.0e+09 * base_prec)
    {
      s++;
      printf ("  SIN(0 .. M_PI)  delmax = %g\n", delmax);
    }

  gsl_test (s, "%s, sine [0,pi], max relative error = %g",
	    gsl_odeiv_step_name (stepper), delmax);

  delmax = 0.0;
  for (; t < 100.5 * M_PI; t += h)
    {
      gsl_odeiv_step_apply (stepper, t, h, y, yerr, 0, 0, &rhs_func_sin);
      del = fabs (y[1] - sin (t));
      delmax = GSL_MAX_DBL (del, delmax);
      count++;
    }

  if (del > count * 2.0 * base_prec)
    {
      s++;
      printf ("  SIN(%22.18g)  %22.18g  %22.18g  %10.6g\n", t + h, y[1],
	      sin (t), del);
    }

  gsl_test (s, "%s, sine [pi,100.5*pi], max absolute error = %g",
	    gsl_odeiv_step_name (stepper), delmax);

  gsl_odeiv_step_free (stepper);
}

void
test_stepper_exp (const gsl_odeiv_step_type * T, double h, double base_prec)
{
  int s = 0;
  double y[2];
  double yerr[2];
  double t;
  double del, delmax = 0.0;
  int count = 0;

  gsl_odeiv_step *stepper = gsl_odeiv_step_alloc (T, 2);

  y[0] = 1.0;
  y[1] = 1.0;

  for (t = 0.0; t < 20.0; t += h)
    {
      double ex = exp (t + h);
      gsl_odeiv_step_apply (stepper, t, h, y, yerr, 0, 0, &rhs_func_exp);
      del = fabs ((y[1] - ex) / y[1]);
      delmax = GSL_MAX_DBL (del, delmax);
      if (del > (count + 1.0) * 2.0 * base_prec)
	{
	  printf ("  EXP(%20.17g)  %20.17g  %20.17g  %8.4g\n", t + h, y[1],
		  ex, del);
	  s++;
	}
      count++;
    }

  gsl_test (s, "%s, exponential [0,20], max relative error = %g",
	    gsl_odeiv_step_name (stepper), delmax);

  gsl_odeiv_step_free (stepper);
}

void
test_stepper_stiff (const gsl_odeiv_step_type * T, double h, double base_prec)
{
  int s = 0;
  double y[2];
  double yerr[2];
  double t;
  double del, delmax = 0.0;
  int count = 0;

  gsl_odeiv_step *stepper = gsl_odeiv_step_alloc (T, 2);

  y[0] = 1.0;
  y[1] = 0.0;

  for (t = 0.0; t < 20.0; t += h)
    {
      gsl_odeiv_step_apply (stepper, t, h, y, yerr, NULL, NULL,
			    &rhs_func_stiff);

      if (t > 0.04)
	{
	  double arg = t + h;
	  double e1 = exp (-arg);
	  double e2 = exp (-1000.0 * arg);
	  double u = 2.0 * e1 - e2;
	  /* double v = -e1 + e2; */
	  del = fabs ((y[0] - u) / y[0]);
	  delmax = GSL_MAX_DBL (del, delmax);

	  if (del > (count + 1.0) * 100.0 * base_prec)
	    {
	      printf ("  STIFF(%20.17g)  %20.17g  %20.17g  %8.4g\n", arg,
		      y[0], u, del);
	      s++;
	    }
	}
      count++;
    }

  gsl_test (s, "%s, stiff [0,20], max relative error = %g",
	    gsl_odeiv_step_name (stepper), delmax);

  gsl_odeiv_step_free (stepper);
}

void
test_evolve_system_flat (gsl_odeiv_step * step,
			 const gsl_odeiv_system * sys,
			 double t0, double t1, double hstart,
			 double y[], double yfin[],
			 double err_target, const char *desc)
{
  int s = 0;
  double frac;

  double t = t0;
  double h = hstart;

  gsl_odeiv_evolve *e = gsl_odeiv_evolve_alloc (sys->dimension);

  while (t < t1)
    {
      gsl_odeiv_evolve_apply (e, NULL, step, sys, &t, t1, &h, y);
    }

  frac = fabs ((y[1] - yfin[1]) / yfin[1]) + fabs ((y[0] - yfin[0]) / yfin[0]);

  if (frac > 2.0 * e->count * err_target)
    {
      printf ("FLAT t = %.5e  y0 = %g y1= %g\n", t, y[0], y[1]);
      s++;
    }

  gsl_test (s, "%s, %s, evolve, no control, max relative error = %g",
	    gsl_odeiv_step_name (step), desc, frac);

  gsl_odeiv_evolve_free (e);
}

void
test_evolve_system (const gsl_odeiv_step_type * T,
		    const gsl_odeiv_system * sys,
		    double t0, double t1, double hstart,
		    double y[], double yfin[],
		    double err_target, const char *desc)
{
  int s = 0;
  double frac;

  double t = t0;
  double h = hstart;

  gsl_odeiv_step * step = gsl_odeiv_step_alloc (T, sys->dimension);

  gsl_odeiv_control *c =
    gsl_odeiv_control_standard_new (0.0, err_target, 1.0, 1.0);
  gsl_odeiv_evolve *e = gsl_odeiv_evolve_alloc (sys->dimension);

  while (t < t1)
    {
      gsl_odeiv_evolve_apply (e, c, step, sys, &t, t1, &h, y);
      /* printf ("SYS  t = %.18e h = %g y0 = %g y1= %g\n", t, h, y[0], y[1]); */
    }

  frac = fabs ((y[1] - yfin[1]) / yfin[1]) + fabs ((y[0] - yfin[0]) / yfin[0]);

  if (frac > 2.0 * e->count * err_target)
    {
      printf ("SYS  t = %.5e h = %g y0 = %g y1= %g\n", t, h, y[0], y[1]);
      s++;
    }

  gsl_test (s, "%s, %s, evolve, standard control, relative error = %g",
	    gsl_odeiv_step_name (step), desc, frac);

  gsl_odeiv_evolve_free (e);
  gsl_odeiv_control_free (c);
  gsl_odeiv_step_free (step);
}

void
test_evolve_exp (const gsl_odeiv_step_type * T, double h, double err)
{
  double y[2];
  double yfin[2];

  y[0] = 1.0;
  y[1] = 1.0;
  yfin[0] = exp (10.0);
  yfin[1] = yfin[0];
  test_evolve_system (T, &rhs_func_exp, 0.0, 10.0, h, y, yfin, err, "exp [0,10]");
}

void
test_evolve_sin (const gsl_odeiv_step_type * T, double h, double err)
{
  double y[2];
  double yfin[2];

  y[0] = 1.0;
  y[1] = 0.0;
  yfin[0] = cos (2.0);
  yfin[1] = sin (2.0);
  test_evolve_system (T, &rhs_func_sin, 0.0, 2.0, h, y, yfin, err, "sine [0,2]");
}

void
test_evolve_stiff1 (const gsl_odeiv_step_type * T, double h, double err)
{
  double y[2];
  double yfin[2];

  y[0] = 1.0;
  y[1] = 0.0;
  {
    double arg = 1.0;
    double e1 = exp (-arg);
    double e2 = exp (-1000.0 * arg);
    yfin[0] = 2.0 * e1 - e2;
    yfin[1] = -e1 + e2;
  }
  test_evolve_system (T, &rhs_func_stiff, 0.0, 1.0, h, y, yfin, err,
                      "stiff [0,1]");
}

void
test_evolve_stiff5 (const gsl_odeiv_step_type * T, double h, double err)
{
  double y[2];
  double yfin[2];

  y[0] = 1.0;
  y[1] = 0.0;
  {
    double arg = 5.0;
    double e1 = exp (-arg);
    double e2 = exp (-1000.0 * arg);
    yfin[0] = 2.0 * e1 - e2;
    yfin[1] = -e1 + e2;
  }
  test_evolve_system (T, &rhs_func_stiff, 0.0, 5.0, h, y, yfin, err,
                      "stiff [0,5]");
}

int
main (void)
{
  int i;

  struct ptype
  {
    const gsl_odeiv_step_type *type;
    double h;
  }
  p[20];

  p[0].type = gsl_odeiv_step_rk2;
  p[0].h = 1.0e-03;
  p[1].type = gsl_odeiv_step_rk2imp;
  p[1].h = 1.0e-03;
  p[2].type = gsl_odeiv_step_rk4;
  p[2].h = 1.0e-03;
  p[3].type = gsl_odeiv_step_rk4imp;
  p[3].h = 1.0e-03;
  p[4].type = gsl_odeiv_step_rkf45;
  p[4].h = 1.0e-03;
  p[5].type = gsl_odeiv_step_rk8pd;
  p[5].h = 1.0e-03;
  p[6].type = gsl_odeiv_step_rkck;
  p[6].h = 1.0e-03;
  p[7].type = gsl_odeiv_step_bsimp;
  p[7].h = 0.1;
  p[8].type = gsl_odeiv_step_gear1;
  p[8].h = 1.0e-03;
  p[9].type = gsl_odeiv_step_gear2;
  p[9].h = 1.0e-03;
  p[10].type = 0;

  gsl_ieee_env_setup ();

  for (i = 0; p[i].type != 0; i++)
    {
      test_stepper_linear (p[i].type, p[i].h, GSL_SQRT_DBL_EPSILON);
      test_stepper_sin (p[i].type, p[i].h / 10.0, GSL_SQRT_DBL_EPSILON);
      test_stepper_exp (p[i].type, p[i].h / 10.0, GSL_SQRT_DBL_EPSILON);
      test_stepper_stiff (p[i].type, p[i].h / 10.0, GSL_SQRT_DBL_EPSILON);
    }

  for (i = 0; p[i].type != 0; i++)
    {
      test_evolve_exp (p[i].type, p[i].h, GSL_SQRT_DBL_EPSILON);
      test_evolve_sin (p[i].type, p[i].h, GSL_SQRT_DBL_EPSILON);
      test_evolve_stiff1 (p[i].type, p[i].h, 1e-5);
      test_evolve_stiff5 (p[i].type, p[i].h, 1e-5);
    }

  exit (gsl_test_summary ());
}


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

int
main ()
{
  int lda = 3;

  float A[] = { 0.11, 0.12, 0.13,
                0.21, 0.22, 0.23 };

  int ldb = 2;

  float B[] = { 1011, 1012,
                1021, 1022,
                1031, 1032 };

  int ldc = 2;

  float C[] = { 0.00, 0.00,
                0.00, 0.00 };

  /* Compute C = A B */

  cblas_sgemm (CblasRowMajor, CblasNoTrans, CblasNoTrans, 2, 2, 3,
               1.0, A, lda, B, ldb, 0.0, C, ldc);

  printf("[ %g, %g\n", C[0], C[1]);
  printf("  %g, %g ]\n", C[2], C[3]);

}

ext: c

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
#include <stdio.h>
#include <stdarg.h>
#include <string.h>
#include <math.h>
#include <gbpLib.h>
#include <gbpRNG.h>
#include <gbpMCMC.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_fit.h>
#include <gsl/gsl_interp.h>

void init_MCMC(MCMC_info * MCMC,
               const char *problem_name,
               void *      params,
               int (*f)(double *, MCMC_info *, double **),
               int     n_P,
               double *P_init,
               char ** P_names,
               double *P_limit_min,
               double *P_limit_max,
               int     n_arrays,
               ...) {
    int     i_P;
    int     i_array;
    int     i;
    FILE *  ft, *ft_restart;
    char    test_dir[256], test_restart[256];
    va_list vargs;
    va_start(vargs, n_arrays);

    SID_log("Initializing MCMC structure...", SID_LOG_OPEN);

    // Set defaults to bare minimums
    MCMC->n_avg                 = 100;
    MCMC->n_iterations_burn     = 4;
    MCMC->n_iterations          = 8;
    MCMC->n_thin                = 1;
    MCMC->coverage_size         = 100;
    MCMC->flag_autocor_on       = GBP_FALSE;
    MCMC->flag_integrate_on     = GBP_FALSE;
    MCMC->flag_analysis_on      = GBP_TRUE;
    MCMC->first_map_call        = GBP_TRUE;
    MCMC->first_link_call       = GBP_TRUE;
    MCMC->flag_init_chain       = GBP_TRUE;
    MCMC->first_chain_call      = GBP_TRUE;
    MCMC->first_parameter_call  = GBP_TRUE;
    MCMC->first_likelihood_call = GBP_TRUE;
    MCMC->ln_likelihood_last    = 0.;
    MCMC->ln_likelihood_new     = 0.;
    MCMC->ln_likelihood_chain   = 0.;
    MCMC->P_init                = NULL;
    MCMC->P_new                 = NULL;
    MCMC->P_last                = NULL;
    MCMC->P_chain               = NULL;
    MCMC->P_limit_min           = NULL;
    MCMC->P_limit_max           = NULL;
    MCMC->P_min                 = NULL;
    MCMC->P_max                 = NULL;
    MCMC->P_avg                 = NULL;
    MCMC->dP_avg                = NULL;
    MCMC->P_best                = NULL;
    MCMC->P_peak                = NULL;
    MCMC->P_lo_68               = NULL;
    MCMC->P_hi_68               = NULL;
    MCMC->P_lo_95               = NULL;
    MCMC->P_hi_95               = NULL;
    MCMC->n_M                   = NULL;
    MCMC->M_new                 = NULL;
    MCMC->M_last                = NULL;
    MCMC->DS                    = NULL;
    MCMC->last                  = NULL;
    MCMC->V                     = NULL;
    MCMC->m                     = NULL;
    MCMC->b                     = NULL;
    MCMC->RNG                   = NULL;
    MCMC->params                = NULL;
    MCMC->temperature           = 1.0;
    MCMC->n_DS                  = 0;
    MCMC->n_M_total             = 0;
    MCMC->n_fail                = 0;
    MCMC->n_success             = 0;
    MCMC->n_propositions        = 0;
    MCMC->n_map_calls           = 0;

    // Process the passed arguments
    MCMC->map_P_to_M                 = f;
    MCMC->compute_MCMC_ln_likelihood = compute_MCMC_ln_likelihood_default;
    MCMC->params                     = params;
    MCMC->n_P                        = n_P;
    sprintf(MCMC->problem_name, "%s", problem_name);
    MCMC->P_names       = (char **)SID_malloc(sizeof(char *) * MCMC->n_P);
    MCMC->P_name_length = 0;
    for(i_P = 0; i_P < n_P; i_P++) {
        MCMC->P_names[i_P] = (char *)SID_malloc(sizeof(char) * MCMC_NAME_SIZE);
        sprintf(MCMC->P_names[i_P], "%s", P_names[i_P]);
        MCMC->P_name_length = GBP_MAX((size_t)(MCMC->P_name_length), strlen(MCMC->P_names[i_P]));
    }
    sprintf(MCMC->P_name_format, "%%-%ds", MCMC->P_name_length);

    // Initialize the MCMC mode and set things associated with it
    set_MCMC_mode(MCMC, MCMC_MODE_DEFAULT); // MCMC->my_chain is set here

    // Set parameter arrays and limits
    MCMC->P_init      = (double *)SID_malloc(sizeof(double) * MCMC->n_P);
    MCMC->P_new       = (double *)SID_malloc(sizeof(double) * MCMC->n_P);
    MCMC->P_last      = (double *)SID_malloc(sizeof(double) * MCMC->n_P);
    MCMC->P_chain     = (double *)SID_malloc(sizeof(double) * MCMC->n_P);
    MCMC->P_limit_min = (double *)SID_malloc(sizeof(double) * MCMC->n_P);
    MCMC->P_limit_max = (double *)SID_malloc(sizeof(double) * MCMC->n_P);
    if(P_limit_min == NULL) {
        for(i_P = 0; i_P < n_P; i_P++)
            MCMC->P_limit_min[i_P] = -DBL_MAX * 1e-3;
    } else {
        for(i_P = 0; i_P < n_P; i_P++)
            MCMC->P_limit_min[i_P] = P_limit_min[i_P];
    }
    if(P_limit_max == NULL) {
        for(i_P = 0; i_P < n_P; i_P++)
            MCMC->P_limit_max[i_P] = DBL_MAX * 1e-3;
    } else {
        for(i_P = 0; i_P < n_P; i_P++)
            MCMC->P_limit_max[i_P] = P_limit_max[i_P];
    }

    // Set parameter initial values
    memcpy(MCMC->P_init, P_init, (size_t)MCMC->n_P * sizeof(double));
    memcpy(MCMC->P_new, P_init, (size_t)MCMC->n_P * sizeof(double));
    memcpy(MCMC->P_last, P_init, (size_t)MCMC->n_P * sizeof(double));
    memcpy(MCMC->P_chain, P_init, (size_t)MCMC->n_P * sizeof(double));

    // Set arrays
    MCMC->n_arrays = n_arrays;
    if(n_arrays > 0) {
        MCMC->array      = (double **)SID_malloc(sizeof(double *) * MCMC->n_arrays);
        MCMC->array_name = (char **)SID_malloc(sizeof(char *) * MCMC->n_arrays);
        for(i_array = 0; i_array < n_arrays; i_array++) {
            MCMC->array[i_array]      = (double *)SID_malloc(sizeof(double) * MCMC->n_P);
            MCMC->array_name[i_array] = (char *)SID_malloc(sizeof(char) * MCMC_NAME_SIZE);
            memcpy(MCMC->array[i_array], (double *)va_arg(vargs, double *), (size_t)(MCMC->n_P) * sizeof(double));
            sprintf(MCMC->array_name[i_array], "%s", (char *)va_arg(vargs, char *));
        }
    } else
        MCMC->array = NULL;

    // Set autotune defaults (if needed)
    if(SID_CHECK_BITFIELD_SWITCH(MCMC->mode, MCMC_MODE_AUTOTUNE))
        set_MCMC_autotune(MCMC, -1., -1., -1., -1, -1, -1, -1); // Negatives mean use defaults (set in gbpMCMC.h)

    // Initialize Communicator
    SID_Comm_init(&(MCMC->comm));
    SID_Comm_split(SID_COMM_WORLD, MCMC->my_chain, SID.My_rank, MCMC->comm);

    // Set the base output directory
    if(MCMC->my_chain == SID.My_rank) {
        // Generate directory and stop filenames
        sprintf(test_dir, "./%s_MCMC/", SID.My_binary);
        sprintf(test_restart, "./%s_MCMC/stop", SID.My_binary);

        // analyze_MCMC() needs to know what the filename root is.  We also need to strip trailing '/'s.
        strcpy(MCMC->filename_output_root, test_dir);
        int i_char;
        int flag_continue;
        for(i_char = strlen(MCMC->filename_output_root) - 1, flag_continue = GBP_TRUE; i_char >= 0 && flag_continue; i_char--) {
            if(MCMC->filename_output_root[i_char] != '/')
                flag_continue = GBP_FALSE;
            else
                MCMC->filename_output_root[i_char] = '\0';
        }

        i = 0;
        // Try to open them ...
        ft         = fopen(test_dir, "r");
        ft_restart = fopen(test_restart, "r");
        // If the directory exists & there is no stop file in it then...
        if((ft != NULL) && (ft_restart == NULL)) {
            // Increment the directory suffix until we don't find a directory,
            // whilst always ensuring we don't come across a stop file (which indicates
            // the previous run was stopped prematurely and we want restart the run).
            do {
                fclose(ft);
                i++;
                sprintf(test_dir, "./%s_MCMC.%d/", SID.My_binary, i);
                ft = fopen(test_dir, "r");
                sprintf(test_restart, "./%s_MCMC.%d/stop", SID.My_binary, i);
                ft_restart = fopen(test_restart, "r");
            } while((ft != NULL) && (ft_restart == NULL));
        }
        // Finally, if we did find a stop file then close it and remove it.
        if(ft_restart != NULL) {
            fclose(ft_restart);
            remove(test_restart);
        }
        // Copy the directory name we settled on to the MCMC structure
        strcpy(MCMC->filename_output_dir, test_dir);
        SID_log("ouput_dir set to: %s", SID_LOG_COMMENT, MCMC->filename_output_dir);
    }
    // Broadcast the output directory to all the other cores.
    SID_Bcast(MCMC->filename_output_dir, SID_MAX_FILENAME_LENGTH, SID_CHAR, MCMC->my_chain, MCMC->comm);

    // Initilize the random number generator
    MCMC->RNG = (RNG_info *)SID_malloc(sizeof(RNG_info));
    init_seed_from_clock(&(MCMC->seed));
    SID_Bcast(&(MCMC->seed), 1, SID_INT, SID_MASTER_RANK, MCMC->comm);
    init_RNG(&(MCMC->seed), MCMC->RNG, RNG_DEFAULT);

    SID_log("Done.", SID_LOG_CLOSE);

    va_end(vargs);
}

max_stars_repo_name: gbpoole/gbpCode

max_issues_repo_name: gbpoole/gbpCode

max_forks_repo_name: gbpoole/gbpCode


text:
#ifndef WILCOXONPARALLELTESTS_H
#define WILCOXONPARALLELTESTS_H
#include <fstream>
#include <cmath>
#include <iostream>
#include <sstream>
#include <cstdlib>
#include <vector>
#include <cstring>
#include <gsl/gsl_cdf.h>

using namespace std;

struct approximatePosition {
    int x;
    double y;
};

class WilcoxonTest {
public:
    WilcoxonTest(float * _data, int _dataXsize, int _dataYsize, string _testIndexes, string _controlIndexes);
    WilcoxonTest(float * _data, int _dataXsize, int _dataYsize, vector<int> * _testIndexes, vector<int> * _controlIndexes);
	vector<double> * test();
private:

    //Variables
    int dataYsize;
    int dataXsize;
    float * data;
    vector<double> * pValues;
    string netCdfFileName;
    vector<int> * testIndexes;
    vector<int> * controlIndexes;
    std::vector<std::vector<approximatePosition> * > * approximatePTable;

    //Set up
    void readNetCdfFile(string dataName, string xDimension, string yDimension);
    std::vector<approximatePosition> * getPositions(string positionsLine);
    std::vector<string> * splitLine(string inputString, char lineSplit = ';');
    void readApproximatePtable();
    vector<int> * parseIntString(string input);

    //Wilcoxon Test main methods
    float calculateWValue(int yIndex, vector<float> * absoluteValues, vector<float> * signs);
    double calculatePValue(float w, int numberOfZeroes);
    float calculateZValue(float w, int Nr);

    //Wilcoxon Test helper methods
    double getApproximatePValue(float w, float z);
    double approximateP(float w, float z, approximatePosition beginningPos, approximatePosition endPos);
    double * rankThePairs(int yIndex, vector<float> * absoluteValues);
    int getSign(float value);
    int getNumberOfZeroes(vector<float> * absoluteValues);
    void trim(string& str);

    //sorting
    void quicksort(int m, int n, vector<float> * absoluteValues, vector<float> * signs);
    void swap(float * x, float * y);
    int choose_pivot(int i, int j );
};

#endif

ext: h

max_stars_repo_name: stenver/wilcoxon-test

max_issues_repo_name: stenver/wilcoxon-test

max_forks_repo_name: stenver/wilcoxon-test


text:
/* integration/qcheb.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 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.
 */

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

/* This function computes the 12-th order and 24-th order Chebyshev
   approximations to f(x) on [a,b] */

void
gsl_integration_qcheb (gsl_function * f, double a, double b, double *cheb12, double *cheb24)
{
  size_t i;
  double fval[25], v[12];

  /* These are the values of cos(pi*k/24) for k=1..11 needed for the
     Chebyshev expansion of f(x) */

  const double x[11] = { 0.9914448613738104,
                         0.9659258262890683,
                         0.9238795325112868,
                         0.8660254037844386,
                         0.7933533402912352,
                         0.7071067811865475,
                         0.6087614290087206,
                         0.5000000000000000,
                         0.3826834323650898,
                         0.2588190451025208,
                         0.1305261922200516 };

  const double center = 0.5 * (b + a);
  const double half_length =  0.5 * (b - a);

  fval[0] = 0.5 * GSL_FN_EVAL (f, b);
  fval[12] = GSL_FN_EVAL (f, center);
  fval[24] = 0.5 * GSL_FN_EVAL (f, a);

  for (i = 1; i < 12; i++)
    {
      const size_t j = 24 - i;
      const double u = half_length * x[i-1];
      fval[i] = GSL_FN_EVAL(f, center + u);
      fval[j] = GSL_FN_EVAL(f, center - u);
    }

  for (i = 0; i < 12; i++)
    {
      const size_t j = 24 - i;
      v[i] = fval[i] - fval[j];
      fval[i] = fval[i] + fval[j];
    }

  {
    const double alam1 = v[0] - v[8];
    const double alam2 = x[5] * (v[2] - v[6] - v[10]);

    cheb12[3] = alam1 + alam2;
    cheb12[9] = alam1 - alam2;
  }

  {
    const double alam1 = v[1] - v[7] - v[9];
    const double alam2 = v[3] - v[5] - v[11];
    {
      const double alam = x[2] * alam1 + x[8] * alam2;

      cheb24[3] = cheb12[3] + alam;
      cheb24[21] = cheb12[3] - alam;
    }

    {
      const double alam = x[8] * alam1 - x[2] * alam2;
      cheb24[9] = cheb12[9] + alam;
      cheb24[15] = cheb12[9] - alam;
    }
  }

  {
    const double part1 = x[3] * v[4];
    const double part2 = x[7] * v[8];
    const double part3 = x[5] * v[6];

    {
      const double alam1 = v[0] + part1 + part2;
      const double alam2 = x[1] * v[2] + part3 + x[9] * v[10];

      cheb12[1] = alam1 + alam2;
      cheb12[11] = alam1 - alam2;
    }

    {
      const double alam1 = v[0] - part1 + part2;
      const double alam2 = x[9] * v[2] - part3 + x[1] * v[10];
      cheb12[5] = alam1 + alam2;
      cheb12[7] = alam1 - alam2;
    }
  }

  {
    const double alam = (x[0] * v[1] + x[2] * v[3] + x[4] * v[5]
                + x[6] * v[7] + x[8] * v[9] + x[10] * v[11]);
    cheb24[1] = cheb12[1] + alam;
    cheb24[23] = cheb12[1] - alam;
  }

  {
    const double alam = (x[10] * v[1] - x[8] * v[3] + x[6] * v[5]
                - x[4] * v[7] + x[2] * v[9] - x[0] * v[11]);
    cheb24[11] = cheb12[11] + alam;
    cheb24[13] = cheb12[11] - alam;
  }

  {
    const double alam = (x[4] * v[1] - x[8] * v[3] - x[0] * v[5]
                - x[10] * v[7] + x[2] * v[9] + x[6] * v[11]);
    cheb24[5] = cheb12[5] + alam;
    cheb24[19] = cheb12[5] - alam;
  }

  {
    const double alam = (x[6] * v[1] - x[2] * v[3] - x[10] * v[5]
                + x[0] * v[7] - x[8] * v[9] - x[4] * v[11]);
    cheb24[7] = cheb12[7] + alam;
    cheb24[17] = cheb12[7] - alam;
  }

  for (i = 0; i < 6; i++)
    {
      const size_t j = 12 - i;
      v[i] = fval[i] - fval[j];
      fval[i] = fval[i] + fval[j];
    }

  {
    const double alam1 = v[0] + x[7] * v[4];
    const double alam2 = x[3] * v[2];

    cheb12[2] = alam1 + alam2;
    cheb12[10] = alam1 - alam2;
  }

  cheb12[6] = v[0] - v[4];

  {
    const double alam = x[1] * v[1] + x[5] * v[3] + x[9] * v[5];
    cheb24[2] = cheb12[2] + alam;
    cheb24[22] = cheb12[2] - alam;
  }

  {
    const double alam = x[5] * (v[1] - v[3] - v[5]);
    cheb24[6] = cheb12[6] + alam;
    cheb24[18] = cheb12[6] - alam;
  }

  {
    const double alam = x[9] * v[1] - x[5] * v[3] + x[1] * v[5];
    cheb24[10] = cheb12[10] + alam;
    cheb24[14] = cheb12[10] - alam;
  }

  for (i = 0; i < 3; i++)
    {
      const size_t j = 6 - i;
      v[i] = fval[i] - fval[j];
      fval[i] = fval[i] + fval[j];
    }

  cheb12[4] = v[0] + x[7] * v[2];
  cheb12[8] = fval[0] - x[7] * fval[2];

  {
    const double alam = x[3] * v[1];
    cheb24[4] = cheb12[4] + alam;
    cheb24[20] = cheb12[4] - alam;
  }

  {
    const double alam = x[7] * fval[1] - fval[3];
    cheb24[8] = cheb12[8] + alam;
    cheb24[16] = cheb12[8] - alam;
  }

  cheb12[0] = fval[0] + fval[2];

  {
    const double alam = fval[1] + fval[3];
    cheb24[0] = cheb12[0] + alam;
    cheb24[24] = cheb12[0] - alam;
  }

  cheb12[12] = v[0] - v[2];
  cheb24[12] = cheb12[12];

  for (i = 1; i < 12; i++)
    {
      cheb12[i] *= 1.0 / 6.0;
    }

  cheb12[0] *= 1.0 / 12.0;
  cheb12[12] *= 1.0 / 12.0;

  for (i = 1; i < 24; i++)
    {
      cheb24[i] *= 1.0 / 12.0;
    }

  cheb24[0] *= 1.0 / 24.0;
  cheb24[24] *= 1.0 / 24.0;
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
//----------------------------------------------------------------------------
//  History:    2018-04-30 Dwayne Robinson - Created
//----------------------------------------------------------------------------
#pragma once

#include <iterator>
#include <memory> // For uninitialized_move/copy and std::unique_ptr.
#include <assert.h>
#include <algorithm>

#if USE_CPP_MODULES
//import std.core;
import Common.ArrayRef;
#else
#ifdef USE_GSL_SPAN_INSTEAD_OF_ARRAY_REF
#include <gsl/span>
#define array_ref gsl::span
#else
#include "Common.ArrayRef.h" // gsl::span may be mostly substituted instead of array_ref, just missing intersects().
#endif
#endif

MODULE(Common.FastVector);
EXPORT_BEGIN

#pragma warning(push)
#pragma warning(disable:4127) // Conditional expression is constant. VS can't tell that certain compound conditionals of template parameters aren't always constant when the tempalate parameter is true.

#if defined(_ITERATOR_DEBUG_LEVEL) && _ITERATOR_DEBUG_LEVEL > 0
// For std::uninitialized_copy and std::uninitialized_move.
#define FASTVECTOR_MAKE_UNCHECKED stdext::make_unchecked_array_iterator
#else
#define FASTVECTOR_MAKE_UNCHECKED
#endif

// fast_vector is a dynamic array that can substitute for std::vector, where it:
// (1) avoids heap allocations when the element count fits within the fixed-size capacity.
// (2) avoids unnecessarily initializing elements if ShouldInitializeElements == false.
//     This is useful for large buffers which will just be overwritten soon anyway.
// (3) supports most vector methods except insert/erase/emplace.
//
// Template parameters:
// - DefaultArraySize - passing 0 means it is solely heap allocated. Passing > 0
//   reserves that much element capacity before allocating heap memory.
//
// - ShouldInitializeElements - ensures elements are constructed when resizing, and
//   it must be true for objects with non-trivial constructors, but it can be set
//   false for large buffers to avoid unnecessarily initializing memory which will
//   just be overwritten soon later anyway.
//
// Examples:
//  fast_vector<int, 20> axes        - up to 20 integers before heap allocation.
//  fast_vector<int, 0, false> axes  - always heap allocated but never initialized.
//  fast_vector<int> axes            - basically std::vector.

template<typename T, size_t DefaultArraySize = 0, bool ShouldInitializeElements = true>
class fast_vector;

enum fast_vector_use_memory_buffer_enum
{
    fast_vector_use_memory_buffer // To avoid ambiguous constructor overload resolution.
};

// The base class is separated out for the customization of passing specific
// memory, and to avoid bloating additional template permutations solely due to
// differing array sizes.
template<typename T, bool ShouldInitializeElements>
class fast_vector<T, 0, ShouldInitializeElements>
{
    static_assert(ShouldInitializeElements || std::is_trivial<T>::value);
    using self = fast_vector<T, 0, ShouldInitializeElements>;

public:
    // Standard container type definitions.
    using value_type = T;
    using pointer = T*;
    using reference = T&;
    using iterator = pointer;
    using const_reference = T const&;
    using const_iterator = T const*;
    using reverse_iterator = std::reverse_iterator<iterator>;
    using const_reverse_iterator = std::reverse_iterator<const_iterator>;
    using size_type = size_t;
    using difference_type = ptrdiff_t;
    using mutable_value_type = typename std::remove_const<T>::type;
    using mutable_iterator = mutable_value_type*;

    // If true, the data types needs additional alignment like SSE and AVX types.
    constexpr static bool NeedsTypeAlignmentBeyondMaxAlignT = alignof(T) > alignof(std::max_align_t);
    constexpr static size_t MinimumMemoryBlockSize = NeedsTypeAlignmentBeyondMaxAlignT ? sizeof(void*) : 0;

public:
    constexpr fast_vector() noexcept
    {
    }

    fast_vector(size_t initialSize)
    {
        resize(initialSize);
    }

    fast_vector(array_ref<const T> initialValues)
    {
        assign(initialValues);
    }

    fast_vector(const fast_vector& otherVector)
    {
        assign(otherVector);
    }

    fast_vector(fast_vector&& otherVector) // Throws std::bad_alloc if low on memory.
    {
        transfer_from(otherVector);
    }

    // Construct using explicit fixed size memory buffer.
    // - Used by the derived template specialization which includes a fixed size buffer.
    // - May also be used by the caller to pass an explicit buffer such as a local array,
    //   where the buffer is guaranteed to live for the lifetime of the fast_vector.
    constexpr fast_vector(fast_vector_use_memory_buffer_enum, array_ref<T> initialBackingArray)
    :   data_(initialBackingArray.data()),
        capacity_(initialBackingArray.size())
    {
    }

    fast_vector(fast_vector_use_memory_buffer_enum, array_ref<T> initialBackingArray, array_ref<const T> initialValues)
    :   data_(initialBackingArray.data()),
        capacity_(initialBackingArray.size())
    {
        assign(initialValues);
    }

    fast_vector(fast_vector_use_memory_buffer_enum, array_ref<T> initialBackingArray, size_t initialSize)
    :   data_(initialBackingArray.data()),
        capacity_(initialBackingArray.size())
    {
        resize(initialSize);
    }

    fast_vector(fast_vector_use_memory_buffer_enum, array_ref<T> initialBackingArray, pointer p, size_t elementCount)
    :   data_(initialBackingArray.data()),
        capacity_(initialBackingArray.size())
    {
        assign({p, elementCount});
    }

    template <typename IteratorType>
    fast_vector(fast_vector_use_memory_buffer_enum, array_ref<T> initialBackingArray, IteratorType begin, IteratorType end)
    :   data_(initialBackingArray.data()),
        capacity_(initialBackingArray.size())
    {
        assign(begin, end);
    }

    ~fast_vector() noexcept(std::is_nothrow_destructible<T>::value)
    {
        free();
    }

    // Clear the vector and free all memory. Calling clear() then shrink_to_fit()
    // is a less efficient way to accomplish it too.
    void free()
    {
        clear(); // Destruct objects and zero the size.

        if (dataIsAllocatedMemory_)
        {
            FreeMemoryBlock(data_);
            data_ = nullptr;
            capacity_ = 0;
        }
    }

    // Iterators
    iterator begin() const noexcept                 { return data_; }
    iterator end() const noexcept                   { return data_ + size_; }
    const_iterator cbegin() const noexcept          { return begin(); }
    const_iterator cend() const noexcept            { return end(); };
    reverse_iterator rbegin() const noexcept        { return reverse_iterator(begin()); }
    reverse_iterator rend() const noexcept          { return reverse_iterator(end()); }
    const_reverse_iterator crbegin() const noexcept { return const_reverse_iterator(begin()); }
    const_reverse_iterator crend() const noexcept   { return const_reverse_iterator(end()); }

    // Capacity
    size_type size() const noexcept                 { return size_; }
    size_type size_in_bytes() const noexcept        { return size_ * sizeof(T); }
    size_type capacity() const noexcept             { return capacity_; }
    static constexpr size_type max_size() noexcept  { return SIZE_MAX / sizeof(T); }
    bool empty() const noexcept                     { return size_ == 0; }

    // Element access
    T& operator[](size_t i) const noexcept          { return data_[i]; }
    T& front() const noexcept                       { return data_[0]; }
    T& back() const noexcept                        { return data_[size_ - 1]; }
    T* data() const noexcept                        { return data_; }
    T* data_end() const noexcept                    { return data_ + size_; }

    array_ref<T> data_span() noexcept               { return {data_, size_}; }
    array_ref<T const> data_span() const noexcept   { return {data_, size_}; }

    T& at(size_t i)                                 { return checked_read(i); }
    T& at(size_t i) const                           { return checked_read(i); }

    T& checked_read(size_t i)
    {
        if (i >= size_)
        {
            throw std::out_of_range("fast_vector array index out of range.");
        }

        return data_[i];
    }

    const T& checked_read(size_t i) const
    {
        return const_cast<self&>(*this).checked_read(i);
    }

    fast_vector& operator=(const fast_vector& otherVector)
    {
        assign(otherVector);
        return *this;
    }

    fast_vector& operator=(fast_vector&& otherVector) // Throws std::bad_alloc if low on memory.
    {
        transfer_from(otherVector);
        return *this;
    }

    // Clear any existing elements and copy the new elements from span.
    void assign(array_ref<const T> span)
    {
        #ifndef USE_GSL_SPAN_INSTEAD_OF_ARRAY_REF
        assert(!span.intersects(data_span())); // No self intersection.
        #endif

        clear();

        size_t newSize = span.size();
        reserve(newSize);
        std::uninitialized_copy(FASTVECTOR_MAKE_UNCHECKED(span.data()), FASTVECTOR_MAKE_UNCHECKED(span.data()) + newSize, /*out*/ FASTVECTOR_MAKE_UNCHECKED(data_));
        size_ = newSize;
    }

    // Clear any existing elements and copy the new elements from the iterable range.
    template <typename IteratorType>
    void assign(IteratorType begin, IteratorType end)
    {
        clear();

        size_t newSize = std::distance(begin, end);
        reserve(newSize);
        std::uninitialized_copy(FASTVECTOR_MAKE_UNCHECKED(begin), FASTVECTOR_MAKE_UNCHECKED(end), /*out*/ FASTVECTOR_MAKE_UNCHECKED(data_));
        size_ = newSize;
    }

    // Tranfer all elements from the other vector to this one.
    // This only offers the weak exception guarantee, that no leaks will happen
    // if type T throws in the middle of a copy.
    void transfer_from(array_ref<const T> span)
    {
        #ifndef USE_GSL_SPAN_INSTEAD_OF_ARRAY_REF
        assert(!span.intersects(data_span())); // No self intersection.
        #endif

        clear();

        size_t newSize = span.size();
        reserve(newSize);
        std::uninitialized_move(FASTVECTOR_MAKE_UNCHECKED(span.data()), FASTVECTOR_MAKE_UNCHECKED(span.data()) + newSize, /*out*/ FASTVECTOR_MAKE_UNCHECKED(data_));
        size_ = newSize;
    }

    // Tranfer all elements from the other vector to this one.
    // This only offers the weak exception guarantee, that no leaks will happen
    // if type T throws in the middle of a copy.
    void transfer_from(fast_vector<T, 0, ShouldInitializeElements>& other)
    {
        if (&other == this)
        {
            return; // Nop for self assignment.
        }

        if (other.dataIsAllocatedMemory_)
        {
            free();

            // Take ownership of new data directly.
            data_ = other.data_;
            other.data_ = nullptr;
            size_ = other.size_;
            other.size_ = 0;
            capacity_ = other.capacity_;
            other.capacity_ = 0;
            dataIsAllocatedMemory_ = other.dataIsAllocatedMemory_;
            other.dataIsAllocatedMemory_ = false;
        }
        else
        {
            // Copying from a fixed size buffer; so it's unsafe to simply
            // steal the pointers as that may leave a dangling pointer
            // when the other fast vector disappears.
            transfer_from(make_array_ref(other));
        }
    }

    void clear() noexcept(std::is_nothrow_destructible<T>::value)
    {
        if (ShouldInitializeElements)
        {
            std::destroy(data_, data_ + size_);
        }
        size_ = 0;
        // But do not free heap memory.
    }

    void resize(size_t newSize)
    {
        if (newSize > size_)
        {
            if (newSize > capacity_)
            {
                reserve_at_least(newSize);
            }

            // Grow data to the new size, calling the default constructor on each new item.
            if (ShouldInitializeElements)
            {
                std::uninitialized_value_construct<iterator>(data_ + size_, data_ + newSize);
            }

            size_ = newSize;
        }
        else if (newSize < size_)
        {
            // Shrink the data to the new size, calling the destructor on each item. Capacity remains intact.
            if (ShouldInitializeElements)
            {
                std::destroy(data_ + newSize, data_ + size_);
            }

            size_ = newSize;
        }
    }

    void reserve(size_t newCapacity)
    {
        if (newCapacity <= capacity_)
        {
            return; // Nothing to do.
        }

        if (newCapacity > max_size())
        {
            throw std::bad_alloc(); // Too many elements.
        }

        size_t newByteSize = newCapacity * sizeof(T);
        ReallocateMemory(newByteSize);

        capacity_ = newCapacity;
    }

    void reserve_at_least(size_t newCapacity)
    {
        // Grow with 1.5x factor to avoid frequent reallocations.
        newCapacity = std::max((size_ * 3 / 2), newCapacity);
        reserve(newCapacity);
    }

    void shrink_to_fit()
    {
        if (!dataIsAllocatedMemory_ || capacity_ == size_)
        {
            return; // Nothing to do.
        }

        size_t newByteSize = size_ * sizeof(T);
        ReallocateMemory(newByteSize);

        capacity_ = size_;
    }

    void push_back(const T& newValue)
    {
        reserve(size_ + 1);//
        new(static_cast<void*>(data_ + size_)) T(newValue);
        ++size_;
    }

    void push_back(T&& newValue)
    {
        reserve(size_ + 1);
        new(static_cast<void*>(data_ + size_)) T(std::move(newValue));
        ++size_;
    }

    void append(array_ref<const T> span)
    {
        insert(size_, span);
    }

    void insert(size_t insertionOffset, array_ref<const T> span)
    {
        insert(insertionOffset, span.begin(), span.end());
    }

    void insert(size_t insertionOffset, const_iterator begin, const_iterator end)
    {
        assert(insertionOffset <= size_);

        // Grow array size.
        size_t newSize = std::distance(begin, end) + size_;
        if (newSize < size_)
        {
            throw std::bad_alloc();
        }
        reserve_at_least(newSize);

        // Add empty elements to the end.
        size_t oldSize = size_;
        if (ShouldInitializeElements)
        {
            std::uninitialized_value_construct<iterator>(data_ + size_, data_ + newSize);
        }
        size_ = newSize;

        // Shift elements to the end to make more room.
        std::move_backward(data_ + insertionOffset, data_ + oldSize, data_ + newSize);

        std::copy(begin, end, /*out*/ data_ + insertionOffset);
    }

    void insert(const_iterator position, array_ref<const T> span)
    {
        insert(position - begin(), span.begin(), span.end());
    }

    void erase(const_iterator begin, const_iterator end)
    {
        assert(end >= begin);
        assert(begin >= data_);
        assert(end <= data_ + size_);

        // Shift elements to front.
        size_t oldSize = size_;
        std::copy(iterator(end), data_ + size_, iterator(begin));
        size_ -= std::distance(begin, end);

        // Destroy elements at the end.
        if (ShouldInitializeElements)
        {
            std::destroy(data_ + size_, data_ + oldSize);
        }
    }

    void erase(const_iterator position)
    {
        erase(position, position + 1);
    }

    // Returns underlying malloc()'d memory block.
    // - May be transferred to another fast_vector via attach_memory.
    // - Caller must free() the memory if not transferred.
    // - This is more safely used when T is a POD type, as the caller would also
    //   need to appropriately call any destructors if complex objects.
    // - If the vector is using fixed size memory (no allocations have happened),
    //   the returned span is empty and points to null.
    // - If the vector is using allocated memory, the memory block will be not
    //   empty and non-null, and the vector is then empty.
    // - No object destructors are called.
    // - If NeedsTypeAlignmentBeyondMaxAlignT is true, the data array might start
    //   after the memory block for alignment.
    //
    // fast_vector<int> vector(20);
    // ...
    // auto memory = vector.detach_memory();
    // std::unique_ptr<char, decltype(&std::free)> p(memory.data(), &std::free);
    // ...
    // otherVector.attach_memory(memory);
    // p.release(); // unique_ptr does not own it anymore.
    //
    array_ref<uint8_t> detach_memory() noexcept
    {
        array_ref<uint8_t> data;

        // Only heap allocated memory can be returned, not the fixed size buffer.
        if (dataIsAllocatedMemory_)
        {
            // Return the raw memory block, from the actual beginning of memory to the true end.
            // The data might start beyond the memory block for alignment purposes.
            uint8_t* bytes = reinterpret_cast<uint8_t*>(data_);
            size_t byteCount = size_ * sizeof(T);
            data = {reinterpret_cast<uint8_t*>(GetMemoryBlock(data_)), bytes + byteCount};

            data_ = nullptr;
            size_ = 0;
            capacity_ = 0;
        }

        return data;
    }

    // Take ownership of the memory, which came from malloc or another fast_vector
    // via detach_memory.
    // - If ShouldInitializeElements == true, then the memory is presumed to contain
    //   valid objects, which will be destructed when the fast_vector dies.
    // - Because the memory is raw, it's possible to reuse a different data type.
    // - The memory must be aligned to alignof(std::max_align_t).
    // - If NeedsTypeAlignmentBeyondMaxAlignT is true, the memory must be at least
    //   sizeof(void*) bytes.
    //
    void attach_memory(array_ref<uint8_t> memory) noexcept(std::is_nothrow_destructible<T>::value)
    {
        #ifndef USE_GSL_SPAN_INSTEAD_OF_ARRAY_REF
        assert(!memory.intersects(data_span())); // No self intersection.
        #endif

        assert((reinterpret_cast<size_t>(memory.data()) & (alignof(std::max_align_t) - 1)) == 0); // malloc/realloc should return aligned pointers up to max_align_t.
        assert(memory.data() == nullptr || memory.size() >= MinimumMemoryBlockSize);

        free();

        T* data = AlignMemoryBlock(memory.data());

        // Take ownership of new data.
        data_ = data;
        uint8_t* dataEnd = memory.data() + memory.size();
        size_ = (dataEnd - reinterpret_cast<uint8_t*>(data_)) / sizeof(T);
        capacity_ = size_;
        dataIsAllocatedMemory_ = true;
    }

protected:
    void ReallocateMemory(size_t newByteSize) // Throws std::bad_alloc if low on memory, or if T's move constructor fails.
    {
        assert(newByteSize >= size_ * sizeof(T)); // Shouldn't have been called otherwise, because it's wrong for size_ to be less than actual memory.

        // Handle alignment needs when the type is beyond what malloc/realloc guarantee.
        // This needed for SSE and AVX types (which need 16 and 32 bytes). Any other standard
        // type should have sufficient alignment by default.
        newByteSize = InflateMemorySizeForAlignment(newByteSize);

        // Try to reallocate the memory block directly rather allocate a new block and manually
        // copying. The memory manager can sometimes just extend the existing data block in-place
        // if there is space behind the block. Only trivially moveable types comply though, whereas
        // std::string would not comply due to the small string optimization.

        if (std::is_trivially_move_constructible<T>::value && dataIsAllocatedMemory_)
        {
            // Try to just reallocate the existing memory block.

            void* memory = realloc(GetMemoryBlock(data_), newByteSize);
            if (memory == nullptr && newByteSize > 0)
            {
                throw std::bad_alloc();
            }
            data_ = AlignMemoryBlock(/*modified*/ memory);
        }
        else
        {
            // Allocate a new memory buffer if one isn't allocated yet,
            // or if the data type is complex enough that it's not trivially
            // moveable (e.g. std::string, which contains pointers that point
            // into the class address itself for the small string optimization).

            void* memory = std::malloc(newByteSize);
            std::unique_ptr<void, decltype(std::free)*> newDataHolder(memory, &std::free);

            // Copy an any existing elements from the fixed size buffer.
            T* newData = AlignMemoryBlock(memory);
            std::uninitialized_move(FASTVECTOR_MAKE_UNCHECKED(data_), FASTVECTOR_MAKE_UNCHECKED(data_) + size_, /*out*/ FASTVECTOR_MAKE_UNCHECKED(newData));

            // Release the existing block, and assign the new one.
            if (dataIsAllocatedMemory_)
            {
                FreeMemoryBlock(data_);
            }
            newDataHolder.release();
            data_ = newData;
            dataIsAllocatedMemory_ = true;

            // Caller updates capacity_ and size_. This simply reallocates the memory block.
        }
    }

    // No object destruction, just a direct free.
    static void FreeMemoryBlock(T* data)
    {
        std::free(GetMemoryBlock(data));
    }

    // Get the actual memory block associated with the data.
    static void* GetMemoryBlock(T* data)
    {
        mutable_value_type* mutableData = const_cast<mutable_value_type*>(data);

        if (NeedsTypeAlignmentBeyondMaxAlignT && data != nullptr)
        {
            // Read the pointer to the actual block which is set immediately before the data.
            // [ptr to memory block][data ...]
            return reinterpret_cast<void**>(mutableData)[-1];
        }
        else
        {
            return mutableData;
        }
    }

    // Return aligned data for the given type, storing the raw memory block from malloc.
    // Write the memory block associated with the data, and return the data.
    // The caller has already ensure enough space exists in the block for alignment.
    static T* AlignMemoryBlock(void* memory)
    {
        if (NeedsTypeAlignmentBeyondMaxAlignT && memory != nullptr)
        {
            // Return a pointer to the aligned data, which starts after the pointer to the original
            // memory.
            //
            // [ptr to memory block][data ...]

            // Return  the pointer to the actual block which is set immediately before the data.
            // [ptr to memory block][data ...]
            constexpr size_t alignmentMask = alignof(T) - 1;
            static_assert((alignof(T) & alignmentMask) == 0, "Alignment must be a power of 2.");
            static_assert(NeedsTypeAlignmentBeyondMaxAlignT == false || alignof(T) >= alignof(void*), "Alignment must be at least the size of a pointer for NeedsTypeAlignmentBeyondMaxAlignT to be true.");

            // Align memory. Note calling std::align is awkward given it expects sizes, which are
            // already correct earlier and do not need to be passed to this function.
            size_t memoryAddress = reinterpret_cast<size_t>(memory);
            memoryAddress = (memoryAddress + sizeof(void*) + alignmentMask) & ~alignmentMask;
            mutable_value_type* data = reinterpret_cast<mutable_value_type*>(memoryAddress);

            // Store the pointer to the original memory block.
            reinterpret_cast<void**>(data)[-1] = memory;
            return data;
        }
        else
        {
            assert((reinterpret_cast<size_t>(memory) & (alignof(std::max_align_t) - 1)) == 0); // malloc/realloc should return aligned pointers up to max_align_t.
            return reinterpret_cast<T*>(memory);
        }
    }

    static size_t InflateMemorySizeForAlignment(size_t newByteSize)
    {
        if (NeedsTypeAlignmentBeyondMaxAlignT && newByteSize != 0)
        {
            // Need to increase the allocation size to include enough for full type alignment
            // and the size of a single pointer preceding the data.
            size_t alignedByteSize = newByteSize + alignof(T) - alignof(std::max_align_t) + MinimumMemoryBlockSize;
            if (alignedByteSize < newByteSize)
            {
                throw std::bad_alloc(); // Too much memory requested and overflowed.
            }
            newByteSize = alignedByteSize;
        }
        // else malloc/realloc default alignment will suffice.

        return newByteSize;
    }

protected:
    T* data_ = nullptr;     // May point to fixed size array or allocated memory, depending on dataIsAllocatedMemory_.
    size_t size_ = 0;       // Count in elements.
    size_t capacity_ = 0;   // Count in elements.
    bool dataIsAllocatedMemory_ = false;
};

// Derived specialization which includes a fixed size buffer.
template <typename T, size_t DefaultArraySize, bool ShouldInitializeElements>
class fast_vector : public fast_vector<T, 0, ShouldInitializeElements>
{
public:
    using BaseClass = fast_vector<T, 0, ShouldInitializeElements>;

    constexpr fast_vector() noexcept
    :   BaseClass(fast_vector_use_memory_buffer, GetFixedSizeArrayData())
    {
    }

    fast_vector(array_ref<const T> initialValues)
    :   BaseClass(fast_vector_use_memory_buffer, GetFixedSizeArrayData(), initialValues)
    {
    }

    fast_vector(std::initializer_list<const T> initialValues)
    :   BaseClass(fast_vector_use_memory_buffer, GetFixedSizeArrayData(), make_array_ref(initialValues))
    {
    }

    fast_vector(const fast_vector& otherVector)
    :   BaseClass(otherVector)
    {
    }

    fast_vector(typename BaseClass::pointer p, size_t elementCount)
    :   BaseClass(fast_vector_use_memory_buffer, GetFixedSizeArrayData(), p, elementCount)
    {
    }

    template <typename IteratorType>
    fast_vector(IteratorType begin, IteratorType end)
    :   BaseClass(fast_vector_use_memory_buffer, GetFixedSizeArrayData(), begin, end)
    {
    }

    fast_vector(size_t initialSize)
    :   BaseClass(fast_vector_use_memory_buffer, GetFixedSizeArrayData(), initialSize)
    {
    }

    // Move constructor. Will not throw unless type T throws, since the
    // target is guaranteed to have enough space.
    fast_vector(fast_vector&& otherVector) noexcept(std::is_nothrow_move_assignable<T>::value)
    :   BaseClass(std::move(otherVector))
    {
    }

    // Explicit overload for base class, to avoid the compiler from undesireably
    // choosing the constructor overload which takes a span of initialValues or
    // the copy constructor.
    fast_vector(BaseClass&& otherVector) // Throws std::bad_alloc if low on memory.
    :   BaseClass(std::move(otherVector))
    {
    }

    // Move assignment. Will not throw unless type T throws, since the
    // target vector is guaranteed to have a large enough fixed-size array.
    fast_vector& operator=(fast_vector&& otherVector) noexcept(std::is_nothrow_move_assignable<T>::value && std::is_nothrow_destructible<T>::value)
    {
        BaseClass::transfer_from(otherVector);
        return *this;
    }

    fast_vector& operator=(BaseClass&& otherVector) // Throws std::bad_alloc if low on memory.
    {
        BaseClass::transfer_from(otherVector);
        return *this;
    }

    fast_vector& operator=(const fast_vector& otherVector)
    {
        BaseClass::assign(otherVector);
        return *this;
    }

    fast_vector& operator=(const BaseClass& otherVector)
    {
        BaseClass::assign(otherVector);
        return *this;
    }

private:
    constexpr array_ref<T> GetFixedSizeArrayData() noexcept
    {
        return array_ref<T>(reinterpret_cast<T*>(std::data(fixedSizedArrayData_)), std::size(fixedSizedArrayData_));
    }

private:
    // Uninitialized data to be used by the base class.
    // It's declared as raw bytes rather than an std::array<T> to avoid any initialization cost
    // up front, only initializing the fields which actually exist when resized later.
    //
    // Note Visual Studio 2017 15.8 requires you to declare _ENABLE_EXTENDED_ALIGNED_STORAGE
    // if your type T is > max_align_t.
    //
    std::aligned_storage_t<sizeof(T), alignof(T)> fixedSizedArrayData_[DefaultArraySize];
};

#ifdef USE_GSL_SPAN_INSTEAD_OF_ARRAY_REF
#undef array_ref
#endif

#undef FASTVECTOR_MAKE_UNCHECKED

#pragma warning(pop)

EXPORT_END

max_stars_repo_name: moyogo/TextLayoutSampler

max_issues_repo_name: moyogo/TextLayoutSampler

max_forks_repo_name: moyogo/TextLayoutSampler


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <assert.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>
#include <stdbool.h>

#include "core_allvars.h"
#include "core_proto.h"

#define IMF_CONSTANT IMF_norm/(IMF_slope + 1.0)

// Local Proto-Types //

void  update_SN_stars_array(int p, double stars, double dt, int tree, int ngal);
void update_stellar_tracking(int p, double stars, double dt, int tree, int ngal);

// External Functions //

void update_from_SN_feedback(int p, int centralgal, double reheated_mass, double ejected_mass, double mass_stars_recycled, double mass_metals_new, double NSN, double dt)
{

  double metallicity, metallicityHot, dust_fraction_cold, dust_fraction_hot;
  double dust_fraction_coldhot, dust_fraction_hotcold; // These are the dust fractions defined by M_colddust / (Mcold + Mhot) and M_hotdust / (Mcold + Mhot).
  // These second dust fractions are required because when we annihilate dust from SN we want the TOTAL amount to be NSN * eta * Ms.

  // Dust constants, taken mainly from Dayal, Ferrara and Saro 2011.
  double yd = 0.5 / 1.0e10 * Hubble_h; // Amount of dust mass created from each supernova event. Units of Msun.
  double eta = 0.4; // Coupling between shock energy and annihilation of dust.
  double Ms = 6.8e3 / 1.0e10 * Hubble_h; // Amount of mass that is accelerated to 100km/s by SN blastwave. Units of Msun.
  double dust_added = 0.0, dust_removed = 0.0;

  XASSERT((reheated_mass >= 0.0) && (ejected_mass >= 0.0) && (mass_stars_recycled >= 0.0) && (mass_metals_new >= 0.0), "When trying to update from the supernova feedback we got negative masses.\nReheated Mass = %.4e (1.0e10 Msun/h) \t Ejected Mass = %.4e (1.0e10 Msun/h) \t Mass of Stars Recycled = %.4e (1.0e10 Msun/h) \t Mass of Metals Added = %.4e (1.0e10 Msun/h)\n", reheated_mass, ejected_mass, mass_stars_recycled, mass_metals_new);

  Gal[p].StellarMass -= mass_stars_recycled; // The supernova remnants are instantly moved back into the cold ISM.
  Gal[p].Total_Stellar_Stars -= mass_stars_recycled;

  if (Gal[p].Total_Stellar_Stars < 0.0)
    Gal[p].Total_Stellar_Stars = 0.0;

  Gal[p].ColdGas += mass_stars_recycled;

  Gal[p].DustColdGas += yd * NSN; // Dust is created by supernova (something) on cold gas.
  Gal[p].ColdGas -= yd * NSN; // Hence the created dust is subtracted from the cold gas reservoir.

  dust_added += yd * NSN;

  if(Gal[p].ColdGas > 1.0e-10)
      Gal[p].MetalsColdGas += mass_metals_new; // ISM is enriched by new metals.
  else
      Gal[centralgal].MetalsHotGas += mass_metals_new;

  // Here we do checks to see if we need to rescale the amount of mass that has been ejected/reheated. //

  if(reheated_mass > Gal[p].ColdGas) // Just perform this check again to be sure.
      reheated_mass = Gal[p].ColdGas; // Because SF takes gas out of cold resevoir.

  if(ejected_mass > Gal[centralgal].HotGas)
      ejected_mass = Gal[centralgal].HotGas;

  metallicity = get_metallicity(Gal[p].ColdGas, Gal[p].MetalsColdGas);
  dust_fraction_cold = get_dust_fraction(Gal[p].ColdGas, Gal[p].DustColdGas);
  dust_fraction_coldhot = get_dust_fraction(Gal[p].ColdGas + Gal[p].HotGas + Gal[p].EjectedMass, Gal[p].DustColdGas);

  Gal[p].ColdGas -= reheated_mass;
  Gal[p].MetalsColdGas -= reheated_mass * metallicity;
  Gal[p].DustColdGas -= reheated_mass * dust_fraction_cold; // Dust that has been reheated to the hot reservoir.
  Gal[p].DustColdGas -= NSN * eta * Ms * dust_fraction_coldhot; // Dust that has been annihilated by shocks.

  dust_removed += NSN * eta * Ms * dust_fraction_coldhot;

  Gal[centralgal].HotGas += reheated_mass;
  Gal[centralgal].MetalsHotGas += reheated_mass * metallicity;
  Gal[centralgal].DustHotGas += reheated_mass * dust_fraction_cold;

  metallicityHot = get_metallicity(Gal[centralgal].HotGas, Gal[centralgal].MetalsHotGas);
  dust_fraction_hot = get_dust_fraction(Gal[p].HotGas, Gal[p].DustHotGas);
  dust_fraction_hotcold = get_dust_fraction(Gal[p].ColdGas + Gal[p].HotGas + Gal[p].EjectedMass, Gal[p].DustHotGas);
  // When dust_fraction_hotcold is calculated is irrelevant (i.e., before or after mass is added/subtracted from ColdGas).
  // This is because mass is conserved between ColdGas and HotGas during this step.

  Gal[centralgal].HotGas -= ejected_mass;
  Gal[centralgal].MetalsHotGas -= ejected_mass * metallicityHot;
  Gal[centralgal].DustHotGas -= ejected_mass * dust_fraction_hot; // Dust that has been ejected from the galaxy.
  Gal[centralgal].DustHotGas -= NSN * eta * Ms * dust_fraction_hotcold; // Dust that has been annihilated by shocks.
  dust_removed += NSN * eta * Ms * dust_fraction_hotcold;

  Gal[centralgal].EjectedMass += ejected_mass;
  Gal[centralgal].MetalsEjectedMass += ejected_mass * metallicityHot;
  Gal[centralgal].DustEjectedMass += ejected_mass * dust_fraction_hot;

  Gal[centralgal].EjectedMassSN += ejected_mass;

  Gal[p].GridOutflowRate[Halo[Gal[centralgal].HaloNr].SnapNum] += ejected_mass / dt;

  Gal[p].GrandSum -= mass_stars_recycled;

  if(Gal[p].ColdGas < 0.0) // Some final checks.
    Gal[p].ColdGas = 0.0;
  if(Gal[p].MetalsColdGas < 0.0)
    Gal[p].MetalsColdGas = 0.0;
  if(Gal[p].StellarMass < 0.0)
    Gal[p].StellarMass = 0.0;
  if(Gal[centralgal].HotGas < 0.0)
    Gal[centralgal].HotGas = 0.0;
  if(Gal[p].MetalsHotGas < 0.0)
    Gal[p].MetalsHotGas = 0.0;
  if(Gal[p].GrandSum < 0.0)
    Gal[p].GrandSum = 0.0;

  if (Gal[p].DustColdGas < 0.0)
    Gal[p].DustColdGas = 0.0;
  if (Gal[p].DustHotGas < 0.0)
    Gal[p].DustHotGas = 0.0;
}

void update_from_star_formation(int p, double stars, double dt, int step, bool ismerger, int tree, int ngal)
{

  double metallicity = get_metallicity(Gal[p].ColdGas, Gal[p].MetalsColdGas);
  double dust_fraction_cold = get_dust_fraction(Gal[p].ColdGas, Gal[p].DustColdGas);
  double time_spanned;

  // If the SF episode was from a merger the stars are placed in different reservoirs.
  if(!ismerger)
  {
    Gal[p].SfrDisk[step] += stars / dt;
    Gal[p].SfrDiskColdGas[step] = Gal[p].ColdGas;
    Gal[p].SfrDiskColdGasMetals[step] = Gal[p].MetalsColdGas;
  } else
  {

    Gal[p].SfrBulge[step] += stars / dt;
    Gal[p].SfrBulgeColdGas[step] += Gal[p].ColdGas;
    Gal[p].SfrBulgeColdGasMetals[step] += Gal[p].MetalsColdGas;

    Gal[p].BulgeMass += stars;
    Gal[p].MetalsBulgeMass += stars;

  }

  // If we're doing delayed SN we need to update the tracking
  if (ismerger == true)
  {
    time_spanned = 0.0;
  }
  else
  {
    time_spanned = dt;
  }

  //time_spanned = dt;

  if (IRA == 0)
  {
    update_SN_stars_array(p, stars, time_spanned, tree, ngal);
  }

  if (PhotonPrescription == 1)
  {
    update_stellar_tracking(p, stars, dt, tree, ngal);
  }

  // update gas and metals from star formation
  Gal[p].ColdGas -= stars;
  Gal[p].MetalsColdGas -= metallicity * stars;
  Gal[p].StellarMass += stars;
  Gal[p].MetalsStellarMass += metallicity * stars;
  Gal[p].DustColdGas -= dust_fraction_cold * stars; // When stars form they eat up dust.
  if (Gal[p].ColdGas < 0.0)
    Gal[p].ColdGas = 0.0;
  if (Gal[p].MetalsColdGas < 0.0)
    Gal[p].MetalsColdGas = 0.0;

}

void starformation_and_feedback(int p, int centralgal, double time, double dt, int halonr, int step, int tree, int ngal)
{
  double reff, tdyn, strdot, stars;
  double cold_crit;
  double reheated_mass = 0.0, mass_metals_new = 0.0, mass_stars_recycled = 0.0, ejected_mass = 0.0, NSN = 0.0;

  // Initialise variables
  strdot = 0.0;

  // First we check to see if we need to do delayed supernova feedback.

  if(IRA == 0)
  {
    do_previous_recycling(p, centralgal, step, dt);

    // The total amount that the reservoirs need to be adjusted by for this evolution step (i.e., ALL steps) has already been determined.
    // Within each substep (i.e., each `step` iteration) we ask 'how many iterations are left?' and then adjusts the reservoirs by this amount.
    // After the reservoirs have been updated, we subtract the amount of mass that we added from the running totals.
    //
    // These running totals are required because contemporaneous SN could cause extra mass to be added/removed within the next step iteration.

    update_from_SN_feedback(p, centralgal, Gal[p].reheated_mass / (STEPS-step), Gal[p].ejected_mass / (STEPS-step), Gal[p].mass_stars_recycled / (STEPS-step), Gal[p].mass_metals_new / (STEPS-step), Gal[p].NSN / (STEPS-step), dt);

    Gal[p].reheated_mass -= Gal[p].reheated_mass / (STEPS-step);
    Gal[p].ejected_mass -= Gal[p].ejected_mass / (STEPS-step);
    Gal[p].mass_stars_recycled -= Gal[p].mass_stars_recycled / (STEPS-step);
    Gal[p].mass_metals_new -= Gal[p].mass_metals_new / (STEPS-step);
    Gal[p].NSN -= Gal[p].NSN / (STEPS-step);

  }

  // star formation recipes
  if(SFprescription == 0)
  {
    // we take the typical star forming region as 3.0*r_s using the Milky Way as a guide
    reff = 3.0 * Gal[p].DiskScaleRadius;
    tdyn = reff / Gal[p].Vvir;

    // from Kauffmann (1996) eq7 x piR^2, (Vvir in km/s, reff in Mpc/h) in units of 10^10Msun/h
    cold_crit = 0.19 * Gal[p].Vvir * reff;

    if(Gal[p].ColdGas > cold_crit && tdyn > 0.0)
    	strdot = SfrEfficiency * (Gal[p].ColdGas - cold_crit) / tdyn;
    else
    {
    	strdot = 0.0;
    }
  }
  else
  {
    printf("No star formation prescription selected!\n");
    ABORT(0);
  }

  stars = strdot * dt;
  if(stars < 0.0)
    stars = 0.0;

  if(SupernovaRecipeOn == 1)
  {
    if(IRA == 0)
      do_contemporaneous_SN(p, centralgal, dt, &stars, &reheated_mass, &mass_metals_new, &mass_stars_recycled, &ejected_mass, &NSN);
    else if(IRA == 1)
    {
      do_IRA_SN(p, centralgal, &stars, &reheated_mass, &mass_metals_new, &mass_stars_recycled, &ejected_mass, &NSN);
    }

  }

  if(stars > Gal[p].ColdGas) // we do this check in 'do_current_sn()' but if supernovarecipeon == 0 then we don't do the check.
  {
    double factor = Gal[p].ColdGas / stars;
    stars *= factor;
  }

  update_from_star_formation(p, stars, dt, step, false, tree, ngal);

  update_from_SN_feedback(p, centralgal, reheated_mass, ejected_mass, mass_stars_recycled, mass_metals_new, NSN, dt);

  // check for disk instability
  if(DiskInstabilityOn)
    check_disk_instability(p, centralgal, halonr, time, dt, step, tree, ngal);

}

// This function will calculate the fraction of stars formed in snapshot i that go nova in snapshot j.
// Will also calculate the mass fraction of stars formed in snapshot i that go nova in snapshot j.
// Based off equation (17) and (24) of Mutch et al. (2016).
//
// INPUT: The bounds for the integral (m_low and m_high). ## UNITS: Msun.
// 	: Pointers to modify the values for number and mass fraction of stars that go nova (*Delta_Eta and *Delta_m). ## UNITS: Msun^-1 (Delta_Eta) and unitless (Delta_m).
//
// OUTPUT: None but delta_Eta, the number fraction of stars, and delta_m, the mass fraction of stars, will be modified.

void calculate_Delta_Eta(double m_low, double m_high, double *Delta_Eta, double *Delta_m)
{

  if (m_low < 8.0)
    m_low = 8.0;

  if (m_high < 8.0)
    m_high = 8.0;

  int32_t bin_idx_low = round((m_low - m_IMFbins_low) / (m_IMFbins_delta));
  int32_t bin_idx_high = round((m_high - m_IMFbins_low) / (m_IMFbins_delta));

  if (bin_idx_high == N_massbins)
    --bin_idx_high;

  if (bin_idx_low == N_massbins)
    --bin_idx_low;

  *Delta_Eta = IMF_CONSTANT * (IMF_massgrid_eta[bin_idx_high] - IMF_massgrid_eta[bin_idx_low]);
  *Delta_m = IMF_CONSTANT * (IMF_massgrid_m[bin_idx_high] - IMF_massgrid_m[bin_idx_low]);

}

// This function determines the amount of mass reheated for a supernova event.
//
// INPUT: The number fraction of stars that have gone nova (Delta_Eta) ## UNITS: Msun^-1.
// 	: The mass of stars formed (stars). ## UNITS: 1.0e10Msun/h (code units).
// 	: The maximum velocity of the galaxy disc (Vmax). ## UNITS: km/s.
//
// OUTPUT: The amount of mass reheated by the supernova event. ## UNITS: 1.0e10Msun/h (code units).

double calculate_reheated_mass(double Delta_Eta, double stars, double Vmax)
{

 // if(stars < 1e-10)
 //   return 0.0;

  double epsilon_mass;

  if(RescaleSN == 0)
  {
    epsilon_mass = FeedbackReheatingEpsilon;
  } else
  {
    epsilon_mass = alpha_mass * (0.5 + pow(Vmax/V_mass, -beta_mass));
  }

  if (epsilon_mass > epsilon_mass_max)
  {
    epsilon_mass = epsilon_mass_max; // We enforce a maximum value for the mass loading factor.
  }

  double reheated_mass = Delta_Eta/Eta_SNII * stars * epsilon_mass;

  return reheated_mass;

}

// This function determines the amount of energy injected from a supernova event.
//
// INPUT: The number fraction of stars that have gone nova (Delta_Eta) ## UNITS: Msun^-1.
// 	: The mass of stars formed (stars). ## UNITS: 1.0e10Msun/h (code units).
// 	: The maximum velocity of the galaxy disc (Vmax). ## UNITS: km/s.
//
// OUTPUT: The amount of energy from a supernova event. ## UNITS: erg.

double calculate_reheated_energy(double Delta_Eta, double stars, double Vmax)
{

  double epsilon_energy;

  if(RescaleSN == 0)
  {
    epsilon_energy = FeedbackEjectionEfficiency;
  } else
  {
    epsilon_energy = alpha_energy * (0.5 + pow(Vmax/V_energy, -beta_energy));
  }

  double reheated_energy = 0.5 * Delta_Eta * stars * epsilon_energy * EnergySN; // Delta_Eta was in Msun so need to convert to stars to Msun.

  return reheated_energy;

}

// This function calculates the mass of stars that would have gone supernova in time t.
// Function and fit taken from Portinari et al. (1998).
//
// INPUT: Time (in Myr).
//
// OUTPUT: Mass of stars (in Msun) that would have gone Nova in time t.

double calculate_coreburning(double t)
{

  /*
  double a = 0.7473; // Fits from Portinari et al. (1998).
  double b = -2.6979;
  double c = -4.7659;
  double d = 0.5934;

  double m = pow(10, a/log10(t) + b * exp(c/log10(t)) + d);

  return m;
  */

  int32_t bin_idx = (t - coreburning_tbins_low) / (coreburning_tbins_delta);

  if (bin_idx < 0) // Time is so short that only stars with mass greater than 120Msun can go nova.
    return 120.0;

  if (bin_idx > N_tbins) // Time is so long that all stars within the IMF range can go nova.
    return 8.0;

  return coreburning_times[bin_idx];

}

// If the cold gas that has been reheated has enough energy, it is possible to unbind some of the gas in the hot halo.
// This function determines this threshhold and how much mass is unbound and ejected.
//
// INPUT: Amount of cold gas reheated by supernova events (reheated_mass). ## UNITS: 1.0e10Msun/h (code units).
// 	: Amount of energy injected from the supernova events (reheated_energy). ## UNITS: erg.
// 	: Virial velocity of the background FoF halo (Vvir). ## UNITS: km/s.
//
// OUTPUT: Mass of ejected material. ## UNITS: 1.0e10Msun/h (code units).
//
// Note: Be aware here that a number of unit changes have occurred.
// This is necessary because of how Joules are defined.

double calculate_ejected_mass(double *reheated_mass, double reheated_energy, double Vvir)
{
  double ejected_mass;
  double Delta_Ehot;
  const double Joule_to_erg = 0.5 * SOLAR_MASS / 1.0e3 * 1.0e3 * 1.0e3 * 1.0e7; // Note: This is ergs assuming the mass used is in kg.
                                                                                // As we use 1.0e10 Msun/h our energy isn't quite ergs.
                                                                                // However we have kept this consistent throughout allowing the factor to be cancelled.
  const double inv_Joule_to_erg = 1.0 / Joule_to_erg;

  if(Vvir > 0.0)
  {
    Delta_Ehot = *reheated_mass * Vvir * Vvir * Joule_to_erg; // Change in the thermal energy of the hot resevoir in erg..
    // Note the units here. Want final answer in Ergs (which is 1e-7 Joules) which has SI units of kg m^2 s^-2.
    // We change the mass from 1.0e10Msun/h to kg.  Change virial velocity from km/s to m/s. Finally change Joules to erg.

    if(reheated_energy > Delta_Ehot) // Have enough energy to have some ejected mass.
    {
    	ejected_mass = (reheated_energy - Delta_Ehot) * inv_Joule_to_erg / (Vvir * Vvir); // Balance between the excess thermal energy and the thermal energy of the hot gas.
      // Energy changed from erg to Joules (kg m^2 s^-2) then divided by virial velocity (converted to m/s) giving final units of kg.

    }
    else // Not enough energy to eject mass.
    {
      ejected_mass = 0.0;
      *reheated_mass = reheated_energy * inv_Joule_to_erg / (Vvir * Vvir); // Amount of mass that can be reheated to virial temperature of the hot halo.
    }

  }
  else
    ejected_mass = 0.0;

  if(ejected_mass < 0.0)
    ejected_mass = 0.0;

  return ejected_mass;

}

// This function answers the question, "How many stars formed in the previous time steps will explode during the current supernova timestep."
// Note: Since the time scale on which we calculate SN feedback is variable (defined by the user in 'TimeResolutionSN') we use the terminology 'supernova timestep'.
// Double note: This function is different to "do_contemporaneous_SN" because that looks at the stars formed during the CURRENT STAR FORMATION time step, not previous ones.
void do_previous_SN(int p, int centralgal, double dt)
{
  double reheated_mass = 0.0;
  double reheated_energy = 0.0;
  double mass_stars_recycled = 0.0;
  double mass_metals_new = 0.0;
  double ejected_mass = 0.0;
  double NSN = 0.0; // This is the number of supernova events during this time step.

  int i;
  double m_low = -1.0, m_high = -1.0, t_low = -1.0, t_high = -1.0;
  double Delta_Eta, Delta_m;
  double time_until_next_SN;

  if(dt * UnitTime_in_Megayears / Hubble_h < TimeResolutionSN) // If the star formation time scale is smaller than the time scale on which we do SN feedback
    time_until_next_SN = TimeResolutionSN; // Then the time that we next calculate delayed SN feedback would be given dictated by the time resolution of SN.
  else // Otherwise the star formation time scale is larger than the SN feedback time scale.
  {
    time_until_next_SN = dt * UnitTime_in_Megayears / Hubble_h;  // Then the time that we next calculate delayed SN feedback is in the next SF timestep.
  }

  if(SupernovaRecipeOn == 1) // Calculating the ejected mass due to previous star formation episodes.
  {
    for(i = 1; i < SN_Array_Len; ++i)
    {
      if(Gal[p].SN_Stars[i] < 1e-10)
	    continue;

      // First calculate the smallest star which would have expended its H and He and gone supernova.
      // This defines the mass boundary below which a star will not go nova in the current SN step.

      t_low = (i * TimeResolutionSN) / 2.0 + time_until_next_SN; // (i * TimeResolutionSN) is the number of stars that were formed that many megayears ago.  time_until_next_SN allows us to ask how many of stars will explode in this SN step.
      if(t_low < 2)
  	    t_low = 2;
      if (t_low > 45)
        m_low = 120.0;

      m_low = calculate_coreburning(t_low);

      if (m_low < 8.0) // We enforce that SN-II only occur in stars M > 8Msun.
        m_low = 8.0;

      if (m_low > 120.0)
        continue;

      // Next we calculate the largest stellar mass which would have expended its H and He and gone supernova.
      // This defines the mass boundary beyond which a star would have gone nova BEFORE the current SN step.

      t_high = (i * TimeResolutionSN) / 2.0; // (i * TimeResolutionSN) is the number of stars that were formed that many Megayears ago.  The / 2.0 factor arises from the fact that we assume stars were formed in a single co-eval burst in the middle of this period.
      if(t_high < 2)
        t_high = 2;

      m_high = calculate_coreburning(t_high);

      if (m_high < 8) // In this instance every star that has formed in i Myr ago has already gone supernova and accounted for.
          continue;

      if (m_high > 120.0)
          continue;

      calculate_Delta_Eta(m_low, m_high, &Delta_Eta, &Delta_m); // Calculate the number and mass fraction of stars i Myr ago that go supernova in the current SF step.
      reheated_mass += calculate_reheated_mass(Delta_Eta, Gal[p].SN_Stars[i], Gal[centralgal].Vmax); // Update the amount of mass reheated from previous stars that have gone nova.

      reheated_energy += calculate_reheated_energy(Delta_Eta, Gal[p].SN_Stars[i], Gal[centralgal].Vmax); // Update the energy injected from previous stars that have gone nova.
      mass_stars_recycled += Delta_m * Gal[p].SN_Stars[i]; // Update the amount of stellar mass recycled from previous stars that have gone nova.

      // Since stars have gone supernova, need to update the number of stars remaining at that time.
      if (PhotonPrescription == 1)
      {
        Gal[p].Stellar_Stars[i] -= Delta_m * Gal[p].SN_Stars[i];
        if (Gal[p].Stellar_Stars[i] < 0.0)
          Gal[p].Stellar_Stars[i] = 0.0;
      }

      Gal[p].SN_Stars[i] -= Delta_m * Gal[p].SN_Stars[i];
      if (Gal[p].SN_Stars[i] < 0.0)
        Gal[p].SN_Stars[i] = 0.0;

      mass_metals_new += Delta_m / m_SNII * Yield * Gal[p].SN_Stars[i]; // Update the amount of new metals that the supernova has enriched the ISM with.

      NSN += Delta_Eta * Gal[p].SN_Stars[i] * 1.0e10 / Hubble_h; // The number of supernova events will be simply given by the number fraction of stars that exploded times the mass of stars.

      XASSERT(reheated_mass >= 0.0, "i = %d \t Reheated mass = %.4e \t t_low = %.4e Myr \t m_low = %.4e \t t_high = %.4e Myr \t m_high = %.4e\n", i, reheated_mass, t_low, m_low, t_high, m_high); // Just make sure we're doing this right.

    }

    if(reheated_mass > Gal[p].ColdGas) // Can't reheated more cold gas than we currently have.
        reheated_mass = Gal[p].ColdGas;

    XASSERT(reheated_mass >= 0.0, "Reheated mass = %.4e \t t_low = %.4e Myr \t m_low = %.4e \t t_high = %.4e Myr \t m_high = %.4e\n", reheated_mass, t_low, m_low, t_high, m_high); // Just make sure we're doing this right.
    assert(reheated_energy >= 0.0);
    assert(mass_stars_recycled >= 0.0);
    assert(mass_metals_new >= 0.0);

    ejected_mass = calculate_ejected_mass(&reheated_mass, reheated_energy, Gal[centralgal].Vmax); // Calculate the amount of mass ejected from supernova events.

    XASSERT(ejected_mass >= 0.0, "For galaxy %d the ejected mass was %.4e \t reheated_mass = %.4e\n", p, ejected_mass, reheated_mass);
  }

  Gal[p].reheated_mass = reheated_mass;
  Gal[p].ejected_mass = ejected_mass;
  Gal[p].mass_stars_recycled = mass_stars_recycled;
  Gal[p].mass_metals_new = mass_metals_new;
  Gal[p].NSN = NSN;

}

/*
void do_current_SN(int p, int centralgal, int halonr, double *stars, double *reheated_mass, double *mass_metals_new, double *mass_stars_recycled, double *ejected_mass)
{

  *reheated_mass = 0.0;
  double reheated_energy = 0.0;
  *mass_stars_recycled = 0.0;
  *mass_metals_new = 0.0;
  *ejected_mass = 0.0;

  double Delta_Eta = Eta_SNII;
  double Delta_m = RecycleFraction;
  double mass_fraction = 1.0;

  *mass_stars_recycled = Delta_m * (*stars);
  *reheated_mass = calculate_reheated_mass(Delta_Eta, *stars, Gal[centralgal].Vmax); // Update the amount of mass reheated from stars that have gone nova.

  if(*reheated_mass > Gal[p].ColdGas) // Can only reheat as much cold gas we have available.
      *reheated_mass = Gal[p].ColdGas;
  if((*reheated_mass + *stars > Gal[p].ColdGas) && (*reheated_mass + *stars > 0.0))
  {
    double factor = Gal[p].ColdGas / (*reheated_mass + *stars);
    *reheated_mass *= factor;
    *stars *= factor;
  }

  *mass_metals_new = mass_fraction * Yield * (*stars); // Update the amount of new metals that the supernova has enriched the ISM with.
  reheated_energy = calculate_reheated_energy(Delta_Eta, *stars, Gal[centralgal].Vmax); // Update the energy injected from previous stars that have gone nova.

  assert(*reheated_mass >= 0.0); // Just make sure we're doing this right.
  assert(reheated_energy >= 0.0);
  assert(*mass_stars_recycled >= 0.0);
  assert(*mass_metals_new >= 0.0);

  *ejected_mass = calculate_ejected_mass(&(*reheated_mass), reheated_energy, Gal[centralgal].Vvir); // Calculate the amount of mass ejected from supernova events.

}
*/

void do_previous_recycling(int p, int centralgal, int step, double dt)
{

  if(Gal[p].StellarAge_Numerator > 0.0)
  {
    double t_low, t_high;
    double m_low, m_high;
    double Delta_Eta, Delta_m;
    double mass_stars_recycled;
    double NSN;

    double mwmsa = Gal[p].StellarAge_Numerator / Gal[p].StellarAge_Denominator; // The numerator is weighted by the age of the stars with the denominator simply being the number of stars.
    double time_into_snap = dt * step;

    t_high = (mwmsa - ((Age[Gal[p].SnapNum - 1]) - time_into_snap)) * UnitTime_in_Megayears / Hubble_h;
    if (t_high < 2.0)
      t_high = 2.0;
    if (t_high > 45.0)
      m_high = 8.0;
    else
      m_high = calculate_coreburning(t_high);

    t_low = (mwmsa - (Age[Gal[p].SnapNum] - time_into_snap)) * UnitTime_in_Megayears / Hubble_h;
    if (t_low < 2.0)
      t_low = 2.0;

    if (t_low > 45.0)
      m_low = 8.0;
    else
      m_low = calculate_coreburning(t_low);

    calculate_Delta_Eta(m_low, m_high, &Delta_Eta, &Delta_m); // Calculate the number and mass fraction of stars from snapshot i that go nova in the snapshot we are evolving FROM.
    mass_stars_recycled = Delta_m * Gal[p].StellarAge_Denominator; // Update the amount of stellar mass recycled from previous stars that have gone nova.
    NSN = Delta_Eta * Gal[p].StellarAge_Denominator * 1.0e10 / Hubble_h;
    update_from_SN_feedback(p, centralgal, 0.0, 0.0, mass_stars_recycled, 0.0, NSN, dt);
  }
}

// In this function we answer the question, "How many stars formed in the current star formation time step will explode by the time we next calculate our supernova feedback?"
// Note: This function only concerns itself with stars formed in the CURRENT STAR FORMATION time step, unlike "do_previous_SN" which focuses on calculated SN feedback from PREVIOUS STAR FORMATION time steps.

void do_contemporaneous_SN(int p, int centralgal, double dt, double *stars, double *reheated_mass, double *mass_metals_new, double *mass_stars_recycled, double *ejected_mass, double *NSN)
{

  *reheated_mass = 0.0;
  double reheated_energy = 0.0;
  *mass_stars_recycled = 0.0;
  *mass_metals_new = 0.0;
  *ejected_mass = 0.0;

  double m_low, m_high, t_low;
  double Delta_Eta, Delta_m;
  double fac;

  if(dt * UnitTime_in_Megayears / Hubble_h / 2.0 < TimeResolutionSN) // If the star formation time scale is smaller than the time scale on which we do SN feedback
    t_low = (TimeResolutionSN - Gal[p].Total_SN_SF_Time); // Then our 'sub-grid' SN feedback time will be the time from this SF episode until we next calculate SN feedback (divided by 2 as we assume the stars are formed in the middle of the interval).
  else // Otherwise the star formation time scale is larger than the SN feedback time scale.
    t_low  = (dt * UnitTime_in_Megayears / Hubble_h) / 2.0; // Then the feedback time will be the time from this SF event to the next star formation event. This is because SN feedback is only calculated when star formation occurs regardless of the actual value of 'TimeResolutionSN'.

  if (t_low < 2.0)
    t_low = 2.0;

  m_low = calculate_coreburning(t_low);
  if(m_low < 8.0)
    m_low = 8.0;

  if(m_low > 120.0)
    return;

  m_high = 120.0;

  calculate_Delta_Eta(m_low, m_high, &Delta_Eta, &Delta_m); // Calculate the number and mass fraction of stars i Myr ago that go supernova in the current SF step.
  *reheated_mass = calculate_reheated_mass(Delta_Eta, *stars, Gal[centralgal].Vmax); // Update the amount of mass reheated from previous stars that have gone nova.

  if((*stars + *reheated_mass) > Gal[p].ColdGas && (*stars + *reheated_mass) > 0.0)
  {
    fac = Gal[p].ColdGas / (*stars + *reheated_mass);
    *stars *= fac;
    *reheated_mass *= fac;
  }

  reheated_energy += calculate_reheated_energy(Delta_Eta, *stars, Gal[centralgal].Vmax); // Update the energy injected from previous stars that have gone nova.
  *mass_stars_recycled += Delta_m * (*stars); // Update the amount of stellar mass recycled from previous stars that have gone nova.
  *mass_metals_new += Delta_m / m_SNII * Yield * (*stars); // Update the amount of new metals that the supernova has enriched the ISM with.

  if(*reheated_mass > Gal[p].ColdGas) // Can't reheated more cold gas than we currently have.
    *reheated_mass = Gal[p].ColdGas;

  *ejected_mass = calculate_ejected_mass(&(*reheated_mass), reheated_energy, Gal[centralgal].Vmax); // Calculate the amount of mass ejected from supernova events.
  *NSN = Delta_Eta * (*stars * 1.0e10 / Hubble_h);

  assert(*reheated_mass >= 0.0); // Just make sure we're doing this right.
  assert(reheated_energy >= 0.0);
  assert(*mass_stars_recycled >= 0.0);
  assert(*mass_metals_new >= 0.0);
  assert(*ejected_mass >= 0.0);

}

void do_IRA_SN(int p, int centralgal, double *stars, double *reheated_mass, double *mass_metals_new, double *mass_stars_recycled, double *ejected_mass, double *NSN)
{

  double fac;
  *reheated_mass = 0.0;
  *mass_metals_new = 0.0;
  *mass_stars_recycled = 0.0;
  *ejected_mass = 0.0;

  *reheated_mass = calculate_reheated_mass(Eta_SNII, *stars, Gal[centralgal].Vmax);
  *NSN = Eta_SNII * (*stars);

  if((*stars + *reheated_mass) > Gal[p].ColdGas && (*stars + *reheated_mass) > 0.0)
  {
    fac = Gal[p].ColdGas / (*stars + *reheated_mass);
    *stars *= fac;
    *reheated_mass *= fac;
  }

  double reheated_energy = calculate_reheated_energy(Eta_SNII, *stars, Gal[centralgal].Vmax);

  *ejected_mass = calculate_ejected_mass(&(*reheated_mass), reheated_energy, Gal[centralgal].Vmax);

  if(*ejected_mass < 0.0)
    *ejected_mass = 0.0;

  *mass_metals_new = Yield * (*stars);
  *mass_stars_recycled = RecycleFraction * (*stars);

}

// Local Functions //

// This function updates the array which holds the amount of stars formed in the past 50 Myr.
// As we only need this array for doing delayed SN, if we are using the IRA then this function will never be called.
// If the SF timestep is less than the resolution on which we do supernova feedback then we will keep track of the total stars formed over these small timesteps and then update them all in one bin.
// If the SF timestep is larger than the resolution on which we do supernova feedback we will spread the stars formed evenly over a number of elements (>= 1).
//
// INPUT: The index of the galaxy (p).
// 	: The number of stars formed in the current SF episode (stars). ## UNITS: 1.0e10 Msun/h (Code Units).
// 	: The timestep over which the SF is occuring (dt). ## UNITS: Code Units, multiply by 'UnitTime_in_Megayears / Hubble_h' for Myr.
// 	: The tree currently being used (tree); currently used for debugging purposes.
//
// OUTPUT: None.

void update_SN_stars_array(int p, double stars, double dt, int tree, int ngal)
{

  double time_spanned = dt * UnitTime_in_Megayears / Hubble_h; // The time spanned by this star formation event.

  Gal[p].Total_SN_SF_Time += time_spanned; // How long it has been since we've updated the array?
  Gal[p].Total_SN_Stars += stars; // How many stars we will need to bin once we do update the array?

  if(Gal[p].Total_SN_SF_Time > TimeResolutionSN * SN_Array_Len) // This handles cases in which the time spanned is greater than 50Myr.  In this case we wipe the array clean and push the star formation into a 50Myr bin.
    Gal[p].Total_SN_SF_Time = TimeResolutionSN * SN_Array_Len;

  if(Gal[p].Total_SN_SF_Time < TimeResolutionSN) // If it hasn't been long enough yet, don't update the array.
    return;

  int num_shuffled = round(Gal[p].Total_SN_SF_Time/TimeResolutionSN); // How many cells will this SF event span.

  double stars_spread = Gal[p].Total_SN_Stars/num_shuffled; // We spread the stars evenly over time.

  int i;

  XASSERT(Gal[p].IsMalloced == 1, "We are attempting to update the stars array but this galaxy has already had its arrays freed.\nGalaxy %d \t Halo %d \t Tree %d \t Time spanned %.4eMyr\n", p, Gal[p].HaloNr, tree, time_spanned);

  for(i = SN_Array_Len - 1; i > num_shuffled - 1; --i)
  {
      Gal[p].SN_Stars[i] = Gal[p].SN_Stars[i-num_shuffled]; // Shuffle the current elements of the array far enough along so we can store the new stars.
  }
    XPRINT(p < ngal, "We have the case where the galaxy p is greater than the number of galaxies.  p = %d \t ngal = %d\n", p, ngal);
  for(i = SN_Array_Len - 1; i > (SN_Array_Len - num_shuffled - 1); --i)
  {
    Gal[p].StellarAge_Numerator += Gal[p].SN_Stars[i] * (Age[Gal[p].SnapNum - 4]);
    Gal[p].StellarAge_Denominator += Gal[p].SN_Stars[i];
  }

  for(i = 0; i < num_shuffled; ++i)
  {
    Gal[p].StellarAge_Numerator += Gal[p].SN_Stars[i] * (Age[Gal[p].SnapNum - 4]);
    Gal[p].SN_Stars[i] = stars_spread; // Update the vacated elements with the new stars.
    Gal[p].GrandSum += stars_spread;
  }

  Gal[p].Total_SN_SF_Time = 0.0; // We've updated so reset our variables.
  Gal[p].Total_SN_Stars = 0.0;
}

void update_stellar_tracking(int p, double stars, double dt, int tree, int ngal)
{

  double time_spanned = dt * UnitTime_in_Megayears / Hubble_h; // The time spanned by this star formation event.

  Gal[p].Total_Stellar_SF_Time += time_spanned; // How long it has been since we've updated the array?
  Gal[p].Total_Stellar_Stars += stars; // How many stars we will need to bin once we do update the array?

  if(Gal[p].Total_Stellar_SF_Time > TimeResolutionStellar * StellarTracking_Len)
    Gal[p].Total_Stellar_SF_Time = TimeResolutionStellar * StellarTracking_Len;

  if(Gal[p].Total_Stellar_SF_Time < TimeResolutionStellar) // If it hasn't been long enough yet, don't update the array.
    return;

  int num_shuffled = round(Gal[p].Total_Stellar_SF_Time/TimeResolutionStellar); // How many cells will this SF event span.
  double stars_spread = Gal[p].Total_Stellar_Stars/num_shuffled; // We spread the stars evenly over time.

  int i;

  XASSERT(Gal[p].IsMalloced == 1, "We are attempting to update the Stellar Tracking array but this galaxy has already had its arrays freed.\nGalaxy %d \t Halo %d \t Tree %d \t Time spanned %.4eMyr\n", p, Gal[p].HaloNr, tree, time_spanned);

  for(i = StellarTracking_Len - 1; i > num_shuffled - 1; --i)
  {
      Gal[p].Stellar_Stars[i] = Gal[p].Stellar_Stars[i-num_shuffled]; // Shuffle the current elements of the array far enough along so we can store the new stars.
  }
  XPRINT(p < ngal, "We have the case where the galaxy p is greater than the number of galaxies.  p = %d \t ngal = %d\n", p, ngal);

  for(i = 0; i < num_shuffled; ++i)
  {
    Gal[p].Stellar_Stars[i] = stars_spread; // Update the vacated elements with the new stars.
  }

  Gal[p].Total_Stellar_SF_Time = 0.0; // We've updated so reset our variables.
  Gal[p].Total_Stellar_Stars = 0.0;
}

max_stars_repo_name: jacobseiler/rsage

max_issues_repo_name: jacobseiler/rsage

max_forks_repo_name: jacobseiler/rsage


text:
#pragma once

#include "JsonHandler.h"
#include "Library.h"

#include <gsl/span>

#include <cstddef>
#include <optional>
#include <string>
#include <vector>

namespace rapidjson {
struct MemoryStream;
}

namespace CesiumJsonReader {

/**
 * @brief The result of {@link Reader::readJson}.
 */
template <typename T> struct ReadJsonResult {
  /**
   * @brief The value read from the JSON, or `std::nullopt` on error.
   */
  std::optional<T> value;

  /**
   * @brief Errors that occurred while reading.
   */
  std::vector<std::string> errors;

  /**
   * @brief Warnings that occurred while reading.
   */
  std::vector<std::string> warnings;
};

/**
 * @brief Reads JSON.
 */
class CESIUMJSONREADER_API JsonReader {
public:
  /**
   * @brief Reads JSON from a byte buffer.
   *
   * @param data The buffer from which to read JSON.
   * @param handler The handler to receive the top-level JSON object.
   * @return The result of reading the JSON.
   */
  template <typename T>
  static ReadJsonResult<typename T::ValueType>
  readJson(const gsl::span<const std::byte>& data, T& handler) {
    ReadJsonResult<typename T::ValueType> result;

    result.value.emplace();

    FinalJsonHandler finalHandler(result.warnings);
    handler.reset(&finalHandler, &result.value.value());

    JsonReader::internalRead(
        data,
        handler,
        finalHandler,
        result.errors,
        result.warnings);

    if (!result.errors.empty()) {
      result.value.reset();
    }

    return result;
  }

private:
  class FinalJsonHandler : public JsonHandler {
  public:
    FinalJsonHandler(std::vector<std::string>& warnings);
    virtual void reportWarning(
        const std::string& warning,
        std::vector<std::string>&& context) override;
    void setInputStream(rapidjson::MemoryStream* pInputStream) noexcept;

  private:
    std::vector<std::string>& _warnings;
    rapidjson::MemoryStream* _pInputStream;
  };

  static void internalRead(
      const gsl::span<const std::byte>& data,
      IJsonHandler& handler,
      FinalJsonHandler& finalHandler,
      std::vector<std::string>& errors,
      std::vector<std::string>& warnings);
};

} // namespace CesiumJsonReader

max_stars_repo_name: yieryi/cesium-native

max_issues_repo_name: yieryi/cesium-native

max_forks_repo_name: yieryi/cesium-native