text:
/*
 Copyright (C) 2004-2007 M. Oliveira, F. Nogueira

 This program is free software; you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
 the Free Software Foundation; either version 2, or (at your option)
 any later version.

 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 GNU General Public License for more details.

 You should have received a copy of the GNU General Public License
 along with this program; if not, write to the Free Software
 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
 02111-1307, USA.
*/

#include <config.h>
#include <string_f.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_odeiv.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_multiroots.h>

/* This file contains interface routines that allow Fortran 90/95
routines to call GSL routines. The names of these routines are
compiler dependent: Fortran 90/95 compilers have different
name-mangling schemes.  For a description of the routines please refer
to the GSL documentation. */

          /* Error Handling */

void FC_FUNC_(gsl_strerror, GSL_STRERROR)
     (int *err, STR_F_TYPE res STR_ARG1)
{
  char *c;

  c = gsl_strerror(*err);
  /*TO_F_STR1(c, res);!!not too sure what this does*/
}

          /* Mathematical Functions */

double FC_FUNC_(gsl_asinh, GSL_ASINH)
     (const double *x)
{
  return gsl_asinh(*x);
}

          /* Special Functions */

double FC_FUNC_(gsl_sf_bessel_knu_scaled, GSL_SF_BESSEL_KNU_SCALED)
     (double *nu, double *x)
{
  return gsl_sf_bessel_Knu_scaled(*nu,*x);
}

          /* Vectors and Matrices */

void FC_FUNC_(gsl_vector_alloc, GSL_VECTOR_ALLOC)
     (int *n, void **v)
{
  *v = (void*) gsl_vector_alloc (*n);
}

void FC_FUNC_(gsl_vector_free, GSL_VECTOR_FREE)
     (void **v)
{
  gsl_vector_free ((gsl_vector *)(*v));
}

double FC_FUNC_(gsl_vector_get, GSL_VECTOR_GET) (void **v, int *i)
{
  return gsl_vector_get ((gsl_vector *)(*v), *i);
}

void FC_FUNC_(gsl_vector_set, GSL_VECTOR_SET)
     (void **v, int *i, double *x)
{
  gsl_vector_set ((gsl_vector *)(*v), *i, *x);
}

void FC_FUNC_(gsl_matrix_alloc, GSL_MATRIX_ALLOC)
     (int *n1, int *n2, void **m)
{
  *m = (void*) gsl_matrix_alloc (*n1, *n2);
}

void FC_FUNC_(gsl_matrix_calloc, GSL_MATRIX_CALLOC)
     (int *n1, int *n2, void **m)
{
  *m = (void*) gsl_matrix_calloc (*n1, *n2);
}

void FC_FUNC_(gsl_matrix_free, GSL_MATRIX_FREE)
     (void **m)
{
  gsl_matrix_free ((gsl_matrix *)(*m));
}

double FC_FUNC_(gsl_matrix_get, GSL_MATRIX_GET)
     (void **m, int *i, int *j)
{
  return gsl_matrix_get ((gsl_matrix *)(*m), *i, *j);
}

void FC_FUNC_(gsl_matrix_set, GSL_MATRIX_SET)
     (void **m, int *i, int *j, double *x)
{
  gsl_matrix_set ((gsl_matrix *)(*m), *i, *j, *x);
}

          /* Permutations */

void FC_FUNC_(gsl_permutation_alloc, GSL_PERMUTATION_ALLOC)
     (int *n, void **p)
{
  *p = (void*) gsl_permutation_alloc (*n);
}

void FC_FUNC_(gsl_permutation_free, GSL_PERMUTATION_FREE)
     (void **p)
{
  gsl_permutation_free ((gsl_permutation *)(*p));
}

          /* Linear Algebra */

int FC_FUNC_(gsl_linalg_lu_decomp, GSL_LINALG_LU_DECOMP)
     (void **m, void **p, int *signum)
{
  int s,ierr;
  ierr = gsl_linalg_LU_decomp ((gsl_matrix *)(*m) , (gsl_permutation *)(*p), &s);
  *signum = s;
  return (ierr);
}

int FC_FUNC_(gsl_linalg_lu_invert, GSL_LINALG_LU_INVERT)
     (void **LU, void **p, void **inverse)
{

  return gsl_linalg_LU_invert ((gsl_matrix *)(*LU),(gsl_permutation *)(*p), (gsl_matrix *)(*inverse));
}

int FC_FUNC_(gsl_linalg_lu_solve, GSL_LINALG_LU_SOLVE)
     (void **LU, void **p, void **b, void **x)
{
  return gsl_linalg_LU_solve ((gsl_matrix *)(*LU), (gsl_permutation *)(*p), (gsl_vector *)(*b), (gsl_vector *)(*x));
}

          /* Ordinary Differential Equations */

void FC_FUNC_(gsl_odeiv_step_alloc, GSL_ODEIV_STEP_ALLOC)
     (const int *stepping_func, const int *dim, void **stp)
{
  /* WARNING: This routine has a small modification: the user needs to
   pass an integer in order to choose the stepping function instead of
   a variable of type gsl_odeiv_step_type. */

  switch(*stepping_func) {
  case 1 : // Embedded 2nd order Runge-Kutta with 3rd order error estimate
    *stp = (void *)gsl_odeiv_step_alloc(gsl_odeiv_step_rk2, *dim);
    break;
  case 2 : // 4th order (classical) Runge-Kutta
    *stp = (void *)gsl_odeiv_step_alloc(gsl_odeiv_step_rk4, *dim);
    break;
  case 3 : // Embedded 4th order Runge-Kutta-Fehlberg method with 5th order error estimate. This method is a good general-purpose integrator
    *stp = (void *)gsl_odeiv_step_alloc(gsl_odeiv_step_rkf45, *dim);
    break;
  case 4 : // Embedded 4th order Runge-Kutta Cash-Karp method with 5th order error estimate
    *stp = (void *)gsl_odeiv_step_alloc(gsl_odeiv_step_rkck, *dim);
    break;
  case 5 : // Embedded 8th order Runge-Kutta Prince-Dormand method with 9th order error estimate
    *stp = (void *)gsl_odeiv_step_alloc(gsl_odeiv_step_rk8pd, *dim);
    break;
  case 6 : // Implicit 2nd order Runge-Kutta at Gaussian points
    *stp = (void *)gsl_odeiv_step_alloc(gsl_odeiv_step_rk2imp, *dim);
    break;
  case 7 : // Implicit 4th order Runge-Kutta at Gaussian points
    *stp = (void *)gsl_odeiv_step_alloc(gsl_odeiv_step_rk4imp, *dim);
    break;
  case 8 : // Implicit Bulirsch-Stoer method of Bader and Deuflhard. This algorithm requires the Jacobian
    *stp = (void *)gsl_odeiv_step_alloc(gsl_odeiv_step_bsimp, *dim);
    break;
  case 9 : // M=1 implicit Gear method
    *stp = (void *)gsl_odeiv_step_alloc(gsl_odeiv_step_gear1, *dim);
    break;
  case 10 : // M=2 implicit Gear method
    *stp = (void *)gsl_odeiv_step_alloc(gsl_odeiv_step_gear2, *dim);
    break;
  }
}

void FC_FUNC_(gsl_odeiv_step_free, GSL_ODEIV_STEP_FREE)
     (void **stp)
{
  gsl_odeiv_step_free((gsl_odeiv_step *)(*stp));
}

void FC_FUNC_(gsl_odeiv_evolve_alloc, GSL_ODEIV_EVOLVE_ALLOC)
     (const int *dim, void **evl)
{
  *evl = (void *)gsl_odeiv_evolve_alloc (*dim);
}

int FC_FUNC_(gsl_odeiv_evolve_reset, GSL_ODEIV_EVOLVE_RESET)
     (void **evl)
{
  return gsl_odeiv_evolve_reset ((gsl_odeiv_evolve *)(*evl));
}

void FC_FUNC_(gsl_odeiv_evolve_free, GSL_ODEIV_EVOLVE_FREE)
     (void **evl)
{
  gsl_odeiv_evolve_free((gsl_odeiv_evolve *)(*evl));
}

void FC_FUNC_(gsl_odeiv_control_standart_new, GSL_ODEIV_CONTROL_STANDART_NEW)
     (void **ctrl, double *epsabs, double *epsrel, double *ay, double *adydt)
{
  *ctrl = (void*) gsl_odeiv_control_standard_new (*epsabs, *epsrel, *ay, *adydt);
}

void FC_FUNC_(gsl_odeiv_control_free, GSL_ODEIV_CONTROL_FREE)
     (void **ctrl)
{
  gsl_odeiv_control_free((gsl_odeiv_control *)(*ctrl));
}

          /* Interpolation */

void FC_FUNC_(gsl_interp_accel_alloc, GSL_INTERP_ACCEL_ALLOC)
     (void **acc)
{
     *acc = (void *)gsl_interp_accel_alloc();
}

void FC_FUNC_(gsl_spline_alloc, GSL_SPLINE_ALLOC)
     (void **spl, int *n, short int *opt)
{
  /* WARNING: This routine has a small modification: the user needs to
     pass an integer in order to choose the interpolation type instead
     of a variable of type gsl_interp_type. */

  switch (*opt) {
  case 1 : // linear interpolation
    *spl = (void *)gsl_spline_alloc(gsl_interp_linear, *n);
    break;
  case 2 : // polynomial interpolation
    *spl = (void *)gsl_spline_alloc(gsl_interp_polynomial, *n);
    break;
  case 3 : // cubic spline with natural boundary conditions
    *spl = (void *)gsl_spline_alloc(gsl_interp_cspline, *n);
    break;
  case 4 : // cubic spline with periodic boundary conditions
    *spl = (void *)gsl_spline_alloc(gsl_interp_cspline_periodic, *n);
    break;
  case 5 : // Akima spline with natural boundary conditions
    *spl = (void *)gsl_spline_alloc(gsl_interp_akima, *n);
    break;
  case 6 : // Akima spline with periodic boundary conditions
    *spl = (void *)gsl_spline_alloc(gsl_interp_akima_periodic, *n);
    break;
  }
}

void FC_FUNC_(gsl_spline_init, GSL_SPLINE_INIT)
     (void **spl, int *n, double *x, double *f)
{
  gsl_spline_init((gsl_spline *)(*spl), x, f, *n);
}

double FC_FUNC_(gsl_spline_eval, GSL_SPLINE_EVAL)
     (double *x, void **spl, void **acc)
{
  return gsl_spline_eval((gsl_spline *)(*spl), *x, (gsl_interp_accel *)(*acc));
}

double FC_FUNC_(gsl_spline_eval_deriv, GSL_SPLINE_EVAL_DERIV)
     (double *x, void **spl, void **acc)
{
  return gsl_spline_eval_deriv((gsl_spline *)(*spl), *x, (gsl_interp_accel *)(*acc));
}

double FC_FUNC_(gsl_spline_eval_deriv2, GSL_SPLINE_EVAL_DERIV2)
     (double *x, void **spl, void **acc)
{
  return gsl_spline_eval_deriv2((gsl_spline *)(*spl), *x, (gsl_interp_accel *)(*acc));
}

void FC_FUNC_(gsl_spline_free, GSL_SPLINE_FREE)
     (void **spl)
{
  gsl_spline_free((gsl_spline *)(*spl));

}

void FC_FUNC_(gsl_interp_accel_free, GSL_INTERP_ACCEL_FREE)
     (void **acc)
{
  gsl_interp_accel_free((gsl_interp_accel *)(*acc));
}

double FC_FUNC_(gsl_spline_eval_integ, GSL_SPLINE_EVAL_INTEG)
     (double *a, double *b, void **spl, void **acc)
{
  return gsl_spline_eval_integ((gsl_spline *)(*spl), *a, *b, (gsl_interp_accel *)(*acc));
}

          /* Multidimensional Root-Finding */

 void FC_FUNC_(gsl_multiroot_fsolver_alloc, GSL_MULTIROOT_FSOLVER_ALLOC)
     (void **s, const int *solver_type, const int *n)
{
  switch(*solver_type) {
  case 1 : // hybrid algorithm with internal scaling.
    *s = (void *)gsl_multiroot_fsolver_alloc(gsl_multiroot_fsolver_hybrids, *n);
    break;
  case 2 : // hybrid algorithm without internal scaling.
    *s = (void *)gsl_multiroot_fsolver_alloc(gsl_multiroot_fsolver_hybrid, *n);
    break;
  case 3 : // discrete Newton algorithm
    *s = (void *)gsl_multiroot_fsolver_alloc(gsl_multiroot_fsolver_dnewton, *n);
    break;
  case 4 : // Broyden algorithm
    *s = (void *)gsl_multiroot_fsolver_alloc(gsl_multiroot_fsolver_broyden, *n);
    break;
  }
}

void FC_FUNC_(gsl_multiroot_fsolver_free, GSL_MULTIROOT_FSOLVER_FREE)
     (void **s)
{
  gsl_multiroot_fsolver_free((gsl_multiroot_fsolver *)(*s));
}

max_stars_repo_name: cathcart/CCPG

max_issues_repo_name: cathcart/CCPG

max_forks_repo_name: cathcart/CCPG


text:
#include <jni.h>
#include "jade_reconstruct_area_PJNI.h"
#include <gsl/gsl_linalg.h>
#include <stdio.h>

JNIEXPORT jdoubleArray JNICALL Java_jade_reconstruct_area_PJNI_matrixExp
  (JNIEnv *jenv, jobject obj, jdoubleArray arr,jint size){
	jdoubleArray ret = (*jenv)->NewDoubleArray(jenv,(jint)( size * size));
	double *OutData = (*jenv)->GetDoubleArrayElements(jenv,ret,JNI_FALSE);
	//double OutData [size*size];
	jdouble * at_data = (*jenv)->GetDoubleArrayElements(jenv, arr, NULL);
	int i,j;
	double a_data [size*size];
	for (i = 0; i < size*size; i++){
		a_data[i] = (double)at_data[i];
	}
	gsl_matrix_view m = gsl_matrix_view_array (a_data, (int)size, (int)size);
	gsl_mode_t mt = 0;
	gsl_matrix *ma = gsl_matrix_alloc (size,size);
	gsl_linalg_exponential_ss(&m.matrix, ma, mt);//m input, ma output

	int x = 0;
	for(i = 0; i < size; i++){
		for(j=0;j<size;j++){
			OutData[x] = gsl_matrix_get(ma,i,j);
			x++;
		}
	}
	//gsl_permutation_free (ma);
 	//(*jenv)->SetDoubleArrayRegion(jenv,ret,(jsize)0,(jsize)size*size,OutData);
	(*jenv)->ReleaseDoubleArrayElements(jenv,ret,OutData,0);
	(*jenv)->ReleaseDoubleArrayElements(jenv,arr,at_data,0);
	gsl_matrix_free (ma);
	return ret;
}

ext: c


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

#include <math.h>
#include <complex.h>
#include <float.h>

#ifdef _MACOSX
#include <Accelerate/Accelerate.h>
#else
#include <cblas.h>
#endif

#include <gsl/gsl_sf_legendre.h>

#include "util.h"
#include "spbessel.h"

/* Use the modified Gram-Schmidt process to compute (in place) the portion of
 * the n-dimensional vector v orthogonal to each of the nv vectors s. The
 * projection * of the vector onto each of the basis vectors is stored in the
 * length-nv array c. */
int cmgs (complex double *v, complex double *c, complex double *s, int n, int nv) {
	int i, j, k;
	complex double *sv, cv;

	for (i = 0, sv = s; i < nv; ++i, sv += n) {
		c[i] = 0;
		k = 0;

		do {
			cv = pardot (sv, v, n);
			c[i] += cv;
#pragma omp parallel for default(shared) private(j)
			for (j = 0; j < n; ++j) v[j] -= cv * sv[j];
		} while (cabs(cv / c[i]) > IMGS_TOL && ++k < IMGS_ITS);

	}

	return n;
}

complex double pardot (complex double *x, complex double *y, int n) {
	complex double dp;

	/* Compute the local portion. */
	cblas_zdotc_sub (n, x, 1, y, 1, &dp);

	return dp;
}

/* The MSE between two vectors. */
double rmserror (complex double *v, complex double *r, int n) {
	double err = 0, errd = 0, e;
	int i;

	for (i = 0; i < n; ++i) {
		e = cabs (v[i] - r[i]);
		err += e * e;
		e = cabs (r[i]);
		errd += e * e;
	}

	return sqrt (err / errd);
}

/* Open a file or die. */
FILE *critopen (char *fname, char *mode) {
	FILE *fptr;

	fptr = fopen (fname, mode);

	if (!fptr) {
		fprintf (stderr, "ERROR: Could not open input file %s.\n", fname);
		exit (EXIT_FAILURE);

		return NULL;
	}

	return fptr;
}

/* Compute the number of harmonics required using the excess bandwidth
 * formula. */
int exband (complex double kr, double ndig) {
	double akr;
	int l;

	akr = creal (kr);
	ndig *= ndig;

	l = ceil (akr + 1.8 * cbrt(ndig * akr));

	return l;
}

/* Computes the Legendre polynomials up to order n for the argument x. The
 * values are stored in the array v. */
int legpoly (int n, double x, double *v) {
	int i;

	/* Don't bother computing anything for order less than zero. */
	if (n < 0) return -1;

	/* Make sure the argument is within [-1,1] as expected. */
	if (fabs(x) > 1.0 + DBL_EPSILON) return -2;

	/* The first function value. */
	v[0] = 1.0;

	/* If the order happens to be zero, there is no going further. */
	if (n < 1) return 0;

	/* The second function value. */
	v[1] = x;

	/* The recursion formula for Legendre polynomials. */
	for (i = 1; i < n - 1; ++i)
		v[i + 1] = ((2 * i + 1) * x * v[i] - i * v[i - 1]) / (i + 1);

	return 0;
}

/* Computes the wave number from the relative sound speed and
 * unitless (dB) attenuation coefficient. */
complex double wavenum (double cr, double alpha) {
	double kr, ki;

	ki = log(10) * alpha / 20;
	kr = 2 * M_PI / cr;

	return kr + I * ki;
}

/* Computes the inverse wave number from the relative sound speed and
 * unitless (dB) attenuation coefficient. */
complex double invwavenum (double cr, double alpha) {
	double kr, ki, mag, csq;

	csq = cr * cr;

	kr = 2 * M_PI * cr;
	ki = log(10) * alpha / 20;
	mag = 4 * M_PI * M_PI + ki * ki * csq;

	return (kr / mag) - I * (ki * csq / mag);
}

/* Copy the spherical harmonic representation from one location to another. */
int copysh (int deg, complex double *out, int ldo, complex double *in, int ldi) {
	int i, j, offo, offi;

	/* Copy the spherical harmonic coefficients into the right place. */
	for (i = 0; i < deg; ++i) {
		offo = i * ldo;
		offi = i * ldi;
		/* Copy the zero-order coefficients. */
		out[offo] = in[offi];

		/* Copy the other coefficients. */
		for (j = 1; j <= i; ++j) {
			out[offo + j] = in[offi + j];
			out[offo + ldo - j] = in[offi + ldi - j];
		}
	}

	return deg;
}

/* Multiply in a radial component to the spherical harmonic. */
int shradial (int deg, complex double *a, int lda, complex double k, double r) {
	complex double kr, *hlkr;
	int i, j, off;

	hlkr = malloc (deg * sizeof(complex double));

	/* Compute the Hankel functions for all degrees. */
	kr = k * r;
	spbesh (hlkr, kr, deg);

	for (i = 0; i < deg; ++i) {
		off = i * lda;

		/* The zero-order coefficients. */
		a[off] *= hlkr[i];

		/* The other coefficients. */
		for (j = 1; j <= i; ++j) {
			a[off + j] *= hlkr[i];
			a[off + lda - j] *= hlkr[i];
		}
	}

	free (hlkr);

	return deg;
}

/* Compute the harmonic coefficients for an incident plane wave. */
int shincident (int deg, complex double *a, int lda,
		complex double mag, double theta, double phi) {
	double cth, *lgvals;
	int l, m, dm1, npm, off;
	long lgi, lgn;
	complex double scale[4] = { 1, -I, -1, I }, cx, fx;

	/* The coefficient has a 4 pi factor that must be included. */
	mag *= 4 * M_PI;

	dm1 = deg - 1;
	cth = cos(theta);

	lgn = gsl_sf_legendre_array_n (dm1);
	lgvals = malloc (lgn * sizeof(double));

	/* Compute the Legendre polynomials of zero order for all degrees. */
	gsl_sf_legendre_array_e (GSL_SF_LEGENDRE_SPHARM, dm1, cth, 1., lgvals);

	/* Compute the zero-order coefficients for all degrees. */
	for (l = 0; l < deg; ++l) {
		lgi = gsl_sf_legendre_array_index(l, 0);
		a[l * lda] += scale[l % 4] * mag * lgvals[lgi];
	}

	for (m = 1; m < deg; ++m) {
		npm = lda - m;

		/* Compute the phi variation. */
		cx = cexp (I * m * phi);

#pragma omp critical(incaug)
		for (l = m; l < deg; ++l) {
			off = l * lda;
			lgi = gsl_sf_legendre_array_index(l, m);
			fx = scale[l % 4] * mag * lgvals[lgi];
			/* The positive-order coefficient. */
			a[off + m] += fx * conj (cx);
			/* The negative-order coefficient. */
			a[off + npm] += fx * cx;
		}
	}

	free (lgvals);

	return deg;
}

static int lgval (double *p, double *dp, double t, int m) {
	double p0 = 1.0, p1 = t;
	int k;

	/* Set values explicitly for low orders. */
	if (m < 1) {
		*p = 1.0;
		*dp = 0.0;
		return m;
	} else if (m < 2) {
		*p = t;
		*dp = 1.0;
		return m;
	}

	/* Otherwise compute the values explicitly. */
	for (k = 1; k < m; ++k) {
		*p = ((2.0 * k + 1.0) * t * p1 - k * p0) / (1.0 + k);
		p0 = p1; p1 = *p;
	}

	*dp = m * (p0 - t * p1) / (1.0 - t * t);

	return m;
}

int gauleg (int m, double *nodes, double *weights) {
	int i, j, nroots = (m + 1) / 2;
	double t, p, dp, dt;
	const int maxit = 100;
	const double tol = 1.0e-14;

#pragma omp parallel for default(shared) private(i,j,t,p,dp,dt)
	for (i = 0; i < nroots; ++i) {
		t = cos (M_PI * (i + 0.75) / (m + 0.5));

		for (j = 0; j < maxit; ++j) {
			/* Compute the value of the Legendre polynomial. */
			lgval (&p, &dp, t, m);

			/* Perform a Newton-Raphson update. */
			dt = -p / dp;
			t += dt;

			/* Break if convergence has been achieved. */
			if (fabs(dt) < tol) break;
		}

		/* Update the nodes and weights. */
		nodes[i] = t;
		nodes[m - i - 1] = -t;

		weights[i] = 2.0 / (1.0 - t * t) / (dp * dp);
		weights[m - i - 1] = weights[i];
	}

	return 0;
}

max_stars_repo_name: ahesford/fastsphere

max_issues_repo_name: ahesford/fastsphere

max_forks_repo_name: ahesford/fastsphere


text:
/* -*- Mode: C++; tab-width: 4; c-basic-offset: 4; indent-tabs-mode: nil -*- */
/*
 *     Copyright 2016-Present Couchbase, Inc.
 *
 *   Use of this software is governed by the Business Source License included
 *   in the file licenses/BSL-Couchbase.txt.  As of the Change Date specified
 *   in that file, in accordance with the Business Source License, use of this
 *   software will be governed by the Apache License, Version 2.0, included in
 *   the file licenses/APL2.txt.
 */
#pragma once

#include <gsl/gsl-lite.hpp>
#include <nlohmann/json_fwd.hpp>
#include <cstdint>
#include <string>

namespace cb::crypto {
enum class Algorithm { MD5, SHA1, SHA256, SHA512 };

bool isSupported(Algorithm algorithm);

const int MD5_DIGEST_SIZE = 16;
const int SHA1_DIGEST_SIZE = 20;
const int SHA256_DIGEST_SIZE = 32;
const int SHA512_DIGEST_SIZE = 64;

/**
 * Generate a HMAC digest of the key and data by using the given
 * algorithm
 *
 * @throws std::invalid_argument - unsupported algorithm
 *         std::runtime_error - Failures generating the HMAC
 */
std::string HMAC(Algorithm algorithm,
                 std::string_view key,
                 std::string_view data);

/**
 * Generate a PBKDF2_HMAC digest of the key and data by using the given
 * algorithm
 *
 * @throws std::invalid_argument - unsupported algorithm
 *         std::runtime_error - Failures generating the HMAC
 */
std::string PBKDF2_HMAC(Algorithm algorithm,
                        const std::string& pass,
                        std::string_view salt,
                        unsigned int iterationCount);

/**
 * Generate a digest by using the requested algorithm
 */
std::string digest(Algorithm algorithm, std::string_view data);

enum class Cipher { AES_256_cbc };

Cipher to_cipher(const std::string& str);

/**
 * Encrypt the specified data by using a given cipher
 *
 * @param cipher The cipher to use
 * @param key The key used for encryption
 * @param ivec The IV to use for encryption
 * @param data The Pointer to the data to encrypt
 * @return The encrypted data
 */
std::string encrypt(Cipher cipher,
                    std::string_view key,
                    std::string_view iv,
                    std::string_view data);

/**
 * Encrypt the specified data
 *
 * {
 *    "cipher" : "name of the cipher",
 *    "key" : "base64 of the raw key",
 *    "iv" : "base64 of the raw iv"
 * }
 *
 * @param meta the json description of the encryption to use
 * @param data Pointer to the data to encrypt
 * @param length the length of the data to encrypt
 * @return The encrypted data
 */
std::string encrypt(const nlohmann::json& json, std::string_view data);

/**
 * Decrypt the specified data by using a given cipher
 *
 * @param cipher The cipher to use
 * @param key The key used for encryption
 * @param ivec The IV to use for encryption
 * @param data The data to decrypt
 * @return The decrypted data
 */
std::string decrypt(Cipher cipher,
                    std::string_view key,
                    std::string_view iv,
                    std::string_view data);

/**
 * Decrypt the specified data
 *
 * {
 *    "cipher" : "name of the cipher",
 *    "key" : "base64 of the raw key",
 *    "iv" : "base64 of the raw iv"
 * }
 *
 * @param meta the json description of the encryption to use
 * @param data Pointer to the data to decrypt
 * @param length the length of the data to decrypt
 * @return The decrypted data
 */
std::string decrypt(const nlohmann::json& json, std::string_view data);

} // namespace cb::crypto

ext: h


text:
#ifndef AVERAGER_H_
#define AVERAGER_H_

#include <gsl/gsl_statistics_double.h>
#include <string.h>
#include <iostream>

class Averager
{
public:
    virtual ~Averager() {};
    virtual void reset() = 0;
    virtual void feed(const double x) = 0;
    virtual double get_mean() const = 0;
    virtual double get_var() const = 0;
    virtual unsigned int get_n() const = 0;
};

class SimpleAverager: public Averager
{
public:
    SimpleAverager():
        _n(0), _sum(0), _sum_of_squares(0) {}
    virtual void reset() { _n = 0; _sum = 0; _sum_of_squares = 0; }
    virtual void feed(const double x)
    {
        _n++; _sum += x;
        _sum_of_squares += x * x;
    }
    virtual double get_mean() const { return _sum / _n; }
    virtual double get_var() const
    {
        return 1./(_n - 1) * (_sum_of_squares - _n * get_mean()*get_mean());
    }
    virtual unsigned int get_n() const { return _n; }
private:
    unsigned int _n;
    double _sum;
    double _sum_of_squares;
};

class BufferAverager: public Averager
{
public:
    BufferAverager(const size_t bufsize):
        bufsize_(bufsize), n_(0), i_(0)
    {
        buf_ = new double[bufsize];
    }
    virtual ~BufferAverager()
    {
        delete[] buf_;
    }
    virtual void feed(const double x)
    {
        if (i_ == bufsize_)
        {
            i_ = 0;
        }
        if (i_ < bufsize_)
        {
            buf_[i_] = x;
            i_++;
            if (i_ > n_)
                n_++;
        }
    }
    virtual double get_mean() const
    {
        return gsl_stats_mean(buf_, 1, n_);
    }
    virtual double get_var() const
    {
        return gsl_stats_variance(buf_, 1, n_);
    }
    virtual unsigned int get_n() const
    {
        return n_;
    }
    virtual void reset()
    {
        memset(buf_, 0, bufsize_ * sizeof(double));
        n_ = 0;
        i_ = 0;
    }
    void dump(std::ostream& out)
    {
        for (int j = 0; j < n_; j++)
            out << buf_[j] << "\n";
    }

private:
    const size_t bufsize_;
    double* buf_;
    unsigned int n_;
    unsigned int i_;
};

class FancyAverager: public Averager
{
public:
    FancyAverager():
        _mean_avg(), _sq_avg(), _cub_avg() {}
    virtual void reset()
    {
        _mean_avg.reset();
        _sq_avg.reset();
    }
    virtual void feed(const double x)
    {
        _mean_avg.feed(x);
        _sq_avg.feed(x*x);
        _cub_avg.feed(x*x*x);
    }
    virtual unsigned int get_n() const { return _mean_avg.get_n(); }
    virtual double get_mean() const { return _mean_avg.get_mean(); }
    virtual double get_var() const { return _mean_avg.get_var(); }
    double get_mean_of_square() const { return _sq_avg.get_mean(); }
    double get_var_of_square() const { return _sq_avg.get_var(); }
    double get_mean_of_cube() const { return _cub_avg.get_mean(); }
    double get_var_of_cube() const { return _cub_avg.get_var(); }
private:
    SimpleAverager _mean_avg;
    SimpleAverager _sq_avg;
    SimpleAverager _cub_avg;
};

#endif // AVERAGER_H_


text:
/* vector/test.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Gerard Jungman, Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>

#if defined( _MSC_VER ) && defined( GSL_DLL )
#undef inline
#define inline __forceinline
#endif

#if (!GSL_RANGE_CHECK) && defined(HAVE_INLINE)
#undef GSL_RANGE_CHECK
#define GSL_RANGE_CHECK 1
#endif

#include <stdlib.h>
#include <unistd.h>
#include <fcntl.h>
#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>

int status = 0;

#ifndef DESC
#define DESC ""
#endif

#define BASE_GSL_COMPLEX_LONG
#include "templates_on.h"
#include "test_complex_source.c"
#include "templates_off.h"
#undef  BASE_GSL_COMPLEX_LONG

#define BASE_GSL_COMPLEX
#include "templates_on.h"
#include "test_complex_source.c"
#include "templates_off.h"
#undef  BASE_GSL_COMPLEX

#define BASE_GSL_COMPLEX_FLOAT
#include "templates_on.h"
#include "test_complex_source.c"
#include "templates_off.h"
#undef  BASE_GSL_COMPLEX_FLOAT

#define BASE_LONG_DOUBLE
#include "templates_on.h"
#include "test_source.c"
#include "templates_off.h"
#undef  BASE_LONG_DOUBLE

#define BASE_DOUBLE
#include "templates_on.h"
#include "test_source.c"
#include "templates_off.h"
#undef  BASE_DOUBLE

#define BASE_FLOAT
#include "templates_on.h"
#include "test_source.c"
#include "templates_off.h"
#undef  BASE_FLOAT

#define BASE_ULONG
#include "templates_on.h"
#include "test_source.c"
#include "templates_off.h"
#undef  BASE_ULONG

#define BASE_LONG
#include "templates_on.h"
#include "test_source.c"
#include "templates_off.h"
#undef  BASE_LONG

#define BASE_UINT
#include "templates_on.h"
#include "test_source.c"
#include "templates_off.h"
#undef  BASE_UINT

#define BASE_INT
#include "templates_on.h"
#include "test_source.c"
#include "templates_off.h"
#undef  BASE_INT

#define BASE_USHORT
#include "templates_on.h"
#include "test_source.c"
#include "templates_off.h"
#undef  BASE_USHORT

#define BASE_SHORT
#include "templates_on.h"
#include "test_source.c"
#include "templates_off.h"
#undef  BASE_SHORT

#define BASE_UCHAR
#include "templates_on.h"
#include "test_source.c"
#include "templates_off.h"
#undef  BASE_UCHAR

#define BASE_CHAR
#include "templates_on.h"
#include "test_source.c"
#include "templates_off.h"
#undef  BASE_CHAR

void my_error_handler (const char *reason, const char *file,
                       int line, int err);

int
main (void)
{
  size_t stride, ostride, N;

  gsl_ieee_env_setup ();

  for (N = 10; N < 1024; N = 2*N + 1)
    {
      for (stride = 1; stride < 5 ; stride++)
        {
          test_func (stride, N);
          test_float_func (stride, N);
          test_long_double_func (stride, N);
          test_ulong_func (stride, N);
          test_long_func (stride, N);
          test_uint_func (stride, N);
          test_int_func (stride, N);
          test_ushort_func (stride, N);
          test_short_func (stride, N);
          test_uchar_func (stride, N);
          test_char_func (stride, N);

          test_complex_func (stride, N);
          test_complex_float_func (stride, N);
          test_complex_long_double_func (stride, N);

          for (ostride = 1; ostride < 5 ; ostride++)
            {
              test_ops (stride, ostride, N);
              test_float_ops (stride, ostride, N);
              test_long_double_ops (stride, ostride, N);
              test_ulong_ops (stride, ostride, N);
              test_long_ops (stride, ostride, N);
              test_uint_ops (stride, ostride, N);
              test_int_ops (stride, ostride, N);
              test_ushort_ops (stride, ostride, N);
              test_short_ops (stride, ostride, N);
              test_uchar_ops (stride, ostride, N);
              test_char_ops (stride, ostride, N);
              test_complex_ops (stride, ostride, N);
              test_complex_float_ops (stride, ostride, N);
              test_complex_long_double_ops (stride, ostride, N);
            }

          test_text (stride, N);
          test_float_text (stride, N);
#if HAVE_PRINTF_LONGDOUBLE
          test_long_double_text (stride, N);
#endif
          test_ulong_text (stride, N);
          test_long_text (stride, N);
          test_uint_text (stride, N);
          test_int_text (stride, N);
          test_ushort_text (stride, N);
          test_short_text (stride, N);
          test_uchar_text (stride, N);
          test_char_text (stride, N);

          test_complex_text (stride, N);
          test_complex_float_text (stride, N);
#if HAVE_PRINTF_LONGDOUBLE
          test_complex_long_double_text (stride, N);
#endif

          test_file (stride, N);
          test_float_file (stride, N);
          test_long_double_file (stride, N);
          test_ulong_file (stride, N);
          test_long_file (stride, N);
          test_uint_file (stride, N);
          test_int_file (stride, N);
          test_ushort_file (stride, N);
          test_short_file (stride, N);
          test_uchar_file (stride, N);
          test_char_file (stride, N);
          test_complex_file (stride, N);
          test_complex_float_file (stride, N);
          test_complex_long_double_file (stride, N);
        }
    }

#if GSL_RANGE_CHECK
  gsl_set_error_handler (&my_error_handler);

  for (N = 1; N < 1024; N *=2)
    {
      for (stride = 1; stride < 5 ; stride++)
        {
          test_trap (stride, N);
          test_float_trap (stride, N);
          test_long_double_trap (stride, N);
          test_ulong_trap (stride, N);
          test_long_trap (stride, N);
          test_uint_trap (stride, N);
          test_int_trap (stride, N);
          test_ushort_trap (stride, N);
          test_short_trap (stride, N);
          test_uchar_trap (stride, N);
          test_char_trap (stride, N);
          test_complex_trap (stride, N);
          test_complex_float_trap (stride, N);
          test_complex_long_double_trap (stride, N);
        }
    }
#endif

  exit (gsl_test_summary ());
}

void
my_error_handler (const char *reason, const char *file, int line, int err)
{
  if (0)
    printf ("(caught [%s:%d: %s (%d)])\n", file, line, reason, err);
  status = 1;
}

max_stars_repo_name: utdsimmons/ohpc

max_issues_repo_name: utdsimmons/ohpc

max_forks_repo_name: utdsimmons/ohpc


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

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

#include "ccl.h"
#include "ccl_params.h"

// Global variable to hold the neutrino phase-space spline
gsl_spline* nu_spline=NULL;

/* ------- ROUTINE: nu_integrand ------
INPUTS: x: dimensionless momentum, *r: pointer to a dimensionless mass / temperature
TASK: Integrand of phase-space massive neutrino integral
*/
static double nu_integrand(double x, void *r)
{
  double rat=*((double*)(r));
  return sqrt(x*x+rat*rat)/(exp(x)+1.0)*x*x;
}

/* ------- ROUTINE: ccl_calculate_nu_phasespace_spline ------
TASK: Get the spline of the result of the phase-space integral required for massive neutrinos.
*/

gsl_spline* calculate_nu_phasespace_spline(int *status) {
  int N=CCL_NU_MNUT_N;
  double *mnut = ccl_linear_spacing(log(CCL_NU_MNUT_MIN),log(CCL_NU_MNUT_MAX),N);
  double *y=malloc(sizeof(double)*CCL_NU_MNUT_N);
  if (y ==NULL) {
    // Not setting a status_message here because we can't easily pass a cosmology to this function - message printed in ccl_error.c.
    *status = CCL_ERROR_NU_INT;
  }

  // Check whether ccl_splines and ccl_gsl exist; exit gracefully if they
  // can't be loaded
  if(ccl_splines==NULL || ccl_gsl==NULL) ccl_cosmology_read_config();
  if(ccl_splines==NULL || ccl_gsl==NULL) {
    ccl_raise_exception(CCL_ERROR_MISSING_CONFIG_FILE,
                        "ccl_neutrinos.c: Failed to read config file.");
    *status = CCL_ERROR_MISSING_CONFIG_FILE;
    return NULL;
  }

  int stat=0, gslstatus;
  gsl_integration_cquad_workspace * workspace = \
                gsl_integration_cquad_workspace_alloc(ccl_gsl->N_ITERATION);
  gsl_function F;
  F.function = &nu_integrand;
  for (int i=0; i<CCL_NU_MNUT_N; i++) {
    double mnut_=exp(mnut[i]);
    F.params = &(mnut_);
    gslstatus = gsl_integration_cquad(&F, 0, 1000.0,
                                      ccl_gsl->INTEGRATION_NU_EPSABS,
                                      ccl_gsl->INTEGRATION_NU_EPSREL,
                                      workspace,&y[i], NULL, NULL);
    if(gslstatus != GSL_SUCCESS) {
      ccl_raise_gsl_warning(gslstatus, "ccl_neutrinos.c: calculate_nu_phasespace_spline():");
      stat |= gslstatus;
    }
  }
  gsl_integration_cquad_workspace_free(workspace);
  double renorm=1./y[0];
  for (int i=0; i<CCL_NU_MNUT_N; i++)
    y[i]*=renorm;
  gsl_spline* spl = gsl_spline_alloc(A_SPLINE_TYPE, CCL_NU_MNUT_N);
  stat |= gsl_spline_init(spl, mnut, y, CCL_NU_MNUT_N);

  // Check for errors in creating the spline
  if (stat) {
    // Not setting a status_message here because we can't easily pass a cosmology to this function - message printed in ccl_error.c.
    *status = CCL_ERROR_NU_INT;
    free(mnut);
    free(y);
    gsl_spline_free(spl);
  }
  free(mnut);
  free(y);
  return spl;
}

/* ------- ROUTINE: ccl_nu_phasespace_intg ------
INPUTS: accel: pointer to an accelerator which will evaluate the neutrino phasespace spline if defined, mnuOT: the dimensionless mass / temperature of a single massive neutrino
TASK: Get the value of the phase space integral at mnuOT
*/

double nu_phasespace_intg(gsl_interp_accel* accel, double mnuOT, int* status)
{
  // Check if the global variable for the phasespace spline has been defined yet:
  if (nu_spline==NULL) nu_spline =calculate_nu_phasespace_spline(status);
  ccl_check_status_nocosmo(status);

  double integral_value =0.;

  // First check the cases where we are in the limits.
  if (mnuOT<CCL_NU_MNUT_MIN) {
    return 7./8.;
  }
  else if (mnuOT>CCL_NU_MNUT_MAX) {
    return 0.2776566337*mnuOT;
  }

  // Evaluate the spline - this will use the accelerator if it has been defined.
  int gslstatus = gsl_spline_eval_e(nu_spline, log(mnuOT),accel, &integral_value);
  if(gslstatus != GSL_SUCCESS) {
    ccl_raise_gsl_warning(gslstatus, "ccl_neutrinos.c: nu_phasespace_intg():");
    *status |= gslstatus;
  }
  return integral_value*7./8.;
}

/* -------- ROUTINE: Omeganuh2 ---------
INPUTS: a: scale factor, Nnumass: number of massive neutrino species, mnu: total mass in eV of neutrinos, T_CMB: CMB temperature, accel: pointer to an accelerator which will evaluate the neutrino phasespace spline if defined, status: pointer to status integer.
TASK: Compute Omeganu * h^2 as a function of time.
!! To all practical purposes, Neff is simply N_nu_mass !!
*/

double ccl_Omeganuh2 (double a, int N_nu_mass, double* mnu, double T_CMB, gsl_interp_accel* accel, int* status)
{
  double Tnu, a4, prefix_massless, mnuone, OmNuh2;
  double Tnu_eff, mnuOT, intval, prefix_massive;
  double total_mass; // To check if this is the massless or massive case.

  // First check if N_nu_mass is 0
  if (N_nu_mass==0) return 0.0;

  Tnu=T_CMB*pow(4./11.,1./3.);
  a4=a*a*a*a;
  // Check if mnu=0. We assume that in the massless case mnu is a pointer to a single element and that element is 0. This should in principle never be called.
  if (mnu[0] < 0.00017) {  // Limit taken from Lesgourges et al. 2012
    prefix_massless = NU_CONST  * Tnu * Tnu * Tnu * Tnu;
    return N_nu_mass*prefix_massless*7./8./a4;
  }

  // And the remaining massive case.
  // Tnu_eff is used in the massive case because CLASS uses an effective temperature of nonLCDM components to match to mnu / Omeganu =93.14eV. Tnu_eff = T_ncdm * T_CMB = 0.71611 * T_CMB
  Tnu_eff = Tnu * TNCDM / (pow(4./11.,1./3.));

  // Define the prefix using the effective temperature (to get mnu / Omega = 93.14 eV) for the massive case:
  prefix_massive = NU_CONST * Tnu_eff * Tnu_eff * Tnu_eff * Tnu_eff;

  OmNuh2 = 0.; // Initialize to 0 - we add to this for each massive neutrino species.
  for(int i=0; i<N_nu_mass; i++){
	// Get mass over T (mass (eV) / ((kb eV/s/K) Tnu_eff (K))
	// This returns the density normalized so that we get nuh2 at a=0
	mnuOT = mnu[i] / (Tnu_eff/a) * (EV_IN_J / (KBOLTZ));

	// Get the value of the phase-space integral
	intval=nu_phasespace_intg(accel,mnuOT, status);
	OmNuh2 = intval*prefix_massive/a4 + OmNuh2;
  }

  return OmNuh2;
}

/* -------- ROUTINE: Omeganuh2_to_Mnu ---------
INPUTS: OmNuh2: neutrino mass density today Omeganu * h^2, label: how you want to split up the masses, see ccl_neutrinos.h for options, T_CMB: CMB temperature, accel: pointer to an accelerator which will evaluate the neutrino phasespace spline if defined, status: pointer to status integer.
TASK: Given Omeganuh2 today, the method of splitting into masses, and the temperature of the CMB, output a pointer to the array of neutrino masses (may be length 1 if label asks for sum)
*/

double* ccl_nu_masses(double OmNuh2, ccl_neutrino_mass_splits mass_split, double T_CMB,  int* status){

  double sumnu;

  sumnu = 93.14 * OmNuh2;

  // Now split the sum up into three masses depending on the label given:

  if(mass_split==ccl_nu_normal){

     double *mnu;
	 mnu = malloc(3*sizeof(double));

     // See CCL note for how we get these expressions for the neutrino masses in normal and inverted hierarchy.
	 mnu[0] = 2./3.* sumnu - 1./6. * pow(-6. * DELTAM12_sq + 12. * DELTAM13_sq_pos + 4. * sumnu*sumnu, 0.5)
	         - 0.25 * DELTAM12_sq / (2./3.* sumnu - 1./6. * pow(-6. * DELTAM12_sq + 12. * DELTAM13_sq_pos + 4. * sumnu*sumnu, 0.5));
	 mnu[1] = 2./3.* sumnu - 1./6. * pow(-6. * DELTAM12_sq + 12. * DELTAM13_sq_pos + 4. * sumnu*sumnu, 0.5)
	         + 0.25 * DELTAM12_sq / (2./3.* sumnu - 1./6. * pow(-6. * DELTAM12_sq + 12. * DELTAM13_sq_pos + 4. * sumnu*sumnu, 0.5));
	 mnu[2] = -1./3. * sumnu + 1./3 * pow(-6. * DELTAM12_sq + 12. * DELTAM13_sq_pos + 4. * sumnu*sumnu, 0.5);
	 if (mnu[0]<0 || mnu[1]<0 || mnu[2]<0){
	    // The user has provided a sum that is below the physical limit.
	    if (sumnu<1e-14){
			mnu[0] = 0.;
			mnu[1] = 0.;
			mnu[2] = 0.;
	    }else{
		 *status = CCL_ERROR_MNU_UNPHYSICAL;
	 }
     }
	 ccl_check_status_nocosmo(status);
	 return mnu;

  } else if (mass_split==ccl_nu_inverted){

	double *mnu;
	mnu = malloc(3*sizeof(double));
	mnu[0] = 2./3.* sumnu - 1./6. * pow(-6. * DELTAM12_sq + 12. * DELTAM13_sq_neg + 4. * sumnu * sumnu, 0.5)
	         - 0.25 * DELTAM12_sq / (2./3.* sumnu - 1./6. * pow(-6. * DELTAM12_sq + 12. * DELTAM13_sq_neg + 4. * sumnu * sumnu, 0.5));
	mnu[1] = 2./3.* sumnu - 1./6. * pow(-6. * DELTAM12_sq + 12. * DELTAM13_sq_neg + 4. * sumnu*sumnu, 0.5)
	         + 0.25 * DELTAM12_sq / (2./3.* sumnu - 1./6. * pow(-6. * DELTAM12_sq + 12. * DELTAM13_sq_neg + 4. * sumnu*sumnu, 0.5));
	mnu[2] = -1./3. * sumnu + 1./3 * pow(-6. * DELTAM12_sq + 12. * DELTAM13_sq_neg + 4. * sumnu*sumnu, 0.5);
	if(mnu[0]<0 || mnu[1]<0 || mnu[2]<0){
	// The user has provided a sum that is below the physical limit.
	    if (sumnu<1e-14){
			mnu[0] = 0.;
			mnu[1] = 0.;
			mnu[2] = 0.;;
        }else{
		*status = CCL_ERROR_MNU_UNPHYSICAL;
	    }
	}

	ccl_check_status_nocosmo(status);
	return mnu;

  } else if (mass_split==ccl_nu_equal){
	  double *mnu;
	  mnu = malloc(3*sizeof(double));
	  mnu[0] = sumnu/3.;
	  mnu[1] = sumnu/3.;
	  mnu[2] = sumnu/3.;

      return mnu;

  } else if (mass_split == ccl_nu_sum){

	  double *mnu;
	  mnu = malloc(sizeof(double));

      mnu[0] = sumnu;

      return mnu;

  } else{

	  printf("WARNING:  mass option = %d not yet supported\n continuing with normal hierarchy\n", mass_split);
	  double *mnu;
	  mnu = malloc(3*sizeof(double));

	 if (sumnu>0.59){
		 mnu[0] = (sumnu - 0.59) / 3.;
		 mnu[1] = mnu[0] + 0.009;
		 mnu[2] = mnu[0] + 0.05;
	 }else{
		 *status = CCL_ERROR_MNU_UNPHYSICAL;
	 }

	 ccl_check_status_nocosmo(status);
	 return mnu;
  }
}

max_stars_repo_name: Russell-Jones-OxPhys/CCL

max_issues_repo_name: Russell-Jones-OxPhys/CCL

max_forks_repo_name: Russell-Jones-OxPhys/CCL


text:
#pragma once

#include <gsl/span>
#include <utility>

namespace hz::range {
    template<typename T>
    using contiguous_view = gsl::span<T>;

    template<typename Container>
    auto make_contiguous_view(Container && c) noexcept -> contiguous_view<std::remove_reference_t<decltype(std::data(*std::forward<Container>(c)))>> {
        return gsl::make_span(std::forward<Container>(c));
    }
}

max_stars_repo_name: KABoissonneault/AGEA

max_issues_repo_name: KABoissonneault/AGEA

max_forks_repo_name: KABoissonneault/AGEA


text:
/*
* -----------------------------------------------------------------
*  Utilities Library --- util_lib.c
*  Version: 1.6180
*  Date: Mar 12, 2010
* -----------------------------------------------------------------
*  Programmer: Americo Barbosa da Cunha Junior
*              americo.cunhajr@gmail.com
* -----------------------------------------------------------------
*  Copyright (c) 2010 by Americo Barbosa da Cunha Junior
*
*  This program is free software: you can redistribute it and/or
*  modify it under the terms of the GNU General Public License as
*  published by the Free Software Foundation, either version 3 of
*  the License, or (at your option) any later version.
*
*  This program is distributed in the hope that it will be useful,
*  but WITHOUT ANY WARRANTY; without even the implied warranty of
*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
*  GNU General Public License for more details.
*
*  A copy of the GNU General Public License is available in
*  LICENSE.txt or http://www.gnu.org/licenses/.
* -----------------------------------------------------------------
*  This is the implementation file of a library
*  with miscelaneous functions.
* -----------------------------------------------------------------
*/

#include <stdio.h>
#include <gsl/gsl_errno.h>

#include "../include/util_lib.h"

/*
*------------------------------------------------------------
*   gsl_fprint_vector
*
*   This function prints a vector in a file.
*
*   Input:
*   file - output file
*   v    - GSL vector
*
*   Output:
*
*   Return:
*   success or error
*
*   last update: Nov 6, 2008
*------------------------------------------------------------
*/

int gsl_fprint_vector(FILE* file,gsl_vector *v)
{
    if( v == NULL )
	return GSL_EINVAL;

    unsigned int i;

    for( i = 0; i < v->size; i++ )
	fprintf(file," %+.3e", v->data[i]);

    return GSL_SUCCESS;
}
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   gsl_print_vector
*
*   This function prints a vector on the screen
*
*   Input:
*   v  - GSL vector
*
*   Output:
*
*   Return:
*   success or error
*
*   last update: Dec 26, 2008
*------------------------------------------------------------
*/

int gsl_print_vector(gsl_vector *v)
{
    if( v == NULL )
	return GSL_EINVAL;

    unsigned int i;

    for( i = 0; i < v->size; i++ )
	printf(" %+.3e", v->data[i]);

    return GSL_SUCCESS;
}
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   gsl_print_matrix_line
*
*   This function prints a matrix line on the screen.
*
*   Input:
*   i  - matrix line to be printed
*   M  - GSL matrix
*
*   Output:
*
*   Return:
*   success or error
*
*   last update: Dec 26, 2008
*------------------------------------------------------------
*/

int gsl_print_matrix_line(unsigned int i,gsl_matrix *M)
{
    if( M == NULL )
	return GSL_EINVAL;

    unsigned int j;

    printf("\n");
    for( j = 0; j < M->size2; j++ )
	printf(" %+.3e ", M->data[i*M->tda+j]);

    return GSL_SUCCESS;
}
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   gsl_print_matrix_column
*
*   This function prints a matrix column on the screen
*
*   Input:
*   M  - GSL matrix
*   i  - matrix column to be printed
*
*   Output:
*
*   Return:
*   success or error
*
*   last update: Dec 26, 2008
*------------------------------------------------------------
*/

int gsl_print_matrix_column(unsigned int j,gsl_matrix *M)
{
    if( M == NULL )
    	return GSL_EINVAL;

    unsigned int i;

    printf("\n");
    for( i = 0; i < M->size1; i++ )
	printf(" %+.3e\n", M->data[i*M->tda+j]);

    return GSL_SUCCESS;
}
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   gsl_print_matrix
*
*   This function prints a matrix on the screen
*
*   Input:
*   M  - GSL matrix
*
*   Outpur:
*
*   Return:
*   success or error
*
*   last update: Dec 26, 2008
*------------------------------------------------------------
*/

int gsl_print_matrix(gsl_matrix *M)
{
    if( M == NULL )
    	return GSL_EINVAL;

    unsigned int i;

    for( i = 0; i < M->size1; i++ )
	gsl_print_matrix_line(i,M);

    return GSL_SUCCESS;
}
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   gsl_rand_vector
*
*   This function creates a random vector.
*
*   Input:
*   seed - random number generator seed
*
*   Output:
*   v  - GSL vector
*
*   Return:
*   success or error
*
*   last update: Feb 19, 2010
*------------------------------------------------------------
*/

int gsl_rand_vector(gsl_rng *r,gsl_vector *v)
{
    if( v == NULL )
    	return GSL_EINVAL;

    unsigned int i;

    for( i = 0; i < v->size; i++ )
	v->data[i] = gsl_rng_uniform(r) + 1.0;

    return GSL_SUCCESS;
}
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   gsl_rand_matrix
*
*   This function creates a radom matrix.
*
*   Input:
*   seed - random number generator seed
*
*   Output:
*   M    - GSL matrix
*
*   Return:
*   success or error
*
*   last update: Feb 27, 2009
*------------------------------------------------------------
*/

int gsl_rand_matrix(gsl_rng *r,gsl_matrix *M)
{
    if( M == NULL )
    	return GSL_EINVAL;

    unsigned int i, j;

    for( i = 0; i < M->size1; i++ )
        for( j = 0; j < M->size2; j++ )
            M->data[i*M->tda+j] = gsl_rng_uniform(r) + 1.0;

    return GSL_SUCCESS;
}
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   gsl_diagonal_matrix
*
*   This function creates a diagonal matrix with the elements
*   of the vector v.
*
*   Input:
*   v - GSL vector
*
*   Output:
*   D  - GSL matrix
*
*   Return:
*   success or error
*
*   last update: Feb 26, 2009
*------------------------------------------------------------
*/

int gsl_diagonal_matrix(gsl_vector *v,gsl_matrix *D)
{
    if( v == NULL || D == NULL )
    	return GSL_EINVAL;

    unsigned int i, j;

    for( i = 0; i < D->size1; i++ )
        for( j = 0; j < D->size2; j++ )
            if( i == j )
                D->data[i*D->tda+j] = v->data[i];
            else
                D->data[i*D->tda+j] = 0.0;

    return GSL_SUCCESS;
}
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   util_time_used
*
*   This function computes the CPU time and Wall time
*   used in a simulation.
*
*   Input:
*   cpu_start  - cpu clock start flag
*   cpu_end    - cpu clock end flag
*   wall_start - wall clock start flag
*   wall_end   - wall clock end flag
*
*   last update: Nov 10, 2009
*------------------------------------------------------------
*/

void util_time_used(clock_t cpu_start,clock_t cpu_end,
			    time_t wall_start,time_t wall_end)
{
    double cpu_used;
    double wall_used;

    cpu_used  = ((double) (cpu_end - cpu_start)) / CLOCKS_PER_SEC;
    wall_used = difftime(wall_end,wall_start);

    printf("\n");
    printf("\n CPU  time (s): %+.6e",cpu_used);
    printf("\n Wall time (s): %+.6e",wall_used);
    printf("\n");

    return;
}
/*------------------------------------------------------------*/

max_stars_repo_name: americocunhajr/CRFlowLib

max_issues_repo_name: americocunhajr/CRFlowLib

max_forks_repo_name: americocunhajr/CRFlowLib


text:
#ifndef H_GPM_STRIPPED_DOWN
#define H_GPM_STRIPPED_DOWN

#include <gsl/gsl_vector.h>
#include <gsl/gsl_histogram.h>
#include "../csm_all.h"

void ght_find_theta_range(LDP laser_ref, LDP laser_sens,
                const double*x0, double max_linear_correction,
        double max_angular_correction_deg, int interval, gsl_histogram*hist, int*num_correspondences);

void ght_one_shot(LDP laser_ref, LDP laser_sens,
                const double*x0, double max_linear_correction,
                double max_angular_correction_deg, int interval, double*x, int*num_correspondences) ;

#endif


text:
#pragma once

#include <gsl\gsl>
#include <winrt\Windows.Foundation.h>
#include <d3d11.h>
#include "DrawableGameComponent.h"
#include "MatrixHelper.h"
#include "DirectionalLight.h"

namespace Library
{
	class Texture2D;
	class ProxyModel;
}

namespace Rendering
{
	class NormalMappingMaterial;

	class NormalMappingDemo final : public Library::DrawableGameComponent
	{
	public:
		NormalMappingDemo(Library::Game& game, const std::shared_ptr<Library::Camera>& camera);
		NormalMappingDemo(const NormalMappingDemo&) = delete;
		NormalMappingDemo(NormalMappingDemo&&) = default;
		NormalMappingDemo& operator=(const NormalMappingDemo&) = default;
		NormalMappingDemo& operator=(NormalMappingDemo&&) = default;
		~NormalMappingDemo();

		bool RealNormalMapEnabled() const;
		void SetRealNormalMapEnabled(bool enabled);
		void ToggleRealNormalMap();

		float AmbientLightIntensity() const;
		void SetAmbientLightIntensity(float intensity);

		float DirectionalLightIntensity() const;
		void SetDirectionalLightIntensity(float intensity);

		const DirectX::XMFLOAT3& LightDirection() const;
		void RotateDirectionalLight(DirectX::XMFLOAT2 amount);

		virtual void Initialize() override;
		virtual void Update(const Library::GameTime& gameTime) override;
		virtual void Draw(const Library::GameTime& gameTime) override;

	private:
		std::shared_ptr<NormalMappingMaterial> mMaterial;
		DirectX::XMFLOAT4X4 mWorldMatrix{ Library::MatrixHelper::Identity };
		winrt::com_ptr<ID3D11Buffer> mVertexBuffer;
		std::uint32_t mVertexCount{ 0 };
		Library::DirectionalLight mDirectionalLight;
		std::unique_ptr<Library::ProxyModel> mProxyModel;
		std::shared_ptr<Library::Texture2D> mRealNormalMap;
		std::shared_ptr<Library::Texture2D> mDefaultNormalMap;
		bool mUpdateMaterial{ true };
		bool mRealNormalMapEnabled{ true };
	};
}


text:
#ifndef OPENMC_TALLIES_FILTER_SPH_HARM_H
#define OPENMC_TALLIES_FILTER_SPH_HARM_H

#include <string>

#include <gsl/gsl>

#include "openmc/tallies/filter.h"

namespace openmc {

enum class SphericalHarmonicsCosine {
  scatter, particle
};

//==============================================================================
//! Gives spherical harmonics expansion moments of a tally score
//==============================================================================

class SphericalHarmonicsFilter : public Filter
{
public:
  //----------------------------------------------------------------------------
  // Constructors, destructors

  ~SphericalHarmonicsFilter() = default;

  //----------------------------------------------------------------------------
  // Methods

  std::string type() const override {return "sphericalharmonics";}

  void from_xml(pugi::xml_node node) override;

  void get_all_bins(const Particle* p, TallyEstimator estimator, FilterMatch& match)
  const override;

  void to_statepoint(hid_t filter_group) const override;

  std::string text_label(int bin) const override;

  //----------------------------------------------------------------------------
  // Accessors

  int order() const { return order_; }

  void set_order(int order);

  SphericalHarmonicsCosine cosine() const { return cosine_; }

  void set_cosine(gsl::cstring_span cosine);

private:
  //----------------------------------------------------------------------------
  // Data members

  int order_;

  //! The type of angle that this filter measures when binning events.
  SphericalHarmonicsCosine cosine_ {SphericalHarmonicsCosine::particle};
};

} // namespace openmc
#endif // OPENMC_TALLIES_FILTER_SPH_HARM_H


text:
#pragma once
#include <gsl/gsl>
#include <map>
#include <optional>

#include "aseprite_file.h"
#include "halley/text/halleystring.h"
#include "halley/maths/vector4.h"

namespace Halley
{
	class Path;
	class Image;
	struct ImageData;

	class AsepriteReader
	{
	public:
		static std::map<String, Vector<ImageData>> importAseprite(const String& spriteName, const Path& filename, gsl::span<const gsl::byte> fileData, bool trim, int padding, bool groupSeparated, bool sequenceSeparated);
		static void addImageData(int tagFrameNumber, int origFrameNumber, Vector<ImageData>& frameData, std::unique_ptr<Image> frameImage,
		                  const AsepriteFile& aseFile, const String& baseName, const String& sequence, const String& direction,
						  int duration, bool trim, int padding, bool hasFrameNumber, std::optional<String> group, bool firstImage,
		                  const String& spriteName, const Path& filename);
	};
}

max_stars_repo_name: amrezzd/halley

max_issues_repo_name: amrezzd/halley

max_forks_repo_name: amrezzd/halley


text:
/*
 * Copyright 2017 RIS.
 * Authored by Ron Dahlgren <ron@red83.net>
 */
#ifndef GLOBALS_H_
#define GLOBALS_H_

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <GL/glew.h>
#include <GLFW/glfw3.h>

#include <unistd.h>
#include <ctype.h>
#include <math.h>
#include <getopt.h>
#include <unistd.h>
#include <stdio.h>
#include <stdint.h>
#include <string.h>
#include <stdlib.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>

typedef int sd_result;

static const sd_result SD_OK = 0;
static const sd_result SD_RENDER_INIT_FAILURE = 1;
static const sd_result SD_WINDOW_CREATION_FAILURE = 2;
static const sd_result SD_NO_FREE_HANDLES = 3;
static const sd_result SD_CHECK_ERRNO = 4;
static const sd_result SD_GEN_ERROR = 5;
static const sd_result SD_SCRIPT_SYNTAX = 6;
static const sd_result SD_SCRIPT_TYPE_ERROR = 7;
static const sd_result SD_SCRIPT_GENERAL_ERROR = 8;
static const sd_result SD_FAILED_TO_READ_NUMBER = 9;

// Transform component indication flags
static const uint8_t XFORM_X_AXIS =         (1 << 0);
static const uint8_t XFORM_Y_AXIS =         (1 << 1);
static const uint8_t XFORM_Z_AXIS =         (1 << 2);
static const uint8_t XFORM_TRANSLATION =    (1 << 3);
static const uint8_t XFORM_TRANSLATION_X =  (1 << 4);
static const uint8_t XFORM_TRANSLATION_Y =  (1 << 5);
static const uint8_t XFORM_TRANSLATION_Z =  (1 << 6);

#endif


text:
/*
 * splinepdf.c: Routines to efficiently sample from a dimension of an
 *    N-dimensional spline surface. Uses a Metropolis-Hastings Markov Chain
 *    sampler.
 */

#include <math.h>
#include <gsl/gsl_rng.h>

#include "photospline/splinetable.h"
#include "photospline/bspline.h"
#include "photospline/splinepdf.h"

void logsplinepdf_n_sample(double *result, int results, int burnin,
    double *coords, int dim, struct splinetable *table, int derivatives,
    double (* proposal)(void*), double (* proposal_pdf)(double, double, void*),
    void* proposal_info, const gsl_rng *rng)
{
	int i;
	int centers[table->ndim];
	double val, lastval;
	double lastlogpdf = -INFINITY, logpdf;
	double lastproppdf, proppdf;
	double odds;
	double mint, maxt;
#ifdef DEBUG
	int accepted=0;
#endif

	/* Find the boundaries by choosing the first and last points
	 * with full support */
	mint = table->extents[dim][0];
	maxt = table->extents[dim][1];

	do {
		/*
		 * Get a starting point from the proposal distribution,
		 * making sure the PDF at the starting point is finite
		 * to avoid numerical problems.
		 */

		coords[dim] = lastval = (*proposal)(proposal_info);

		lastproppdf = (*proposal_pdf)(lastval,lastval,proposal_info);
		if (tablesearchcenters(table, coords, centers) != 0) {
			lastval = mint-1;
			continue;
		}
		lastlogpdf = ndsplineeval(table, coords, centers, 0);
		if (derivatives)
			lastlogpdf += log(ndsplineeval(table, coords, centers,
			    derivatives));
	} while (isnan(lastlogpdf) || lastval < mint || lastval > maxt);

	for (i = -burnin; i/burnin < results; i++) {
		coords[dim] = val = (*proposal)(proposal_info);

		/*
		 * If we ended up outside the table, reject the sample
		 * and try again
		 */
		if (val > maxt || val < mint ||
		    tablesearchcenters(table, coords, centers) != 0)
			goto reject;

		logpdf = ndsplineeval(table, coords, centers, 0);
		if (derivatives)
			logpdf += log(ndsplineeval(table, coords, centers,
			    derivatives));
		proppdf = (*proposal_pdf)(val,lastval,proposal_info);
		odds = exp(logpdf - lastlogpdf);
		odds *= lastproppdf/proppdf;

		if (odds > 1. || gsl_rng_uniform(rng) < odds) {
			/* Accept this value */
			lastval = val;
			lastlogpdf = logpdf;
			lastproppdf = proppdf;
		#ifdef DEBUG
			accepted++;
		#endif
		}

		/*
		 * Lastval has whatever we decided on now. If the burn-in
		 * period has elapsed, write it to output array.
		 *
		 * NB: This code is used in both the accept and reject
		 * cases.
		 */

	     reject:
		if (i >= 0 && (i % burnin) == 0)
			result[i/burnin] = lastval;
	}
#ifdef DEBUG
	printf("Efficiency: %e\n", (double)(accepted)/(double)(i));
#endif
}

void splinepdf_n_sample(double *result, int results, int burnin,
    double *coords, int dim, struct splinetable *table, int derivatives,
    double (* proposal)(void*), double (* proposal_pdf)(double, double, void*),
    void *proposal_info, const gsl_rng *rng)
{
	int i;
	int centers[table->ndim];
	double val, lastval;
	double lastpdf, pdf;
	double lastproppdf, proppdf;
	double odds;
	double mint, maxt;
#ifdef DEBUG
	int accepted=0;
#endif

	/* Find the boundaries by choosing the first and last points
	 * with full support */
	mint = table->extents[dim][0];
	maxt = table->extents[dim][1];

	/* Get a fully-supported starting point from the proposal distribution. */
	do {
		coords[dim] = lastval = (*proposal)(proposal_info);
	} while (lastval < mint || lastval > maxt ||
	    tablesearchcenters(table, coords, centers) != 0);

	lastproppdf = (*proposal_pdf)(lastval,lastval,proposal_info);
	lastpdf = ndsplineeval(table, coords, centers, derivatives);

	for (i = -burnin; i/burnin < results; i++) {
		coords[dim] = val = (*proposal)(proposal_info);

		/*
		 * If we ended up outside the table, reject the sample
		 * and try again
		 */
		if (val > maxt || val < mint ||
		    tablesearchcenters(table, coords, centers) != 0)
			goto reject;

		pdf = ndsplineeval(table, coords, centers, derivatives);
		proppdf = (*proposal_pdf)(val,lastval,proposal_info);
		odds = pdf/lastpdf;
		odds *= lastproppdf/proppdf;

		if (odds > 1. || gsl_rng_uniform(rng) < odds) {
			/* Accept this value */
			lastval = val;
			lastpdf = pdf;
			lastproppdf = proppdf;
			#ifdef DEBUG
				accepted++;
			#endif
		}

		/*
		 * Lastval has whatever we decided on now. If the burn-in
		 * period has elapsed, write it to output array.
		 *
		 * NB: This code is used in both the accept and reject
		 * cases.
		 */

	     reject:
		if (i >= 0 && (i % burnin) == 0)
			result[i / burnin] = lastval;
	}
	#ifdef DEBUG
		printf("Efficiency: %e\n", (double)(accepted)/(double)(i));
	#endif
}


text:
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multiroots.h>
#include "ray.h"

#define pi M_PI

//#############################################################################
//#############################################################################
//Mode finder

//Data structure for julia cavity data and mode parameters
typedef struct {
    //Julia cavity info and functions -----------------------------------------
    //pointer to Boundary object
    void *bnd;
    //pointer to radius function
    double (*rfunc_p)(void *bnd,double theta);
    //pointer to (mutating) radius and normal vector angle function
    void (*rsys_p)(void *bnd, double theta, double results[]);
    //pointer to RefractiveIndex object
    void *idx;
    //pointer to refractive index value function
    double (*nfunc_p)(void *idx, double r, double theta);
    //pointer to (mutating) refractive index value and derivative function
    void (*nderiv_p)(void *idx, double r, double theta, double results[]);

    //mode-specific parameters ------------------------------------------------
    long order; //number of bounces in a period
    double rtol; //stopping criterion tolerance

    //preallocated arrays -----------------------------------------------------
    double *raypath_r;
    double *raypath_theta;
    long *bounceindices;
    double *bouncepts_chi;
    long *lengths;
    double *modebounces;
} modeinfo;

//#############################################################################
//Function for rootfinding
int rootfunc(const gsl_vector *x, void *params, gsl_vector *f){

    //params is pointer to modeinfo struct containing both cavity info and
    //simulation parameters
    modeinfo *mip = (modeinfo *)params;
    const long order = (*mip).order;
    double *raypath_r = (double *)(*mip).raypath_r;
    double *raypath_theta = (double *)(*mip).raypath_theta;
    long *bounceindices = (long *)(*mip).bounceindices;
    double *bouncepts_chi = (double *)(*mip).bouncepts_chi;
    long *lengths = (long *)(*mip).lengths;

    //Extract initial condition data (print for debugging)
    const double theta0 = gsl_vector_get(x,0), chi0 = gsl_vector_get(x,1);
    //printf("theta0 = %.8f, chi0 = %.8f\n",theta0,chi0);

    //Check for sensible values
    if(fabs(chi0) > pi/2){
        //Out of domain
        return GSL_EDOM;
    }

    //Setup initial ODE coordinate vector components
    double results[2];
    (*(*mip).rsys_p)((*mip).bnd,theta0,results);
    //result[0] = r0, results[1] = normang0
    const double n = (*(*mip).nfunc_p)((*mip).idx,results[0],theta0);
    const double phi0 = results[1] + chi0;
    const double pr0 = n*cos(phi0-theta0);
    const double ptheta0 = n*results[0]*sin(phi0-theta0);

    //Compute ray
    rayevolve(
        //Storage arrays
        raypath_r,raypath_theta,bounceindices,bouncepts_chi,lengths,
        //Initial conditions
        results[0],theta0,pr0,ptheta0,
        //Simulation parameters
        200.0,order+1,1e-12,1e-12,
        //Cavity parameters
        (*mip).bnd,(*mip).rfunc_p,(*mip).rsys_p,
        (*mip).idx,(*mip).nfunc_p,(*mip).nderiv_p);

    //Store difference in initial and final bounce locations (print for debugging)
    const double dtheta = fmod(raypath_theta[bounceindices[order]-1] - raypath_theta[bounceindices[0]-1] + pi, 2*pi) - pi;
    const double dchi = bouncepts_chi[order] - bouncepts_chi[0];
    gsl_vector_set(f,0,dtheta);
    gsl_vector_set(f,1,dchi);
    //printf("dtheta = %.8f, dchi = %.8f\n",dtheta,dchi);

    //Test stopping criterion
    if(gsl_multiroot_test_residual(f,(*mip).rtol) != GSL_CONTINUE){
        //Record mode bounce data in results array (to avoid running simulation again)
        double *modebounces = (double *)(*mip).modebounces;
        //theta values stored in first half of array (first column after reshaping)
        //chi values stored in second half of array (second column after reshaping)
        int i;
        for(i=0; i<order; i+=1){
            modebounces[i] = raypath_theta[bounceindices[i]-1];
            modebounces[order+i] = bouncepts_chi[i];
        }
    }

    return GSL_SUCCESS;
}

//#############################################################################
//Solver function
int findmode(
    //mode-specific variables
    long order, double modebounces[],
    //initial position on PSS
    double theta0, double chi0,
    //root-finding algorithm parameters
    double rtol, long maxiter,
    //Cavity parameters
    void *bnd, double (*rfunc_p)(void *bnd,double theta),
    void (*rsys_p)(void *bnd, double theta,double results[]),
    void *idx, double (*nfunc_p)(void *idx, double r, double theta),
    void (*nderiv_p)(void *idx, double r, double theta,double results[])){

        //Preallocate arrays for reuse
        long bounceindices[30], lengths[2];
        double raypath_r[50000],raypath_theta[50000],bouncepts_chi[30];

        //Condense input mode data into modeinfo structs
        modeinfo mi = {bnd,rfunc_p,rsys_p,idx,nfunc_p,nderiv_p,order,rtol,
            raypath_r,raypath_theta,bounceindices,bouncepts_chi,lengths,modebounces};

        //Store initial guess
        gsl_vector *x0 = gsl_vector_alloc(2);
        gsl_vector_set(x0,0,theta0);
        gsl_vector_set(x0,1,chi0);

        //Initiate solver
        const gsl_multiroot_fsolver_type *fs_t = gsl_multiroot_fsolver_dnewton;
        gsl_multiroot_fsolver *fs = gsl_multiroot_fsolver_alloc(fs_t,2);
        gsl_multiroot_function mrf = {rootfunc,2,&mi};
        gsl_multiroot_fsolver_set(fs,&mrf,x0);
        gsl_set_error_handler_off();

        int status, count = 0;

        //Iterate solver
        while(count != maxiter){
            status = gsl_multiroot_fsolver_iterate(fs);
            if(status != GSL_SUCCESS) break; //Some error has occurred

            gsl_vector *f = gsl_multiroot_fsolver_f(fs);
            if(gsl_multiroot_test_residual(f,rtol) != GSL_CONTINUE) break; //Done!

            count += 1;
        }

        //Free memory
        gsl_multiroot_fsolver_free(fs);
        gsl_vector_free(x0);

        //In case of error or aborted computation, indicate no result
        if((status != GSL_SUCCESS) || (count == maxiter)) return 1;

        return 0;
}