text:
#include <stdlib.h>
#include <math.h>
#include <complex.h>
#include <gsl/gsl_sf_result.h>
#include <gsl/gsl_sf_gamma.h>
#include <fftw3.h>
//#include "fftlog.h"

#define PI 3.14159265358979323846
#define PI_2 1.5707963267948966 /// pi/2
#define PI_times_2 6.28318530718 // 2*pi
#define inv_PI 0.31830988618 // 1/pi
#define log2 0.6931471805599453 //ln(2)
#define PI2cubed 248.05021344239853 //(2*pi)^3

/* Compute the discrete Hankel transform of the function a(r).  See the FFTLog
 * documentation (or the Fortran routine of the same name in the FFTLog
 * sources) for a description of exactly what this function computes.
 * If u is NULL, the transform coefficients will be computed anew and discarded
 * afterwards.  If you plan on performing many consecutive transforms, it is
 * more efficient to pre-compute the u coefficients. */
void fht(int N,  const double r[],  const double complex a[], double k[],
	 double complex b[], double mu,
         double q, double kcrc, int noring, double complex* u);
//         double q = 0, double kcrc = 1, bool noring = true,
//         double complex* u = NULL);

/* Pre-compute the coefficients that appear in the FFTLog implementation of
 * the discrete Hankel transform.  The parameters N, mu, and q here are the
 * same as for the function fht().  The parameter L is defined (for whatever
 * reason) to be N times the logarithmic spacing of the input array, i.e.
 *   L = N * log(r[N-1]/r[0])/(N-1) */
void compute_u_coefficients(int N, double mu, double q, double L, double kcrc,
			    double complex u[]);

/* This code is FFTLog, which is described in arXiv:astro-ph/9905191 */

static double complex lngamma_fftlog(double complex z)
{
  gsl_sf_result lnr, phi;
  gsl_sf_lngamma_complex_e(creal(z), cimag(z), &lnr, &phi);
  return lnr.val + I*phi.val;
}

/* Unused function -- commenting out to prevent clang warnings
static double complex gamma_fftlog(double complex z)
{
  return cexp(lngamma_fftlog(z));
  }
*/

static double complex polar (double r, double phi)
{
  return (r*cos(phi) +I*(r*sin(phi)));
}

static void lngamma_4(double x, double y, double* lnr, double* arg)
{
  double complex w = lngamma_fftlog(x+y*I);
  if(lnr) *lnr = creal(w);
  if(arg) *arg = cimag(w);
}

static double goodkr(int N, double mu, double q, double L, double kr)
{
  double N_L = N/L;
  double xp = (mu+1+q)/2;
  double xm = (mu+1-q)/2;
  double y = PI_2*N_L;
  double lnr, argm, argp;
  lngamma_4(xp, y, &lnr, &argp);
  lngamma_4(xm, y, &lnr, &argm);
  double arg = log(2./kr) * N_L + (argp + argm)*inv_PI;
  double iarg = round(arg);
  if(arg != iarg)
    kr *= exp((arg - iarg)/N_L);
  return kr;
}

void compute_u_coefficients(int N, double mu, double q, double L,
			    double kcrc, double complex u[])
{
  double y = PI/L;
  double k0r0 = kcrc * exp(-L);
  double t = -2*y*log(k0r0*0.5);

  if(q == 0) {
    double x = (mu+1)*0.5;
    double lnr, phi;
    for(int m = 0; m <= N/2; m++) {
      lngamma_4(x, m*y, &lnr, &phi);
      u[m] = polar(1.0,m*t + 2*phi);
    }
  }
  else {
    double xp = (mu+1+q)*0.5;
    double xm = (mu+1-q)*0.5;
    double lnrp, phip, lnrm, phim;
    for(int m = 0; m <= N/2; m++) {
      lngamma_4(xp, m*y, &lnrp, &phip);
      lngamma_4(xm, m*y, &lnrm, &phim);
      u[m] = polar(exp(q*log2 + lnrp - lnrm), m*t + phip - phim);
    }
  }

  for(int m = N/2+1; m < N; m++)
    u[m] = conj(u[N-m]);
  if((N % 2) == 0)
    u[N/2] = (creal(u[N/2]) + I*0.0);
}

void fht(int N, const double r[], const double complex a[], double k[],
	 double complex b[], double mu, double q, double kcrc,
	 int noring, double complex* u)
{
  double L = log(r[N-1]/r[0]) * N/(N-1.);
  double complex* ulocal = NULL;
  if(u == NULL) {
    if(noring)
      kcrc = goodkr(N, mu, q, L, kcrc);
    ulocal = malloc(sizeof(complex double)*N);
    compute_u_coefficients(N, mu, q, L, kcrc, ulocal);
    u = ulocal;
  }

  /* Compute the convolution b = a*u using FFTs */
  fftw_plan forward_plan = fftw_plan_dft_1d(N, (fftw_complex*) a,
					    (fftw_complex*) b,  -1,
					    FFTW_ESTIMATE);
  fftw_plan reverse_plan = fftw_plan_dft_1d(N, (fftw_complex*) b,
					    (fftw_complex*) b, +1,
					    FFTW_ESTIMATE);
  fftw_execute(forward_plan);
  for(int m = 0; m < N; m++)
    b[m] *= u[m] / (double)(N);       // divide by N since FFTW doesn't normalize the inverse FFT
  fftw_execute(reverse_plan);
  fftw_destroy_plan(forward_plan);
  fftw_destroy_plan(reverse_plan);

  /* Reverse b array */
  double complex tmp;
  for(int n = 0; n < N/2; n++) {
    tmp = b[n];
    b[n] = b[N-n-1];
    b[N-n-1] = tmp;
  }

  /* Compute k's corresponding to input r's */
  double k0r0 = kcrc * exp(-L);
  k[0] = k0r0/r[0];
  for(int n = 1; n < N; n++)
    k[n] = k[0] * exp(n*L/N);

  free(ulocal);
}

void fftlog_ComputeXi2D(double bessel_order,int N,const double l[],
			const double cl[], double th[], double xi[])
{
  double complex* a = malloc(sizeof(complex double)*N);
  double complex* b = malloc(sizeof(complex double)*N);

  for(int i=0;i<N;i++)
    a[i]=l[i]*cl[i];
  fht(N,l,a,th,b,bessel_order,0,1,1,NULL);
  for(int i=0;i<N;i++)
    xi[i]=creal(b[i]/(2*PI*th[i]));

  free(a);
  free(b);
}

void fftlog_ComputeXiLM(double l, double m, int N, const double k[],
			const double pk[], double r[], double xi[],
			double q)
{
  static int init_flag = 0;
  static int N_old = -99999;
  static double complex* a = NULL;
  static double complex* b = NULL;
  if ((init_flag == 0) || (N != N_old)){
    if (a!=NULL) free(a);
    if (b!=NULL) free(b);
    a = malloc(sizeof(complex double)*N);
    b = malloc(sizeof(complex double)*N);
    N_old = N;
    init_flag = 1;
  }

  for(int i = 0; i < N; i++)
    a[i] = pow(k[i], m-0.5) * pk[i];
  fht(N, k, a, r, b, l + 0.5, q, 1, 1, NULL);
  for(int i = 0; i < N; i++)
    xi[i] = creal(pow(PI_times_2*r[i], 0.5-m) * b[i]);

}

void pk2xi(int N, const double k[], const double pk[], double r[], double xi[])
{
  fftlog_ComputeXiLM(0, 2, N, k, pk, r, xi, 0);
}

void xi2pk(int N, const double r[], const double xi[], double k[], double pk[])
{
  //static const double TwoPiCubed = 8*PI*PI*PI;
  fftlog_ComputeXiLM(0, 2, N, r, xi, k, pk, 0);
  for(int j = 0; j < N; j++)
    pk[j] *= PI2cubed;
}

max_forks_repo_name: tmcclintock/FastCorr


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

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_multiset.h>

gsl_multiset *
gsl_multiset_alloc (const size_t n, const size_t k)
{
  gsl_multiset * c;

  if (n == 0)
    {
      GSL_ERROR_VAL ("multiset parameter n must be positive integer",
                        GSL_EDOM, 0);
    }
  if (k > n)
    {
      GSL_ERROR_VAL ("multiset length k must be an integer less than or equal to n",
                        GSL_EDOM, 0);
    }
  c = (gsl_multiset *) malloc (sizeof (gsl_multiset));

  if (c == 0)
    {
      GSL_ERROR_VAL ("failed to allocate space for multiset struct",
                        GSL_ENOMEM, 0);
    }

  if (k > 0)
    {
      c->data = (size_t *) malloc (k * sizeof (size_t));

      if (c->data == 0)
        {
          free (c);             /* exception in constructor, avoid memory leak */

          GSL_ERROR_VAL ("failed to allocate space for multiset data",
                         GSL_ENOMEM, 0);
        }
    }
  else
    {
      c->data = 0;
    }

  c->n = n;
  c->k = k;

  return c;
}

gsl_multiset *
gsl_multiset_calloc (const size_t n, const size_t k)
{
  size_t i;

  gsl_multiset * c =  gsl_multiset_alloc (n, k);

  if (c == 0)
    return 0;

  /* initialize multiset to repeated first element */

  for (i = 0; i < k; i++)
    {
      c->data[i] = 0;
    }

  return c;
}

void
gsl_multiset_init_first (gsl_multiset * c)
{
  const size_t k = c->k ;
  size_t i;

  /* initialize multiset to repeated first element */

  for (i = 0; i < k; i++)
    {
      c->data[i] = 0;
    }
}

void
gsl_multiset_init_last (gsl_multiset * c)
{
  const size_t k = c->k ;
  size_t i;
  size_t n = c->n;

  /* initialize multiset to repeated last element */

  for (i = 0; i < k; i++)
    {
      c->data[i] = n - 1;
    }
}

void
gsl_multiset_free (gsl_multiset * c)
{
  RETURN_IF_NULL (c);
  if (c->k > 0) free (c->data);
  free (c);
}

ext: c


text:
#pragma once

#include <winrt\Windows.Foundation.h>
#include <d3d11.h>
#include <gsl/gsl>
#include "RTTI.h"

namespace Library
{
	class Shader : public RTTI
	{
		RTTI_DECLARATIONS(Shader, RTTI)

	public:
		Shader(const Shader&) = default;
		Shader& operator=(const Shader&) = default;
		Shader(Shader&&) = default;
		Shader& operator=(Shader&&) = default;
		virtual ~Shader() = default;

		static winrt::com_ptr<ID3D11ClassLinkage> CreateClassLinkage(gsl::not_null<ID3D11Device*> device);

	protected:
		Shader() = default;
	};
}

ext: h


text:
#include <stdio.h>
#include <stdlib.h>
#include <gbpLib.h>
#include <gbpMisc.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_roots.h>

double bisect_array_function(double x_interp, void *params) {
    double interp;
    double r_val;
    interp = interpolate(((bisect_af_params *)params)->interp, x_interp);
    r_val  = interp - ((bisect_af_params *)params)->value;
    return (r_val);
}

double bisect_array(interp_info *interp, double value, double threshold) {
    double *                     x;
    double *                     y;
    double *                     d2y;
    double                       x_lo;
    double                       x_hi;
    double                       r_val;
    int                          i;
    int                          alloc_flag = GBP_FALSE;
    int                          iter       = 0;
    int                          max_iter   = 500;
    int                          status     = GSL_CONTINUE;
    bisect_af_params             params;
    gsl_function                 F;
    const gsl_root_fsolver_type *T;
    gsl_root_fsolver *           s;

    params.interp = interp;
    params.value  = value;
    F.function    = bisect_array_function;
    F.params      = (void *)(&params);

    x_lo = interp->x[0];
    x_hi = interp->x[interp->n - 1];

    if(fabs(bisect_array_function(x_lo, &params)) < threshold)
        return (x_lo);
    if(fabs(bisect_array_function(x_hi, &params)) < threshold)
        return (x_hi);

    T = gsl_root_fsolver_brent;
    s = gsl_root_fsolver_alloc(T);
    gsl_root_fsolver_set(s, &F, x_lo, x_hi);

    while(status == GSL_CONTINUE && iter < max_iter) {
        status = gsl_root_fsolver_iterate(s);
        r_val  = gsl_root_fsolver_root(s);
        x_lo   = gsl_root_fsolver_x_lower(s);
        x_hi   = gsl_root_fsolver_x_upper(s);
        status = gsl_root_test_interval(x_lo, x_hi, 0, threshold);
        iter++;
    }

    // Clean-up
    gsl_root_fsolver_free(s);

    return (r_val);
}

max_stars_repo_name: gbpoole/gbpCode

max_issues_repo_name: gbpoole/gbpCode

max_forks_repo_name: gbpoole/gbpCode


text:
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <nlopt.h>
#include <unistd.h>
#include <sys/time.h>

void quadbezier(const double *a, const double *b, const double *c,
		double t, double *out)
{
  double onemt = 1.0 - t;
  out[0] = onemt * (onemt * a[0] + t * b[0]) + t * (onemt * b[0] + t * c[0]);
  out[1] = onemt * (onemt * a[1] + t * b[1]) + t * (onemt * b[1] + t * c[1]);
}

void linearinterp(const double *a, const double *b,
		  double t, double *out)
{
  out[0] = (b[0] - a[0]) * t + a[0];
  out[1] = (b[1] - a[1]) * t + a[1];
}

typedef struct {
  int N;
  int steps;
  double *starts;
  double *stops;
  double *scratch;
} problem_t;

double randd()
{
  int n = rand() % (1<<24);
  return ((double)n) / ((double)(1<<24));
}

double distance(const double *a, const double *b)
{
  return sqrt((b[0] - a[0]) * (b[0] - a[0]) +
	      (b[1] - a[1]) * (b[1] - a[1]));
}

#define RADIUS 5.0

double errorf(unsigned n, const double *pts, double *grad,
	      void *data)
{
  problem_t *setup = (problem_t *)data;
  double error = 0.0;

  for(int step=0; step<setup->steps; step++) {
    double t = (1.0+step) / (1.0+setup->steps);
    for(int i=0; i<setup->N; i++)
      quadbezier(&setup->starts[2*i],
		 &pts[2*i],
		 &setup->stops[2*i],
		 t,
		 &setup->scratch[2*i]);

    for(int i=0; i<setup->N; i++)
      for(int j=i+1; j<setup->N; j++) {
	double d = distance(&setup->scratch[2*i],
			    &setup->scratch[2*j]);
	d /= RADIUS;
	error += (1.0/d) * (1.0/d);
      }
  }

  return error / setup->steps;
}

int main(int argc, char *argv[])
{
  int N = 12;

  double *starts = (double *)malloc(sizeof(double)*N*2);
  double *stops = (double *)malloc(sizeof(double)*N*2);

  srand(2);
  nlopt_srand(1);

  double scaling = 100.0;
  for(int i=0; i<N; i++) {
    starts[2*i+0] = scaling * randd();
    starts[2*i+1] = scaling * randd();
    stops[2*i+0]  = scaling * randd();
    stops[2*i+1]  = scaling * randd();
  }

  double hs = scaling / 2.0;
  for(int i=0; i<N; i++) {
    double theta = ((double)i)/N * 2 * 3.141592;
    starts[2*i+0] = hs * cos(theta) + hs;
    starts[2*i+1] = hs * sin(theta) + hs;
    stops[2*i+0]  = hs * cos(theta + 3.141592) + hs;
    stops[2*i+1]  = hs * sin(theta + 3.141592) + hs;
  }

  struct timeval start, stop;

  double *guesses = (double *)malloc(sizeof(double)*N*2);

  gettimeofday(&start, NULL);
  for(int i=0; i<N; i++) {
    linearinterp(&starts[2*i], &stops[2*i], 0.5, &guesses[2*i]);
  }

  nlopt_opt global = nlopt_create(NLOPT_GN_ISRES, 2*N);
  nlopt_opt local = nlopt_create(NLOPT_LN_SBPLX, 2*N);

  problem_t *setup = (problem_t *)malloc(sizeof(problem_t));
  setup->N = N;
  setup->starts = starts;
  setup->stops = stops;
  setup->steps = 50;
  setup->scratch = (double *)malloc(sizeof(double)*N*2);

  nlopt_result res;

  res = nlopt_set_min_objective(global, errorf, setup);
  if(res!=1) { exit(-1); }

  res = nlopt_set_lower_bounds1(global, -25.0);
  if(res!=1) { exit(-1); }

  res = nlopt_set_upper_bounds1(global, 125.0);
  if(res!=1) { exit(-1); }

  res = nlopt_set_stopval(global, 1.0);
  if(res!=1) { exit(-1); }

  res = nlopt_set_maxeval(global, 100);
  if(res!=1) { exit(-1); }

  double final_error;
  res = nlopt_optimize(global, guesses, &final_error);

  res = nlopt_set_min_objective(local, errorf, setup);
  if(res!=1) { exit(-1); }

  res = nlopt_set_lower_bounds1(local, -25.0);
  if(res!=1) { exit(-1); }

  res = nlopt_set_upper_bounds1(local, 125.0);
  if(res!=1) { exit(-1); }

  res = nlopt_set_stopval(local, 1.0);
  if(res!=1) { exit(-1); }

  res = nlopt_set_maxeval(local, 1000);
  if(res!=1) { exit(-1); }

  res = nlopt_optimize(local, guesses, &final_error);
  gettimeofday(&stop, NULL);

  fprintf(stderr, "final: %f\n", final_error);

  for(int i=0; i<N; i++) {
    printf("%f %f %f %f %f %f\n",
	   starts[2*i+0], starts[2*i+1],
	   guesses[2*i+0], guesses[2*i+1],
	   stops[2*i+0], stops[2*i+1]);
  }

  double startd, stopd;
  startd = start.tv_sec;
  startd += ((double)(start.tv_usec)) / 1e6;
  stopd = stop.tv_sec;
  stopd += ((double)(stop.tv_usec)) / 1e6;

  //printf("%f\n", stopd - startd);
}

max_stars_repo_name: jkominek/discs

max_issues_repo_name: jkominek/discs

max_forks_repo_name: jkominek/discs


text:
//popGen stuff for Next Gen Sequencing Data

#include "pgSummaryStats.h"
#include <math.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include "vector.h"
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_matrix.h>

//Sampling Probability stuff

//prob_i returns the probability under SNM of a snp freq i
double prob_i(int i, int n){
	double tmp;
	tmp = 1.0 / (double) i;
	return(tmp/ a1f(n));
}

double prob_km(int k, int m, int n){
	double sum = 0.0;
	int i;
	double tmpFreq;

	for(i=1;i<n;i++){
		tmpFreq = (double) i / (double) n;
		sum += gsl_ran_binomial_pdf(k,tmpFreq,m) * prob_i(i,n);
	}
	return(sum);
}

//achaz system stuff
double achazSFS(gsl_matrix *sfsCount, gsl_matrix *weightMat){
	int i, j;
	double theta = 0.0;
	double weightSum, tmpTheta;

	//rows = sampleSizes
	for(i = 2; i < sfsCount->size1; i++){
		weightSum = 0.0;
		tmpTheta = 0.0;
		for(j = 1; j < i; j++){
			//size i freq j
			tmpTheta += gsl_matrix_get(sfsCount, i, j) * j * gsl_matrix_get(weightMat,i,j);
			weightSum += gsl_matrix_get(weightMat,i,j);
		}
		if(weightSum > 0){ theta += tmpTheta / weightSum; }
	}
	return(theta);
}

double achazSFSCorrected(gsl_matrix *sfsCount, gsl_matrix *weightMat,int nPool){
	int i, j;
	double theta = 0.0;
	double weightSum, tmpTheta;

	//rows = sampleSizes
	for(i = 2; i < sfsCount->size1; i++){
		weightSum = 0.0;
		tmpTheta = 0.0;
		for(j = 1; j < i; j++){
			//size i freq j
			tmpTheta += gsl_matrix_get(sfsCount, i, j) * gsl_matrix_get(weightMat,i,j) * (1.0/prob_km(j,i,nPool));
			weightSum += gsl_matrix_get(weightMat,i,j);
		}
		if(weightSum > 0){

			theta += tmpTheta / weightSum / a1f(nPool);
		}
	}
	return(theta);
}

// fillPiWeights -- returns the weight matrix for theta_pi
void fillPiWeights(gsl_matrix *weightMat, int minFreq){
	int i, j;
	gsl_matrix_set_zero(weightMat);
	for(i=minFreq+1; i < weightMat->size1; i++){
		for(j=minFreq;j<i;j++){
			gsl_matrix_set(weightMat,i,j,i-j);
		}
	}
}

// fillHWeights -- returns the weight matrix for theta_H
void fillHWeights(gsl_matrix *weightMat, int minFreq){
	int i, j;
	gsl_matrix_set_zero(weightMat);
	for(i=minFreq+1; i < weightMat->size1; i++){
		for(j=minFreq;j<i;j++){
			gsl_matrix_set(weightMat,i,j,j);
		}
	}
}

// fillWWeights -- returns the weight matrix for theta_w
void fillWWeights(gsl_matrix *weightMat, int minFreq){
	int i, j;
	gsl_matrix_set_zero(weightMat);
	for(i=minFreq+1; i < weightMat->size1; i++){
		for(j=minFreq;j<i;j++){
			gsl_matrix_set(weightMat,i,j,1/(double)j);
		}
	}
}

max_stars_repo_name: andrewkern/pgLib

max_issues_repo_name: andrewkern/pgLib

max_forks_repo_name: andrewkern/pgLib


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

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

/*! \file predict.c
 *  \brief drift particles by a small time interval
 *
 *  This function contains code to implement a drift operation on all the
 *  particles, which represents one part of the leapfrog integration scheme.
 */

/*! This function drifts all particles from the current time to the future:
 *  time0 - > time1
 *
 *  If there is no explicit tree construction in the following timestep, the
 *  tree nodes are also drifted and updated accordingly. Note: For periodic
 *  boundary conditions, the mapping of coordinates onto the interval
 *  [0,All.BoxSize] is only done before the domain decomposition, or for
 *  outputs to snapshot files.  This simplifies dynamic tree updates, and
 *  allows the domain decomposition to be carried out only every once in a
 *  while.
 */
void move_particles(int time0, int time1)
{
  int i, j;
  double dt_drift, dt_gravkick, dt_hydrokick, dt_entr;
  double t0, t1;

  t0 = second();

  if(All.ComovingIntegrationOn)
    {
      dt_drift = get_drift_factor(time0, time1);
      dt_gravkick = get_gravkick_factor(time0, time1);
      dt_hydrokick = get_hydrokick_factor(time0, time1);
    }
  else
    {
      dt_drift = dt_gravkick = dt_hydrokick = (time1 - time0) * All.Timebase_interval;
    }

  for(i = 0; i < NumPart; i++)
    {
      for(j = 0; j < 3; j++)
	P[i].Pos[j] += P[i].Vel[j] * dt_drift;

      if(P[i].Type == 0)
	{
#ifdef PMGRID
	  for(j = 0; j < 3; j++)
	    SphP[i].VelPred[j] +=
	      (P[i].GravAccel[j] + P[i].GravPM[j]) * dt_gravkick + SphP[i].HydroAccel[j] * dt_hydrokick;
#else
	  for(j = 0; j < 3; j++)
	    SphP[i].VelPred[j] += P[i].GravAccel[j] * dt_gravkick + SphP[i].HydroAccel[j] * dt_hydrokick;
#endif
	  SphP[i].Density *= exp(-SphP[i].DivVel * dt_drift);
	  SphP[i].Hsml *= exp(0.333333333333 * SphP[i].DivVel * dt_drift);

	  if(SphP[i].Hsml < All.MinGasHsml)
	    SphP[i].Hsml = All.MinGasHsml;

	  dt_entr = (time1 - (P[i].Ti_begstep + P[i].Ti_endstep) / 2) * All.Timebase_interval;

	  SphP[i].Pressure = (SphP[i].Entropy + SphP[i].DtEntropy * dt_entr) * pow(SphP[i].Density, GAMMA);

#ifdef MORRIS97VISC
          SphP[i].Alpha += SphP[i].DAlphaDt * dt_drift;
#endif
	}
    }

  /* if domain-decomp and tree are not going to be reconstructed, update dynamically.  */
  if(All.NumForcesSinceLastDomainDecomp < All.TotNumPart * All.TreeDomainUpdateFrequency)
    {
      for(i = 0; i < Numnodestree; i++)
	for(j = 0; j < 3; j++)
	  Nodes[All.MaxPart + i].u.d.s[j] += Extnodes[All.MaxPart + i].vs[j] * dt_drift;

      force_update_len();

      force_update_pseudoparticles();
    }

  t1 = second();

  All.CPU_Predict += timediff(t0, t1);
}

/*! This function makes sure that all particle coordinates (Pos) are
 *  periodically mapped onto the interval [0, BoxSize].  After this function
 *  has been called, a new domain decomposition should be done, which will
 *  also force a new tree construction.
 */
#ifdef PERIODIC
void do_box_wrapping(void)
{
  int i, j;
  double boxsize[3];

  for(j = 0; j < 3; j++)
    boxsize[j] = All.BoxSize;

#ifdef LONG_X
  boxsize[0] *= LONG_X;
#endif
#ifdef LONG_Y
  boxsize[1] *= LONG_Y;
#endif
#ifdef LONG_Z
  boxsize[2] *= LONG_Z;
#endif

  for(i = 0; i < NumPart; i++)
    for(j = 0; j < 3; j++)
      {
	while(P[i].Pos[j] < 0)
	  P[i].Pos[j] += boxsize[j];

	while(P[i].Pos[j] >= boxsize[j])
	  P[i].Pos[j] -= boxsize[j];
      }
}
#endif

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

int
main (void)
{
  int i;
  double x = 0, y = 0, dx, dy;

  const gsl_rng_type * T;
  gsl_rng * r;

  gsl_rng_env_setup();
  T = gsl_rng_default;
  r = gsl_rng_alloc (T);

  printf ("%g %g\n", x, y);

  for (i = 0; i < 10; i++)
    {
      gsl_ran_dir_2d (r, &dx, &dy);
      x += dx; y += dy;
      printf ("%g %g\n", x, y);
    }

  gsl_rng_free (r);
  return 0;
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


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

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

#include "ccl.h"

#ifdef HAVE_ANGPOW
#include "Angpow/angpow_ccl.h"
#endif

#define CCL_FRAC_RELEVANT 5E-4
//#define CCL_FRAC_RELEVANT 1E-3
//Gets the x-interval where the values of y are relevant
//(meaning, that the values of y for those x are at least above a fraction frac of its maximum)
static void get_support_interval(int n,double *x,double *y,double frac,
				 double *xmin_out,double *xmax_out)
{
  int ix;
  double ythr=-1000;

  //Initialize as the original edges in case we don't find an interval
  *xmin_out=x[0];
  *xmax_out=x[n-1];

  //Find threshold
  for(ix=0;ix<n;ix++) {
    if(y[ix]>ythr) ythr=y[ix];
  }
  ythr*=frac;

  //Find minimum
  for(ix=0;ix<n;ix++) {
    if(y[ix]>=ythr) {
      *xmin_out=x[ix];
      break;
    }
  }

  //Find maximum
  for(ix=n-1;ix>=0;ix--) {
    if(y[ix]>=ythr) {
      *xmax_out=x[ix];
      break;
    }
  }
}

//Wrapper around spline_eval with GSL function syntax
static double speval_bis(double x,void *params)
{
  return ccl_spline_eval(x,(SplPar *)params);
}

void ccl_cl_workspace_free(CCL_ClWorkspace *w)
{
  free(w->l_arr);
  free(w);
}

CCL_ClWorkspace *ccl_cl_workspace_new(int lmax,int l_limber,
				      double l_logstep,int l_linstep,int *status)
{
  int i_l,l0,increment;
  CCL_ClWorkspace *w=(CCL_ClWorkspace *)malloc(sizeof(CCL_ClWorkspace));
  if(w==NULL)
    *status=CCL_ERROR_MEMORY;

  if(*status==0) {
    //Set params
    w->lmax=lmax;
    w->l_limber=l_limber;
    w->l_logstep=l_logstep;
    w->l_linstep=l_linstep;

    //Compute number of multipoles
    i_l=0; l0=0;
    increment=CCL_MAX(((int)(l0*(w->l_logstep-1.))),1);
    while((l0 < w->lmax) && (increment < w->l_linstep)) {
      i_l++;
      l0+=increment;
      increment=CCL_MAX(((int)(l0*(w->l_logstep-1))),1);
    }
    increment=w->l_linstep;
    while(l0 < w->lmax) {
      i_l++;
      l0+=increment;
    }

    //Allocate array of multipoles
    w->n_ls=i_l+1;
    w->l_arr=(int *)malloc(w->n_ls*sizeof(int));
    if(w->l_arr==NULL)
      *status=CCL_ERROR_MEMORY;
  }

  if(*status==0) {
    //Redo the computation above and store values of ell
    i_l=0; l0=0;
    increment=CCL_MAX(((int)(l0*(w->l_logstep-1.))),1);
    while((l0 < w->lmax) && (increment < w->l_linstep)) {
      w->l_arr[i_l]=l0;
      i_l++;
      l0+=increment;
      increment=CCL_MAX(((int)(l0*(w->l_logstep-1))),1);
    }
    increment=w->l_linstep;
    while(l0 < w->lmax) {
      w->l_arr[i_l]=l0;
      i_l++;
      l0+=increment;
    }
    //Don't go further than lmaw
    w->l_arr[w->n_ls-1]=w->lmax;
  }

  return w;
}

CCL_ClWorkspace *ccl_cl_workspace_new_limber(int lmax,double l_logstep,int l_linstep,int *status)
{
  return ccl_cl_workspace_new(lmax,-1,l_logstep,l_linstep,status);
}

//Params for lensing kernel integrand
typedef struct {
  double chi;
  SplPar *spl_pz;
  ccl_cosmology *cosmo;
  int *status;
} IntLensPar;

//Integrand for lensing kernel
static double integrand_wl(double chip,void *params)
{
  IntLensPar *p=(IntLensPar *)params;
  double chi=p->chi;
  double a=ccl_scale_factor_of_chi(p->cosmo,chip, p->status);
  double z=1./a-1;
  double pz=ccl_spline_eval(z,p->spl_pz);
  double h=p->cosmo->params.h*ccl_h_over_h0(p->cosmo,a, p->status)/CLIGHT_HMPC;

  if(chi==0)
    return h*pz;
  else
    return h*pz*ccl_sinn(p->cosmo,chip-chi,p->status)/ccl_sinn(p->cosmo,chip,p->status);
}

//Integral to compute lensing window function
//chi     -> comoving distance
//cosmo   -> ccl_cosmology object
//spl_pz  -> normalized N(z) spline
//chi_max -> maximum comoving distance to which the integral is computed
//win     -> result is stored here
static int window_lensing(double chi,ccl_cosmology *cosmo,SplPar *spl_pz,double chi_max,double *win)
{
  int gslstatus =0, status =0;
  double result,eresult;
  IntLensPar ip;
  gsl_function F;
  gsl_integration_workspace *w=gsl_integration_workspace_alloc(ccl_gsl->N_ITERATION);

  ip.chi=chi;
  ip.cosmo=cosmo;
  ip.spl_pz=spl_pz;
  ip.status = &status;
  F.function=&integrand_wl;
  F.params=&ip;
  // This conputes the lensing kernel:
  //   w_L(chi) = Integral[ dN/dchi(chi') * f(chi'-chi)/f(chi') , chi < chi' < chi_horizon ]
  // Where f(chi) is the comoving angular distance (which is just chi for zero curvature).
  gslstatus=gsl_integration_qag(&F, chi, chi_max, 0,
                                ccl_gsl->INTEGRATION_EPSREL, ccl_gsl->N_ITERATION,
                                ccl_gsl->INTEGRATION_GAUSS_KRONROD_POINTS,
                                w, &result, &eresult);
  *win=result;
  gsl_integration_workspace_free(w);
  if(gslstatus!=GSL_SUCCESS || *ip.status) {
    ccl_raise_gsl_warning(gslstatus, "ccl_cls.c: window_lensing():");
    return 1;
  }
  //TODO: chi_max should be changed to chi_horizon
  //we should precompute this quantity and store it in cosmo by default

  return 0;
}

//Params for lensing kernel integrand
typedef struct {
  double chi;
  SplPar *spl_pz;
  SplPar *spl_sz;
  ccl_cosmology *cosmo;
  int *status;
} IntMagPar;

//Integrand for magnification kernel
static double integrand_mag(double chip,void *params)
{
  IntMagPar *p=(IntMagPar *)params;
  double chi=p->chi;
  double a=ccl_scale_factor_of_chi(p->cosmo,chip, p->status);
  double z=1./a-1;
  double pz=ccl_spline_eval(z,p->spl_pz);
  double sz=ccl_spline_eval(z,p->spl_sz);
  double h=p->cosmo->params.h*ccl_h_over_h0(p->cosmo,a, p->status)/CLIGHT_HMPC;

  if(chi==0)
    return h*pz*(1-2.5*sz);
  else
    return h*pz*(1-2.5*sz)*ccl_sinn(p->cosmo,chip-chi,p->status)/ccl_sinn(p->cosmo,chip,p->status);
}

//Integral to compute magnification window function
//chi     -> comoving distance
//cosmo   -> ccl_cosmology object
//spl_pz  -> normalized N(z) spline
//spl_pz  -> magnification bias s(z)
//chi_max -> maximum comoving distance to which the integral is computed
//win     -> result is stored here
static int window_magnification(double chi,ccl_cosmology *cosmo,SplPar *spl_pz,SplPar *spl_sz,
				double chi_max,double *win)
{
  int gslstatus =0, status =0;
  double result,eresult;
  IntMagPar ip;
  gsl_function F;
  gsl_integration_workspace *w=gsl_integration_workspace_alloc(ccl_gsl->N_ITERATION);

  ip.chi=chi;
  ip.cosmo=cosmo;
  ip.spl_pz=spl_pz;
  ip.spl_sz=spl_sz;
  ip.status = &status;
  F.function=&integrand_mag;
  F.params=&ip;
  // This conputes the magnification lensing kernel:
  //   w_M(chi) = Integral[ dN/dchi(chi') * (1-5/2 * s(chi)) * f(chi'-chi)/f(chi') , chi < chi' < chi_horizon ]
  // Where f(chi) is the comoving angular distance (which is just chi for zero curvature)
  // and s(chi) is the magnification bias parameter.
  gslstatus=gsl_integration_qag(&F, chi, chi_max, 0,
                                ccl_gsl->INTEGRATION_EPSREL, ccl_gsl->N_ITERATION,
                                ccl_gsl->INTEGRATION_GAUSS_KRONROD_POINTS,
                                w, &result, &eresult);
  *win=result;
  gsl_integration_workspace_free(w);
  if(gslstatus!=GSL_SUCCESS || *ip.status) {
    ccl_raise_gsl_warning(gslstatus, "ccl_cls.c: window_magnification():");
    return 1;
  }
  //TODO: chi_max should be changed to chi_horizon
  //we should precompute this quantity and store it in cosmo by default

  return 0;
}

static void clt_init_nz(CCL_ClTracer *clt,ccl_cosmology *cosmo,
			int nz_n,double *z_n,double *n,int *status)
{
  int gslstatus;
  gsl_function F;
  double nz_norm,nz_enorm;
  double *nz_normalized;

  //Find redshift range where the N(z) has support
  get_support_interval(nz_n,z_n,n,CCL_FRAC_RELEVANT,&(clt->zmin),&(clt->zmax));
  clt->chimax=ccl_comoving_radial_distance(cosmo,1./(1+clt->zmax),status);
  clt->chimin=ccl_comoving_radial_distance(cosmo,1./(1+clt->zmin),status);
  clt->spl_nz=ccl_spline_init(nz_n,z_n,n,0,0);
  if(clt->spl_nz==NULL) {
    *status=CCL_ERROR_SPLINE;
    ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: clt_init_nz(): error initializing spline for N(z)\n");
  }

  if(*status==0) {
    //Normalize n(z)
    nz_normalized=(double *)malloc(nz_n*sizeof(double));
    if(nz_normalized==NULL) {
      *status=CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: clt_init_nz(): memory allocation\n");
      return;
    }
  }

  if(*status==0) {
    gsl_integration_workspace *w=gsl_integration_workspace_alloc(ccl_gsl->N_ITERATION);
    F.function=&speval_bis;
    F.params=clt->spl_nz;
    //Here we're just integrating the N(z) to normalize it to unit probability.
    gslstatus=gsl_integration_qag(&F, z_n[0], z_n[nz_n-1], 0,
				  ccl_gsl->INTEGRATION_EPSREL, ccl_gsl->N_ITERATION,
				  ccl_gsl->INTEGRATION_GAUSS_KRONROD_POINTS,
				  w, &nz_norm, &nz_enorm);
    gsl_integration_workspace_free(w);
    if(gslstatus!=GSL_SUCCESS) {
      ccl_raise_gsl_warning(gslstatus, "ccl_cls.c: clt_init_nz():");
      *status=CCL_ERROR_INTEG;
      ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: clt_init_nz(): integration error when normalizing N(z)\n");
    }
  }

  if(*status==0) {
    for(int ii=0;ii<nz_n;ii++)
      nz_normalized[ii]=n[ii]/nz_norm;
    ccl_spline_free(clt->spl_nz);
    clt->spl_nz=ccl_spline_init(nz_n,z_n,nz_normalized,0,0);
    if(clt->spl_nz==NULL) {
      *status=CCL_ERROR_SPLINE;
      ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: clt_init_nz(): error initializing normalized spline for N(z)\n");
    }
  }

  free(nz_normalized);
}

static void clt_init_bz(CCL_ClTracer *clt,ccl_cosmology *cosmo,
			int nz_b,double *z_b,double *b,int *status)
{
  //Initialize bias spline
  clt->spl_bz=ccl_spline_init(nz_b,z_b,b,b[0],b[nz_b-1]);
  if(clt->spl_bz==NULL) {
    *status=CCL_ERROR_SPLINE;
    ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: clt_init_bz(): error initializing spline for b(z)\n");
  }
}

static void clt_init_wM(CCL_ClTracer *clt,ccl_cosmology *cosmo,
			int nz_s,double *z_s,double *s,int *status)
{
  //Compute magnification kernel
  int nchi;
  double *x,*y;
  double dchi_here=5.;
  double zmax=clt->spl_nz->xf;
  double chimax=ccl_comoving_radial_distance(cosmo,1./(1+zmax),status);
  //TODO: The interval in chi (5. Mpc) should be made a macro

  //In this case we need to integrate all the way to z=0. Reset zmin and chimin
  clt->zmin=0;
  clt->chimin=0;
  clt->spl_sz=ccl_spline_init(nz_s,z_s,s,s[0],s[nz_s-1]);
  if(clt->spl_sz==NULL) {
    *status=CCL_ERROR_SPLINE;
    ccl_cosmology_set_status_message(cosmo,
				     "ccl_cls.c: clt_init_wM(): error initializing spline for s(z)\n");
  }

  if(*status==0) {
    nchi=(int)(chimax/dchi_here)+1;
    x=ccl_linear_spacing(0.,chimax,nchi);
    dchi_here=chimax/nchi;
    if(x==NULL || (fabs(x[0]-0)>1E-5) || (fabs(x[nchi-1]-chimax)>1e-5)) {
      *status=CCL_ERROR_LINSPACE;
      ccl_cosmology_set_status_message(cosmo,
				       "ccl_cls.c: clt_init_wM(): Error creating linear spacing in chi\n");
    }
  }

  if(*status==0) {
    y=(double *)malloc(nchi*sizeof(double));
    if(y==NULL) {
      *status=CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: clt_init_wM(): memory allocation\n");
    }
  }

  if(*status==0) {
    int clstatus=0;
    for(int j=0;j<nchi;j++)
      clstatus|=window_magnification(x[j],cosmo,clt->spl_nz,clt->spl_sz,chimax,&(y[j]));
    if(clstatus) {
      *status=CCL_ERROR_INTEG;
      ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: clt_init_wM(): error computing lensing window\n");
    }
  }

  if(*status==0) {
    clt->spl_wM=ccl_spline_init(nchi,x,y,y[0],0);
    if(clt->spl_wM==NULL) {
      *status=CCL_ERROR_SPLINE;
      ccl_cosmology_set_status_message(cosmo,
				       "ccl_cls.c: clt_init_wM(): error initializing spline for lensing window\n");
    }
  }
  free(x); free(y);
}

//CCL_ClTracer initializer for number counts
static void clt_nc_init(CCL_ClTracer *clt,ccl_cosmology *cosmo,
			int has_rsd,int has_magnification,
			int nz_n,double *z_n,double *n,
			int nz_b,double *z_b,double *b,
			int nz_s,double *z_s,double *s,int *status)
{
  clt->has_rsd=has_rsd;
  clt->has_magnification=has_magnification;

  if ( ((cosmo->params.N_nu_mass)>0) && clt->has_rsd){
    *status=CCL_ERROR_NOT_IMPLEMENTED;
    ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: ccl_cl_tracer_new(): Number counts tracers with RSD not yet implemented in cosmologies with massive neutrinos.");
    return;
  }

  clt_init_nz(clt,cosmo,nz_n,z_n,n,status);
  clt_init_bz(clt,cosmo,nz_b,z_b,b,status);
  if(clt->has_magnification)
    clt_init_wM(clt,cosmo,nz_s,z_s,s,status);
}

static void clt_init_wL(CCL_ClTracer *clt,ccl_cosmology *cosmo,
			int *status)
{
  //Compute weak lensing kernel
  int nchi;
  double *x,*y;
  double dchi_here=5.;
  double zmax=clt->spl_nz->xf;
  double chimax=ccl_comoving_radial_distance(cosmo,1./(1+zmax),status);
  //TODO: The interval in chi (5. Mpc) should be made a macro

  //In this case we need to integrate all the way to z=0. Reset zmin and chimin
  clt->zmin=0;
  clt->chimin=0;
  nchi=(int)(chimax/dchi_here)+1;
  x=ccl_linear_spacing(0.,chimax,nchi);
  dchi_here=chimax/nchi;
  if(x==NULL || (fabs(x[0]-0)>1E-5) || (fabs(x[nchi-1]-chimax)>1e-5)) {
    *status=CCL_ERROR_LINSPACE;
    ccl_cosmology_set_status_message(cosmo,
				     "ccl_cls.c: clt_init_wL(): Error creating linear spacing in chi\n");
  }

  if(*status==0) {
    y=(double *)malloc(nchi*sizeof(double));
    if(y==NULL) {
      *status=CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: clt_init_wL(): memory allocation\n");
    }
  }

  if(*status==0) {
    int clstatus=0;
    for(int j=0;j<nchi;j++)
      clstatus|=window_lensing(x[j],cosmo,clt->spl_nz,chimax,&(y[j]));
    if(clstatus) {
      *status=CCL_ERROR_INTEG;
      ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: clt_init_wL(): error computing lensing window\n");
    }
  }

  if(*status==0) {
    clt->spl_wL=ccl_spline_init(nchi,x,y,y[0],0);
    if(clt->spl_wL==NULL) {
      *status=CCL_ERROR_SPLINE;
      ccl_cosmology_set_status_message(cosmo,
				     "ccl_cls.c: clt_init_wL(): error initializing spline for lensing window\n");
    }
  }
  free(x); free(y);
}

static void clt_init_rf(CCL_ClTracer *clt,ccl_cosmology *cosmo,
			int nz_rf,double *z_rf,double *rf,int *status)
{
  //Initialize bias spline
  clt->spl_rf=ccl_spline_init(nz_rf,z_rf,rf,rf[0],rf[nz_rf-1]);
  if(clt->spl_rf==NULL) {
    *status=CCL_ERROR_SPLINE;
    ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: clt_init_rf(): error initializing spline for b(z)\n");
  }
}

static void clt_init_ba(CCL_ClTracer *clt,ccl_cosmology *cosmo,
			int nz_ba,double *z_ba,double *ba,int *status)
{
  //Initialize bias spline
  clt->spl_ba=ccl_spline_init(nz_ba,z_ba,ba,ba[0],ba[nz_ba-1]);
  if(clt->spl_ba==NULL) {
    *status=CCL_ERROR_SPLINE;
    ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: clt_init_ba(): error initializing spline for b(z)\n");
  }
}

static void clt_wl_init(CCL_ClTracer *clt,ccl_cosmology *cosmo,
			int has_intrinsic_alignment,
			int nz_n,double *z_n,double *n,
			int nz_ba,double *z_ba,double *ba,
			int nz_rf,double *z_rf,double *rf,int *status)
{
  clt->has_intrinsic_alignment=has_intrinsic_alignment;

  clt_init_nz(clt,cosmo,nz_n,z_n,n,status);
  clt_init_wL(clt,cosmo,status);
  if(clt->has_intrinsic_alignment) {
    clt_init_rf(clt,cosmo,nz_rf,z_rf,rf,status);
    clt_init_ba(clt,cosmo,nz_ba,z_ba,ba,status);
  }
}

//CCL_ClTracer creator
//cosmo   -> ccl_cosmology object
//tracer_type -> type of tracer. Supported: ccl_number_counts_tracer, ccl_weak_lensing_tracer
//nz_n -> number of points for N(z)
//z_n  -> array of z-values for N(z)
//n    -> corresponding N(z)-values. Normalization is irrelevant
//        N(z) will be set to zero outside the range covered by z_n
//nz_b -> number of points for b(z)
//z_b  -> array of z-values for b(z)
//b    -> corresponding b(z)-values.
//        b(z) will be assumed constant outside the range covered by z_n
static CCL_ClTracer *cl_tracer(ccl_cosmology *cosmo,int tracer_type,
			       int has_rsd,int has_magnification,int has_intrinsic_alignment,
			       int nz_n,double *z_n,double *n,
			       int nz_b,double *z_b,double *b,
			       int nz_s,double *z_s,double *s,
			       int nz_ba,double *z_ba,double *ba,
			       int nz_rf,double *z_rf,double *rf,
			       double z_source, int * status)
{
  int clstatus=0;
  CCL_ClTracer *clt=(CCL_ClTracer *)malloc(sizeof(CCL_ClTracer));
  if(clt==NULL) {
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: ccl_cl_tracer(): memory allocation\n");
  }

  if(*status==0) {
    clt->tracer_type=tracer_type;

    double hub=cosmo->params.h*ccl_h_over_h0(cosmo,1.,status)/CLIGHT_HMPC;
    clt->prefac_lensing=1.5*hub*hub*cosmo->params.Omega_m;

    if(tracer_type==ccl_number_counts_tracer)
      clt_nc_init(clt,cosmo,has_rsd,has_magnification,
		  nz_n,z_n,n,nz_b,z_b,b,nz_s,z_s,s,status);
    else if(tracer_type==ccl_weak_lensing_tracer)
      clt_wl_init(clt,cosmo,has_intrinsic_alignment,
		  nz_n,z_n,n,nz_ba,z_ba,ba,nz_rf,z_rf,rf,status);
    else if(tracer_type==ccl_cmb_lensing_tracer) {
      clt->chi_source=ccl_comoving_radial_distance(cosmo,1./(1+z_source),status);
      clt->chimax=clt->chi_source;
      clt->chimin=0;
    }
    else {
      free(clt);
      *status=CCL_ERROR_INCONSISTENT;
      ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: ccl_cl_tracer(): unknown tracer type\n");
      return NULL;
    }
  }

  if(*status) {
    free(clt);
    clt=NULL;
  }

  return clt;
}

//CCL_ClTracer constructor with error checking
//cosmo   -> ccl_cosmology object
//tracer_type -> type of tracer. Supported: ccl_number_counts_tracer, ccl_weak_lensing_tracer
//nz_n -> number of points for N(z)
//z_n  -> array of z-values for N(z)
//n    -> corresponding N(z)-values. Normalization is irrelevant
//        N(z) will be set to zero outside the range covered by z_n
//nz_b -> number of points for b(z)
//z_b  -> array of z-values for b(z)
//b    -> corresponding b(z)-values.
//        b(z) will be assumed constant outside the range covered by z_n
CCL_ClTracer *ccl_cl_tracer(ccl_cosmology *cosmo,int tracer_type,
				int has_rsd,int has_magnification,int has_intrinsic_alignment,
				int nz_n,double *z_n,double *n,
				int nz_b,double *z_b,double *b,
				int nz_s,double *z_s,double *s,
				int nz_ba,double *z_ba,double *ba,
				int nz_rf,double *z_rf,double *rf,
				double z_source, int * status)
{
  CCL_ClTracer *clt=cl_tracer(cosmo,tracer_type,has_rsd,has_magnification,has_intrinsic_alignment,
			      nz_n,z_n,n,nz_b,z_b,b,nz_s,z_s,s,
			      nz_ba,z_ba,ba,nz_rf,z_rf,rf,z_source,status);
  ccl_check_status(cosmo,status);
  return clt;
}

//CCL_ClTracer destructor
void ccl_cl_tracer_free(CCL_ClTracer *clt)
{
  if((clt->tracer_type==ccl_number_counts_tracer) || (clt->tracer_type==ccl_weak_lensing_tracer))
    ccl_spline_free(clt->spl_nz);

  if(clt->tracer_type==ccl_number_counts_tracer) {
    ccl_spline_free(clt->spl_bz);
    if(clt->has_magnification) {
      ccl_spline_free(clt->spl_sz);
      ccl_spline_free(clt->spl_wM);
    }
  }
  else if(clt->tracer_type==ccl_weak_lensing_tracer) {
    ccl_spline_free(clt->spl_wL);
    if(clt->has_intrinsic_alignment) {
      ccl_spline_free(clt->spl_ba);
      ccl_spline_free(clt->spl_rf);
    }
  }
  free(clt);
}

CCL_ClTracer *ccl_cl_tracer_cmblens(ccl_cosmology *cosmo,double z_source,int *status)
{
  return ccl_cl_tracer(cosmo,ccl_cmb_lensing_tracer,
			   0,0,0,
			   0,NULL,NULL,0,NULL,NULL,0,NULL,NULL,
			   0,NULL,NULL,0,NULL,NULL,z_source,status);
}

CCL_ClTracer *ccl_cl_tracer_number_counts(ccl_cosmology *cosmo,
					      int has_rsd,int has_magnification,
					      int nz_n,double *z_n,double *n,
					      int nz_b,double *z_b,double *b,
					      int nz_s,double *z_s,double *s, int * status)
{
  return ccl_cl_tracer(cosmo,ccl_number_counts_tracer,has_rsd,has_magnification,0,
			   nz_n,z_n,n,nz_b,z_b,b,nz_s,z_s,s,
			   -1,NULL,NULL,-1,NULL,NULL,0, status);
}

CCL_ClTracer *ccl_cl_tracer_number_counts_simple(ccl_cosmology *cosmo,
						     int nz_n,double *z_n,double *n,
						     int nz_b,double *z_b,double *b, int * status)
{
  return ccl_cl_tracer(cosmo,ccl_number_counts_tracer,0,0,0,
			   nz_n,z_n,n,nz_b,z_b,b,-1,NULL,NULL,
			   -1,NULL,NULL,-1,NULL,NULL,0, status);
}

CCL_ClTracer *ccl_cl_tracer_lensing(ccl_cosmology *cosmo,
					int has_alignment,
					int nz_n,double *z_n,double *n,
					int nz_ba,double *z_ba,double *ba,
					int nz_rf,double *z_rf,double *rf, int * status)
{
  return ccl_cl_tracer(cosmo,ccl_weak_lensing_tracer,0,0,has_alignment,
			   nz_n,z_n,n,-1,NULL,NULL,-1,NULL,NULL,
			   nz_ba,z_ba,ba,nz_rf,z_rf,rf,0, status);
}

CCL_ClTracer *ccl_cl_tracer_lensing_simple(ccl_cosmology *cosmo,
					       int nz_n,double *z_n,double *n, int * status)
{
  return ccl_cl_tracer(cosmo,ccl_weak_lensing_tracer,0,0,0,
			   nz_n,z_n,n,-1,NULL,NULL,-1,NULL,NULL,
			   -1,NULL,NULL,-1,NULL,NULL,0, status);
}

static double f_dens(double a,ccl_cosmology *cosmo,CCL_ClTracer *clt, int * status)
{
  double z=1./a-1;
  double pz=ccl_spline_eval(z,clt->spl_nz);
  double bz=ccl_spline_eval(z,clt->spl_bz);
  double h=cosmo->params.h*ccl_h_over_h0(cosmo,a,status)/CLIGHT_HMPC;

  return pz*bz*h;
}

static double f_rsd(double a,ccl_cosmology *cosmo,CCL_ClTracer *clt, int * status)
{
  double z=1./a-1;
  double pz=ccl_spline_eval(z,clt->spl_nz);
  double fg=ccl_growth_rate(cosmo,a,status);
  double h=cosmo->params.h*ccl_h_over_h0(cosmo,a,status)/CLIGHT_HMPC;

  return pz*fg*h;
}

static double f_mag(double a,double chi,ccl_cosmology *cosmo,CCL_ClTracer *clt, int * status)
{
  double wM=ccl_spline_eval(chi,clt->spl_wM);

  if(wM<=0)
    return 0;
  else
    return wM/(a*chi);
}

//Transfer function for number counts
//l -> angular multipole
//k -> wavenumber modulus
//cosmo -> ccl_cosmology object
//w -> CCL_ClWorskpace object
//clt -> CCL_ClTracer object (must be of the ccl_number_counts_tracer type)
static double transfer_nc(int l,double k,
			  ccl_cosmology *cosmo,CCL_ClWorkspace *w,CCL_ClTracer *clt, int * status)
{
  double ret=0;
  double x0=(l+0.5);
  double chi0=x0/k;
  if(chi0<=clt->chimax) {
    double a0=ccl_scale_factor_of_chi(cosmo,chi0,status);
    double f_all=f_dens(a0,cosmo,clt,status);
    if(clt->has_rsd) {
      double x1=(l+1.5);
      double chi1=x1/k;
      if(chi1<=clt->chimax) {
	double a1=ccl_scale_factor_of_chi(cosmo,chi1,status);
	double pk0=ccl_nonlin_matter_power(cosmo,k,a0,status);
	double pk1=ccl_nonlin_matter_power(cosmo,k,a1,status);
	double fg0=f_rsd(a0,cosmo,clt,status);
	double fg1=f_rsd(a1,cosmo,clt,status);
	f_all+=fg0*(1.-l*(l-1.)/(x0*x0))-fg1*2.*sqrt((l+0.5)*pk1/((l+1.5)*pk0))/x1;
      }
    }
    if(clt->has_magnification)
      f_all+=-2*clt->prefac_lensing*l*(l+1)*f_mag(a0,chi0,cosmo,clt,status)/(k*k);
    ret=f_all;
  }

  return ret;
}

static double f_lensing(double a,double chi,ccl_cosmology *cosmo,CCL_ClTracer *clt, int * status)
{
  double wL=ccl_spline_eval(chi,clt->spl_wL);

  if(wL<=0)
    return 0;
  else
    return clt->prefac_lensing*wL/(a*chi);
}

static double f_IA_NLA(double a,double chi,ccl_cosmology *cosmo,CCL_ClTracer *clt, int * status)
{
  if(chi<=1E-10)
    return 0;
  else {
    double a=ccl_scale_factor_of_chi(cosmo,chi, status);
    double z=1./a-1;
    double pz=ccl_spline_eval(z,clt->spl_nz);
    double ba=ccl_spline_eval(z,clt->spl_ba);
    double rf=ccl_spline_eval(z,clt->spl_rf);
    double h=cosmo->params.h*ccl_h_over_h0(cosmo,a,status)/CLIGHT_HMPC;

    return pz*ba*rf*h/(chi*chi);
  }
}

//Transfer function for shear
//l -> angular multipole
//k -> wavenumber modulus
//cosmo -> ccl_cosmology object
//w -> CCL_ClWorskpace object
//clt -> CCL_ClTracer object (must be of the ccl_weak_lensing_tracer type)
static double transfer_wl(int l,double k,
			  ccl_cosmology *cosmo,CCL_ClWorkspace *w,CCL_ClTracer *clt, int * status)
{
  double ret=0;
  double chi=(l+0.5)/k;
  if(chi<=clt->chimax) {
    double a=ccl_scale_factor_of_chi(cosmo,chi,status);
    double f_all=f_lensing(a,chi,cosmo,clt,status);
    if(clt->has_intrinsic_alignment)
      f_all+=f_IA_NLA(a,chi,cosmo,clt,status);

    ret=f_all;
  }

  return sqrt((l+2.)*(l+1.)*l*(l-1.))*ret/(k*k);
  //return (l+1.)*l*ret/(k*k);
}

static double transfer_cmblens(int l,double k,ccl_cosmology *cosmo,CCL_ClTracer *clt,int *status)
{
  double chi=(l+0.5)/k;
  if(chi>=clt->chi_source)
    return 0;

  if(chi<=clt->chimax) {
    double a=ccl_scale_factor_of_chi(cosmo,chi,status);
    double w=1-chi/clt->chi_source;
    return clt->prefac_lensing*l*(l+1.)*w/(a*chi*k*k);
  }
  return 0;
}

//Wrapper for transfer function
//l -> angular multipole
//k -> wavenumber modulus
//cosmo -> ccl_cosmology object
//clt -> CCL_ClTracer object
static double transfer_wrap(int il,double lk,ccl_cosmology *cosmo,
			    CCL_ClWorkspace *w,CCL_ClTracer *clt, int * status)
{
  double transfer_out=0;
  double k=pow(10.,lk);

  if(clt->tracer_type==ccl_number_counts_tracer)
    transfer_out=transfer_nc(w->l_arr[il],k,cosmo,w,clt,status);
  else if(clt->tracer_type==ccl_weak_lensing_tracer)
    transfer_out=transfer_wl(w->l_arr[il],k,cosmo,w,clt,status);
  else if(clt->tracer_type==ccl_cmb_lensing_tracer)
    transfer_out=transfer_cmblens(w->l_arr[il],k,cosmo,clt,status);
  else
    transfer_out=-1;
  return transfer_out;
}

//Params for power spectrum integrand
typedef struct {
  int il;
  ccl_cosmology *cosmo;
  CCL_ClWorkspace *w;
  CCL_ClTracer *clt1;
  CCL_ClTracer *clt2;
  int *status;
} IntClPar;

//Integrand for integral power spectrum
static double cl_integrand(double lk,void *params)
{
  double d1,d2;
  IntClPar *p=(IntClPar *)params;
  d1=transfer_wrap(p->il,lk,p->cosmo,p->w,p->clt1,p->status);
  if(d1==0)
    return 0;
  d2=transfer_wrap(p->il,lk,p->cosmo,p->w,p->clt2,p->status);
  if(d2==0)
    return 0;

  double k=pow(10.,lk);
  double chi=(p->w->l_arr[p->il]+0.5)/k;
  double a=ccl_scale_factor_of_chi(p->cosmo,chi,p->status);
  double pk=ccl_nonlin_matter_power(p->cosmo,k,a,p->status);

  return k*pk*d1*d2;
}

//Figure out k intervals where the Limber kernel has support
//clt1 -> tracer #1
//clt2 -> tracer #2
//l    -> angular multipole
//lkmin, lkmax -> log10 of the range of scales where the transfer functions have support
static void get_k_interval(ccl_cosmology *cosmo,CCL_ClWorkspace *w,
			   CCL_ClTracer *clt1,CCL_ClTracer *clt2,int l,
			   double *lkmin,double *lkmax)
{
  double chimin,chimax;
  int cut_low_1=0,cut_low_2=0;

  //Define a minimum distance only if no lensing is needed
  if((clt1->tracer_type==ccl_number_counts_tracer) && (clt1->has_magnification==0)) cut_low_1=1;
  if((clt2->tracer_type==ccl_number_counts_tracer) && (clt2->has_magnification==0)) cut_low_2=1;

  if(cut_low_1) {
    if(cut_low_2) {
      chimin=fmax(clt1->chimin,clt2->chimin);
      chimax=fmin(clt1->chimax,clt2->chimax);
    }
    else {
      chimin=clt1->chimin;
      chimax=clt1->chimax;
    }
  }
  else if(cut_low_2) {
    chimin=clt2->chimin;
    chimax=clt2->chimax;
  }
  else {
    chimin=0.5*(l+0.5)/ccl_splines->K_MAX;
    chimax=2*(l+0.5)/ccl_splines->K_MIN;
  }

  if(chimin<=0)
    chimin=0.5*(l+0.5)/ccl_splines->K_MAX;

  *lkmax=log10(fmin( ccl_splines->K_MAX  ,2  *(l+0.5)/chimin));
  *lkmin=log10(fmax( ccl_splines->K_MIN  ,0.5*(l+0.5)/chimax));
}

//Compute angular power spectrum between two bins
//cosmo -> ccl_cosmology object
//il -> index in angular multipole array
//clt1 -> tracer #1
//clt2 -> tracer #2
static double ccl_angular_cl_native(ccl_cosmology *cosmo,CCL_ClWorkspace *cw,int il,
				    CCL_ClTracer *clt1,CCL_ClTracer *clt2,int * status)
{
  int clastatus=0, gslstatus;
  IntClPar ipar;
  double result=0,eresult;
  double lkmin,lkmax;
  gsl_function F;
  gsl_integration_workspace *w=gsl_integration_workspace_alloc(ccl_gsl->N_ITERATION);

  ipar.il=il;
  ipar.cosmo=cosmo;
  ipar.w=cw;
  ipar.clt1=clt1;
  ipar.clt2=clt2;
  ipar.status = &clastatus;
  F.function=&cl_integrand;
  F.params=&ipar;
  get_k_interval(cosmo,cw,clt1,clt2,cw->l_arr[il],&lkmin,&lkmax);
  // This computes the angular power spectra in the Limber approximation between two quantities a and b:
  //  C_ell^ab = 2/(2*ell+1) * Integral[ Delta^a_ell(k) Delta^b_ell(k) * P(k) , k_min < k < k_max ]
  // Note that we use log10(k) as an integration variable, and the ell-dependent prefactor is included
  // at the end of this function.
  gslstatus=gsl_integration_qag(&F, lkmin, lkmax, 0,
                                ccl_gsl->INTEGRATION_LIMBER_EPSREL, ccl_gsl->N_ITERATION,
                                ccl_gsl->INTEGRATION_LIMBER_GAUSS_KRONROD_POINTS,
                                w, &result, &eresult);
  gsl_integration_workspace_free(w);

  // Test if a round-off error occured in the evaluation of the integral
  // If so, try another integration function, more robust but potentially slower
  if(gslstatus == GSL_EROUND) {
    ccl_raise_gsl_warning(gslstatus, "ccl_cls.c: ccl_angular_cl_native(): Default GSL integration failure, attempting backup method.");
    gsl_integration_cquad_workspace *w_cquad= gsl_integration_cquad_workspace_alloc (ccl_gsl->N_ITERATION);
    size_t nevals=0;
    gslstatus=gsl_integration_cquad(&F, lkmin, lkmax, 0,
				    ccl_gsl->INTEGRATION_LIMBER_EPSREL,
				    w_cquad, &result, &eresult, &nevals);
    gsl_integration_cquad_workspace_free(w_cquad);
  }
  if(gslstatus!=GSL_SUCCESS || *ipar.status) {
    ccl_raise_gsl_warning(gslstatus, "ccl_cls.c: ccl_angular_cl_native():");
    // If an error status was already set, don't overwrite it.
    if(*status == 0){
        *status=CCL_ERROR_INTEG;
        ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: ccl_angular_cl_native(): error integrating over k\n");
    }
    return -1;
  }
  ccl_check_status(cosmo,status);

  return result*M_LN10/(cw->l_arr[il]+0.5);
}

void ccl_angular_cls(ccl_cosmology *cosmo,CCL_ClWorkspace *w,
		     CCL_ClTracer *clt1,CCL_ClTracer *clt2,
		     int nl_out,int *l_out,double *cl_out,int *status)
{
  int ii,do_angpow;
  double *l_nodes,*cl_nodes;
  SplPar *spcl_nodes;

  //First check if ell range is within workspace
  for(ii=0;ii<nl_out;ii++) {
    if(l_out[ii]>w->lmax) {
      *status=CCL_ERROR_SPLINE_EV;
      ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: ccl_angular_cls(); "
	     "requested l beyond range allowed by workspace\n");
      return;
    }
  }

  if(*status==0) {
    //Allocate array for power spectrum at interpolation nodes
    l_nodes=(double *)malloc(w->n_ls*sizeof(double));
    if(l_nodes==NULL) {
      *status=CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: ccl_angular_cls(); memory allocation\n");
    }
  }

  if(*status==0) {
    cl_nodes=(double *)malloc(w->n_ls*sizeof(double));
    if(cl_nodes==NULL) {
      *status=CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: ccl_cl_angular_cls(); memory allocation\n");
    }
  }

  if(*status==0) {
    for(ii=0;ii<w->n_ls;ii++)
      l_nodes[ii]=(double)(w->l_arr[ii]);

    do_angpow=0;
    //Now check if angpow is needed at all
    if(w->l_limber>0) {
      for(ii=0;ii<w->n_ls;ii++) {
	if(w->l_arr[ii]<=w->l_limber)
	  do_angpow=1;
      }
    }
#ifndef HAVE_ANGPOW
    do_angpow=0;
#endif //HAVE_ANGPOW

    //Resort to Limber if we have lensing (this will hopefully only be temporary)
    if(clt1->tracer_type==ccl_weak_lensing_tracer || clt2->tracer_type==ccl_weak_lensing_tracer ||
       clt1->has_magnification || clt2->has_magnification) {
      do_angpow=0;
    }

    //Use angpow if non-limber is needed
    if(do_angpow)
      ccl_angular_cls_angpow(cosmo,w,clt1,clt2,cl_nodes,status);
    ccl_check_status(cosmo,status);
  }

  if(*status==0) {
    //Compute limber nodes
    for(ii=0;ii<w->n_ls;ii++) {
      if((!do_angpow) || (w->l_arr[ii]>w->l_limber))
	cl_nodes[ii]=ccl_angular_cl_native(cosmo,w,ii,clt1,clt2,status);
    }

    //Interpolate into ells requested by user
    spcl_nodes=ccl_spline_init(w->n_ls,l_nodes,cl_nodes,0,0);
    if(spcl_nodes==NULL) {
      *status=CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: ccl_cl_angular_cls(); memory allocation\n");
    }
  }

  if(*status==0) {
    for(ii=0;ii<nl_out;ii++)
      cl_out[ii]=ccl_spline_eval((double)(l_out[ii]),spcl_nodes);
  }

  //Cleanup
  ccl_spline_free(spcl_nodes);
  free(cl_nodes);
  free(l_nodes);
}

static int check_clt_fa_inconsistency(CCL_ClTracer *clt,int func_code)
{
  if(((func_code==ccl_trf_nz) && (clt->tracer_type==ccl_cmb_lensing_tracer)) || //lensing has no n(z)
     (((func_code==ccl_trf_bz) || (func_code==ccl_trf_sz) || (func_code==ccl_trf_wM)) &&
      (clt->tracer_type!=ccl_number_counts_tracer)) || //bias and magnification only for clustering
     (((func_code==ccl_trf_rf) || (func_code==ccl_trf_ba) || (func_code==ccl_trf_wL)) &&
      (clt->tracer_type!=ccl_weak_lensing_tracer))) //IAs only for weak lensing
    return 1;
  if((((func_code==ccl_trf_sz) || (func_code==ccl_trf_wM)) &&
      (clt->has_magnification==0)) || //correct combination, but no magnification
     (((func_code==ccl_trf_rf) || (func_code==ccl_trf_ba)) &&
      (clt->has_intrinsic_alignment==0))) //Correct combination, but no IAs
    return 1;
  return 0;
}

double ccl_get_tracer_fa(ccl_cosmology *cosmo,CCL_ClTracer *clt,double a,int func_code,int *status)
{
  SplPar *spl;

  if(check_clt_fa_inconsistency(clt,func_code)) {
    *status=CCL_ERROR_INCONSISTENT;
    ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: inconsistent combination of tracer and internal function to be evaluated");
    return -1;
  }

  switch(func_code) {
  case ccl_trf_nz :
    spl=clt->spl_nz;
    break;
  case ccl_trf_bz :
    spl=clt->spl_bz;
    break;
  case ccl_trf_sz :
    spl=clt->spl_sz;
    break;
  case ccl_trf_rf :
    spl=clt->spl_rf;
    break;
  case ccl_trf_ba :
    spl=clt->spl_ba;
    break;
  case ccl_trf_wL :
    spl=clt->spl_wL;
    break;
  case ccl_trf_wM :
    spl=clt->spl_wM;
    break;
  }

  double x;
  if((func_code==ccl_trf_wL) || (func_code==ccl_trf_wM))
    x=ccl_comoving_radial_distance(cosmo,a,status); //x-variable is comoving distance for lensing kernels
  else
    x=1./a-1; //x-variable is redshift by default

  return ccl_spline_eval(x,spl);
}

int ccl_get_tracer_fas(ccl_cosmology *cosmo,CCL_ClTracer *clt,int na,double *a,double *fa,
		       int func_code,int *status)
{
  SplPar *spl;

  if(check_clt_fa_inconsistency(clt,func_code)) {
    *status=CCL_ERROR_INCONSISTENT;
    ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: inconsistent combination of tracer and internal function to be evaluated");
    return -1;
  }

  switch(func_code) {
  case ccl_trf_nz :
    spl=clt->spl_nz;
    break;
  case ccl_trf_bz :
    spl=clt->spl_bz;
    break;
  case ccl_trf_sz :
    spl=clt->spl_sz;
    break;
  case ccl_trf_rf :
    spl=clt->spl_rf;
    break;
  case ccl_trf_ba :
    spl=clt->spl_ba;
    break;
  case ccl_trf_wL :
    spl=clt->spl_wL;
    break;
  case ccl_trf_wM :
    spl=clt->spl_wM;
    break;
  }

  int compchi = (func_code==ccl_trf_wL) || (func_code==ccl_trf_wM);

  int ia;
  for(ia=0;ia<na;ia++) {
    double x;
    if(compchi) //x-variable is comoving distance for lensing kernels
      x=ccl_comoving_radial_distance(cosmo,a[ia],status);
    else //x-variable is redshift by default
      x=1./a[ia]-1;
    fa[ia]=ccl_spline_eval(x,spl);
  }

  return 0;
}

max_stars_repo_name: vrastil/CCL

max_issues_repo_name: vrastil/CCL

max_forks_repo_name: vrastil/CCL


text:
#define MAIN_FILE

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_spline.h>
#include "global.h"
#include "eos_pres.h"
#include "eos_rho.h"
#include "param.h"

int
main (void)
{
  spline = NULL;
  acc = NULL;
  int i, status;
  double r, r1, y[2];
  // bounds for making grid of m(R) vs. rho(0)
  const double rho_min = 1.0e+13, rho_max = 9.0e+14;
  const int MAX = 100;
  // pressure/density arrays from tabulated EOS data
  int n_eos_pts;
  FILE *fp = NULL;
  double (*pres) (double, void *);	// EOS pressure pointer
  double (*rho) (double, void *);	// EOS density pointer
  double tmp;
  // pointers to tabulated EOS data
  double *eos_tab_pres = NULL, *eos_tab_dens = NULL;

  /* this is the "third arm" of the EOS inversion routine. GSL calls
     the function whose root is being found several times during
     initialization and if we leave the pressure here uninitialized,
     it might be 0 or something crazy, which will be out of bounds of
     the "forward" EOS, and since the forward EOS uses interpolation,
     an out-of-bounds value of pressure will make it crash. */
  tmp_pres = 1.0e+25;

  pres = &eos_pres;
  rho = &eos_rho;

  //  fp = fopen("bck.eos", "r");
  //  fp = fopen("eosC", "r");
  //  fp = fopen("timmes.eos", "r");
  fp = fopen ("../EOS/helmholtz.eos", "r");

  if (fp == NULL)
    {
      fprintf (stderr, "Can't open file!\n");
      exit (1);
    }

  // read in # of EOS data points
  fscanf (fp, "%i", &n_eos_pts);

  eos_tab_dens = (double *) malloc (n_eos_pts * sizeof (double));
  eos_tab_pres = (double *) malloc (n_eos_pts * sizeof (double));

  // read in density and pressure data
  for (i = 0; i < n_eos_pts; i++)
    {
      fscanf (fp, "%le %le", &eos_tab_dens[i], &eos_tab_pres[i]);
    }
  fclose (fp);

  // set up interpolation machinery for EOS
  acc = gsl_interp_accel_alloc ();
  spline = gsl_spline_alloc (gsl_interp_cspline, n_eos_pts);
  gsl_spline_init (spline, eos_tab_dens, eos_tab_pres, n_eos_pts);

  // now make grid
  r = 1.0;			// the integrator will go nuts if we start right at r=0
  r1 = 1.0e+10;			// some final 'radius' (much larger than actual radius)

  printf ("%12s %12s %12s\n", "R", "M(r=R)", "rho(r=0)");

  params.single_star = 1;
  for (i = 0; i < MAX; i++)
    {
      // rho(0) evenly spaced in log
      params.rho_init = log10 (rho_min) +
	(double) i *((log10 (rho_max) - log10 (rho_min)) / (double) MAX);
      params.rho_init = pow (10.0, params.rho_init);

      y[1] = pres (params.rho_init, &params);
      y[0] = (4.0 / 3.0) * M_PI * pow (r, 3.0) * rho (y[1], &params);

      params.pinit = y[1];
      // This function is useful if you want to plot, e.g., central
      // pressure vs. total mass. You can also hang onto the run of
      // pressure with radius, which can be interesting when compared to
      // the Newtonian case.
      status = make_grid (params, r, r1, y);
    }

/*
  params.single_star = 0;
  printf("\nnow for a single star!\n");
  printf("%12s %12s %12s %12s\n", "r", "M(r)", "P(r)", "rho(r)");
  // print P and M vs r for Chandrasekhar-mass star
  params.rho_init = 1.0e+13;

  y[1] = pres(params.rho_init, &params); // rho(0) (roughly) for maximum mass
  y[0] = (4.0/3.0) * M_PI * pow(r, 3.0) * rho(y[1], &params);

  params.pinit = y[1];
  status = make_grid(params, r, r1, y);
*/

  gsl_spline_free (spline);
  gsl_interp_accel_free (acc);

  free (eos_tab_pres);
  free (eos_tab_dens);

  eos_tab_dens = NULL;
  eos_tab_pres = NULL;
  fp = NULL;
  pres = NULL;
  rho = NULL;
  return 0;
}


text:
/* Xrecon.h */
/*---------------------------------------------------------------------------*/
/*                                                                           */
/* Xrecon.h: Global includes and defines                                     */
/*                                                                           */
/* Copyright (C) 2012 Paul Kinchesh                                          */
/*               2012 Martyn Klassen                                         */
/*               2012 Margaret Kritzer                                       */
/*                                                                           */
/* This file is part of Xrecon.                                              */
/*                                                                           */
/* Xrecon 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.                                       */
/*                                                                           */
/* Xrecon 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 Xrecon. If not, see <http://www.gnu.org/licenses/>.            */
/*                                                                           */
/*---------------------------------------------------------------------------*/
/**/

/*--------------------------------------------------------*/
/*---- Wrap the header to prevent multiple inclusions ----*/
/*--------------------------------------------------------*/
#ifndef _H_Xrecon_H
#define _H_Xrecon_H

/*------------------------------------------------------------*/
/*---- If C++ compiler is used it needs to know this is C ----*/
/*------------------------------------------------------------*/
//#ifdef __cplusplus
//extern "C" {
//#endif

/*----------------------------------*/
/*---- Include standard headers ----*/
/*----------------------------------*/
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include <sys/param.h>
#include <sys/stat.h>
#include <stddef.h>
#include <sys/time.h>
#include <sys/types.h>
#include <unistd.h>
#include <signal.h>

/*-----------------------------*/
/*---- Include GSL Headers ----*/
/*-----------------------------*/
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_sf_bessel.h>

/*-----------------------------*/
/*---- Include fftw header ----*/
/*-----------------------------*/
#include <fftw3.h>

/*--------------------------------*/
/*---- Include complex header ----*/
/*--------------------------------*/
/* NB If you have a C compiler, such as gcc, that supports the recent C99 standard,
   and you #include <complex.h> before <fftw3.h>, then fftw_complex is the native
   double-precision complex type and you can manipulate it with ordinary arithmetic.
   Otherwise, FFTW defines its own complex type, which is bit-compatible with the
   C99 complex type. */
/* We will wait for C99 standard to become a bit more standard */
//#include <complex.h>

#ifdef DOTIFF
/*---------------------------------*/
/*---- Include tiff I/O header ----*/
/*---------------------------------*/
#include <tiffio.h>
#endif

/*------------------------------------------*/
/*---- Include Varian data file handler ----*/
/*------------------------------------------*/
#include "data.h"

/*-------------------------------------------*/
/*---- Make sure __FUNCTION__ is defined ----*/
/*-------------------------------------------*/
#ifndef __FUNCTION__
#define __FUNCTION__ __func__
#endif

/*------------------------------*/
/*---- Include bool headers ----*/
/*------------------------------*/
#ifdef _STDC_C99
#include <stdbool.h>
#else
#ifndef bool
#define bool int
#endif
#ifndef true
#define true 1
#endif
#ifndef false
#define false 0
#endif
#endif

/*----------------------------*/
/*---- Some basic defines ----*/
/*----------------------------*/
#define TRUE       1
#define FALSE      0
#define ACTIVE     1
#define NACTIVE    0
#define NACQ       2

/* image scaling */
#define FFT_SCALE 0.5e-6

/*-------------------------*/
/*---- Some data types ----*/
/*-------------------------*/
#define INT16      1  /* 16 bit integer */
#define INT32      2  /* 32 bit integer */
#define FLT32      3  /* 32 bit float */
#define DBL64      4  /* 64 bit double */

/*-------------------------*/
/*---- Some data modes ----*/
/*-------------------------*/
#define NONE       0  /* No data */
#define FID        1  /* Raw */
#define IMAGE      2  /* Image */
#define RE         3  /* Real */
#define IM         4  /* Imaginary */
#define MG         5  /* Magnitude */
#define PH         6  /* Phase */
#define CX         7  /* Complex */
#define STD        8  /* Standard */
#define READ       9  /* Read */
#define PHASE     10  /* Phase */
#define PHASE2    11  /* 2nd Phase */
#define REF       12  /* Reference */
#define REFREAD   13  /* Reference Read */
#define REFPHASE  14  /* Reference Phase */
#define REFPHASE2 15  /* Reference Phase 2 */
#define MK        16  /* Mask */
#define MKREAD    17  /* Mask Read */
#define MKPHASE   18  /* Mask Phase */
#define MKPHASE2  19  /* Mask Phase 2 */
#define RMK       20  /* Reverse Mask of magnitude */
#define MKROI     21  /* Mask ROI */
#define SM        22  /* Sensitivity map */
#define SMREAD    23  /* Sensitivity Map Read */
#define SMPHASE   24  /* Sensitivity Map Phase */
#define SMPHASE2  25  /* Sensitivity Map Phase 2 */
#define GF        26  /* Geometry Factor */
#define RS        27  /* Relative SNR */
#define EPIREF    28  /* EPI reference */
#define VJ        29  /* Defined by VnmrJ parameters */
#define SPATIAL   30  /* for csi */
#define SPECTRAL  31  /* for csi */
#define SPATSPECT 32  /* for csi */
#define PHASE3    33  /* for 3D CSI */

/*---------------------------------------------------------*/
/*---- Some data status bits (0-15) for each dimension ----*/
/*---------------------------------------------------------*/
#define DATA     0x1  /* 0 = no data, 1= data */
#define SHIFT    0x2  /* 0 = no shift, 1 = shift */
#define ZEROFILL 0x4  /* 0 = no zerofill, 1 = zerofill */
#define FFT      0x8  /* 0 = no FFT, 1 = FFT */
/* Bits 4-15 are currently unused */

/*-----------------------------*/
/*---- Some sequence modes ----*/
/*-----------------------------*/
/* For 1D */
#define IM1D      100
#define IM0DCSI   101  /* single voxel csi */

/* For 2D multislice 200 < seqmode < 300 */
#define IM2D      200
#define IM2DCC    201  /* seqcon = "*ccnn" */
#define IM2DCS    202  /* seqcon = "*csnn" */
#define IM2DSC    203  /* seqcon = "*scnn" */
#define IM2DSS    204  /* seqcon = "*ssnn" */
#define IM2DCCFSE 205  /* apptype = "im2Dfse", seqcon = "nccnn" */
#define IM2DCSFSE 206  /* apptype = "im2Dfse", seqcon = "ncsnn" */
#define IM2DSCFSE 207  /* apptype = "im2Dfse", seqcon = "nscnn" */
#define IM2DSSFSE 208  /* apptype = "im2Dfse", seqcon = "nssnn" */
#define IM2DEPI   209  /* apptype = "im2Depi" */
#define IM1DCSI   210  /* csi has one extra dimension */

#define IM2DCCLL  211  /* looklocker seqcon = "*ccnn" */
#define IM2DCSLL  212  /* looklocker seqcon = "*csnn" */
#define IM2DSCLL  213  /* looklocker seqcon = "*scnn" */
#define IM2DSSLL  214  /* looklocker seqcon = "*ssnn" */

/* For 3D 300 < seqmode < 400 */
#define IM3D      300
#define IM3DCC    301  /* seqcon = "**ccn" */
#define IM3DCS    302  /* seqcon = "**csn" */
#define IM3DSC    303  /* seqcon = "**scn" */
#define IM3DSS    304  /* seqcon = "**ssn" */
#define IM3DCFSE  305  /* apptype = "im3Dfse", seqcon = "ncccn" */
#define IM3DSFSE  306  /* apptype = "im3Dfse", seqcon = "nccsn" */
#define IM2DCSCSI   307  /* apptype = "im2DCSI", seqcon = "nncsn" */
#define IM2DSCCSI   308  /* apptype = "im2DCSI", seqcon = "nnscn" */
#define IM2DSSCSI   309  /* apptype = "im2DCSI", seqcon = "ncssn" */
#define IM2DCCCSI   310  /* apptype = "im2DCSI", seqcon = "nnccn" */

/* For 4D */
#define IM4D      400
#define IM3DCCCCSI   401  /* csi has one extra dimension */
#define IM3DCCSCSI   402
#define IM3DCSCCSI   403
#define IM3DCSSCSI   404
#define IM3DSCCCSI   405
#define IM3DSCSCSI   406
#define IM3DSSCCSI   407
#define IM3DSSSCSI   408

 /* max dimensions */
#define MAXDIM    500

/*-------------------------*/
/*---- For TIF scaling ----*/
/*-------------------------*/
#define NSCALE     0  /* Don't scale data */
#define SCALE     -1  /* Scale to maximum of data in each volume */
#define NOISE     -2  /* Scale to see noise */

/*--------------------------------*/
/*---- For NIFTI-1/Analyze7.5 ----*/
/*--------------------------------*/
#define NIFTI      1  /* NIFTI-1 format */
#define ANALYZE    2  /* Analyze7.5 format */

/*---------------------------------------------------------------------------*/
/* Modes for DC correction using dbh.lvl and dbh.tlt in getvol(d,index,mode) */
/*---------------------------------------------------------------------------*/
#define NDCC       0  /* No DC correction using dbh.lvl and dbh.tlt */
#define DCC        1  /* DC correction using dbh.lvl and dbh.tlt */

/*-----------------------------------------------------------------------------------------*/
/* Modes for writing raw data in wrawbin3D(struct data *d,int mode,int type,int precision) */
/*-----------------------------------------------------------------------------------------*/
#define D1         1  /* 1D data mode, data[nr][dim3][dim1] */
#define D12        2  /* 2D/3D data mode, data[nr][dim3][dim2*dim1] */
#define D3         3  /* 3D block processing mode, data[nr][dim2*dim1][dim3] */

/*-------------------*/
/*---- Constants ----*/
/*-------------------*/
#ifndef M_PI
#define M_PI 3.14159265358979323846 /* pi */
#endif
#define DEG2RAD 0.017453292         /* convert degrees to radians */
#define MAXRCVRS 2048               /* maximum number of receivers */

/*------------------------------------------------------*/
/*---- Floating point comparison macros (as in SGL) ----*/
/*------------------------------------------------------*/
/* EPSILON is the largest allowable deviation due to floating point storage */
#define EPSILON 1e-9
#define FP_LT(A,B)  (((A)<(B)) && (fabs((A)-(B))>EPSILON))  /* A less than B */
#define FP_GT(A,B)  (((A)>(B)) && (fabs((A)-(B))>EPSILON))  /* A greater than B */
#define FP_EQ(A,B)  (fabs((A)-(B))<=EPSILON)                /* A equal to B */
#define FP_NEQ(A,B) (!FP_EQ(A,B))                           /* A not equal to B */
#define FP_GTE(A,B) (FP_GT(A,B) || FP_EQ(A,B))              /* A greater than or equal to B */
#define FP_LTE(A,B) (FP_LT(A,B) || FP_EQ(A,B))              /* A less than or equal to B */

/*-------------------------------------------------------*/
/*---- So we can simply include Xrecon.h in Xrecon.c ----*/
/*-------------------------------------------------------*/
#ifdef LOCAL
#define EXTERN
#else
#define EXTERN extern
#endif

int interupt;  // for killing the process

/*-------------------------------------------*/
/*---- Structure to hold input filenames ----*/
/*-------------------------------------------*/
struct file
{
  int      nfiles;    /* number of files */
  char     **fid;     /* fid files */
  char     **procpar; /* procpar files */
};

/*------------------------------------------------*/
/*---- Structure to hold 'procpar' parameters ----*/
/*------------------------------------------------*/
struct pars
{
  int      npars;   /* number of parameters in structure */
  char     **name;  /* parameter name */
  int      *type;   /* parameter type */
  int      *active; /* parameter state */
  int      *nvals;  /* number of values of parameters */
  int      **i;     /* integer values */
  double   **d;     /* real values */
  char     ***s;    /* strings */
};

/*---------------------------------------------------------*/
/*---- Structure to hold maximum value and coordinates ----*/
/*---------------------------------------------------------*/
struct max
{
  int      data;    /* Data flag 0/FALSE=no data, or 1/TRUE=data */
  double   Mval;    /* Maximum magnitude value */
  double   Rval;    /* Maximum real value */
  double   Ival;    /* Maximum imaginary value */
  int      np;      /* position of maximum along np */
  int      nv;      /* position of maximum along nv */
  int      nv2;     /* position of maximum along nv2 */
  int      nv3;     /* position of maximum along nv3 */
};

/*----------------------------------------------*/
/*---- Structure to hold noise measurements ----*/
/*----------------------------------------------*/
struct noise
{
  int      data;    /* Data flag 0/FALSE=no data, or 1/TRUE=data */
  int      zero;    /* Zero data flag 0/FALSE=noise data not zero, or 1/TRUE=noise data is zero */
  int      equal;   /* Equalize flag 0/FALSE=not equal, or 1/TRUE=equal */
  int      samples; /* Number of noise samples */
  double   *M;      /* Mean magnitude */
  double   *M2;     /* Mean magnitude^2 */
  double   *Re;     /* Mean real */
  double   *Im;     /* Mean imaginary */
  double   avM;     /* Averaged over all receivers */
  double   avM2;    /* Averaged over all receivers */
  double   avRe;    /* Averaged over all receivers */
  double   avIm;    /* Averaged over all receivers */
  gsl_matrix_complex   **mat;    /* Noise matrix */
};

/*--------------------------------*/
/*---- Structure to hold data ----*/
/*--------------------------------*/
struct data
{
  FILE                    *fp;         /* fid file pointer */
  FILE                    *datafp;     /* data file pointer */
  FILE                    *phasfp;     /* phasefile file pointer */
  FILE                    *fidfp;      /* for writing in fid format */
  char                    *file;       /* fid file */
  char                    *procpar;    /* procpar file */
  struct datafilehead     fh;          /* data file header */
  struct datablockhead    bh;          /* data block header */
  struct hypercmplxbhead  hcbh;        /* hypercomplex block header */
  struct stat             buf;         /* data file info */

  int                     nvols;       /* number of volumes at end of expt */
  int                     startvol;    /* start volume for processing */
  int                     endvol;      /* end volume for processing */
  int                     vol;         /* current volume counter */
  int                     outvol;      /* output volume (not counting reference volumes) */

  int                     nblocks;     /* number of blocks */
  int                     startpos;    /* start position of block for processing */
  int                     endpos;      /* end position of block for processing */
  int                     block;       /* current block counter */

  int                     np,nv,nv2,nv3; /* number of points and views including nv3 for CSI */
  int                     nseg,etl;    /* number of segments and echo train length */
  int                     ne;          /* number of echoes */
  int                     fn,fn1,fn2,fn3;  /* zerofilling */
  int                     nr,ns;       /* number of receivers and slices */
  int                     korder;       /* for csi, indicates elliptical or otherwise acquired */
  int                     synclist[3];  /* indicates which dimensions are part of sync'ed lists such as pelist, pe2list, etc. */
  int                     startd1, startd2, startd3;  // for csi cropping
  int                     cropd1, cropd2, cropd3;     // for csi cropping

  int                     profile;     /* profile flag, 0/FALSE=not profile, 1/TRUE=profile */
  int                     proj2D;      /* 2D projection flag, 0/FALSE=not 2D projection, 1/TRUE=2D projection */

  int                     nav;         /* navigator flag, 0/FALSE=none, 1/TRUE=navigators */
  int                     nnav;        /* number of navigators */
  int                     *navpos;     /* position(s) of navigator data amongst actual data */

  int                     seqmode;     /* sequence mode depends on seqcon and apptype */

  int                     *dim2order;  /* dim2 order in fid file */
  int                     *dim3order;  /* dim3 order in fid file */
  int                     *dim4order;  /* dim4 order in fid file */
  int                     *pssorder;   /* slice order in fid file */

  int                     *dim2orderR; /* Table for reduced dim2 order */
  int                     *dim3orderR; /* Table for reduced dim3 order */
  int                     *dim4orderR; /* Table for reduced dim4 order */
  int                     dim2R,dim3R, dim4R; /* Size of dim2 and dim3 and dim4 BEFORE reduction */
  int					  d1rev, d2rev, d3rev; /* flag for csi dimension reversal */

  fftw_complex            ***data;     /* actual data */
  fftw_complex	          **csi_data; /* different organization */
  int                     **mask;      /* data mask */

  int                     ndim;        /* number of data dimensions */
  int                     *dimstatus;  /* status of each dimension */

  int                     datamode;    /* data mode to flag different types of data, eg. see default EPI recon */

  struct max              max;         /* maximum value and coordinates */
  struct noise            noise;       /* noise measurement data */

  struct pars             p;           /* 'procpar' parameters */
  struct pars             a;           /* 'array' parameters */
  struct pars             s;           /* 'sviblist' parameters */

  char                    *fdfhdr;     /* fdf header */
};

/*-------------------------------------*/
/*---- Structure to hold dimstatus ----*/
/*-------------------------------------*/
struct dimstatus
{
  int ndim;        /* number of data dimensions */
  int *dimstatus;  /* status of each dimension */
};

/*-------------------------------------------------*/
/*---- Structure to hold segment scaling (EPI) ----*/
/*-------------------------------------------------*/
struct segscale
{
  int data;
  double ***value;
};

/*------------------------------*/
/*---- Main recon functions ----*/
/*------------------------------*/
void recon1D(struct data *d);    /* recon1D.c */
void recon2D(struct data *d);    /* recon2D.c */
void reconEPI(struct data *d);   /* reconEPI.c */
void recon3D(struct data *d);    /* recon3D.c */
void recon2DCSI(struct data *d);

/*---------------------------------*/
/*---- Default recon functions ----*/
/*---------------------------------*/
void default1D(struct data *d);     /* default1D.c */
void default2D(struct data *d);     /* default2D.c */
void defaultEPI(struct data *d);    /* defaultEPI.c */
void default3D(struct data *d);     /* default3D.c */

/*-------------------------------*/
/*---- default1D.c functions ----*/
/*-------------------------------*/
void refcorr1D(struct data *d,struct data *ref);
void combine1D(struct data *d);

/*-------------------------------*/
/*---- profile1D.c functions ----*/
/*-------------------------------*/
void profile1D(struct data *d);

/*----------------------------*/
/*---- proj2D.c functions ----*/
/*----------------------------*/
void proj2D(struct data *d);

/*------------------------*/
/*---- CSI  functions ----*/
/*------------------------*/
void reconCSI(struct data *d);
void reconCSI2D(struct data *d);
void reconCSI3D(struct data *d);
void dimorderCSI2D(struct data *d);
void dimorderCSI3D(struct data *d);
void getblockCSI2D(struct data *d,int volindex,int DCCflag);
void getblockCSI3D(struct data *d,int volindex,int DCCflag);
void w2DCSI(struct data *d,int receiver,int d4, int d3, int d2,int fileid);
void wdfhCSI(struct data *d,int fileid);
void zerofill2DCSI(struct data *d, int mode);
void phaseramp2DCSI(struct data *d, int d4index, int mode);
void shiftdim3dataCSI(struct data *d,int shft);
void phaserampdim3CSI(struct data *d, int mode);
void regridEPSI(struct data *d);

/*----------------------------*/
/*---- mask2D.c functions ----*/
/*----------------------------*/
void mask2D(struct data *d);
void add2mask2D(struct data *d);
int read2Dmask(struct data *d,int mode);
void mask2Ddata(struct data *d1,struct data *d2);

/*-----------------------------*/
/*---- sense2D.c functions ----*/
/*-----------------------------*/
void sense2D(struct data *d);
void sense2Dunfold(struct data *d,struct data *ref);
void svd2Dinmem(struct data *ref);
void setref2Dmatrix(struct data *ref,struct data *d);
int checksenseref(struct data *d,struct data *ref);

/*-----------------------------*/
/*---- sensi2D.c functions ----*/
/*-----------------------------*/
void sensibility2D(struct data *d);
void gmap2D(struct data *d);

/*----------------------------*/
/*---- smap2D.c functions ----*/
/*----------------------------*/
void smap2D(struct data *d,int mode);
int smap2Dvcoil(struct data *d,struct data *ref,int mode);
int smap2Dacoil(struct data *d,int mode);
int setvcoil2D(struct data *d,struct data *ref,struct file *fref);
void gen2Dsmapsos(struct data *d);
void gen2Dsmapsuper(struct data *d);
void gen2Dsmapvcoil(struct data *d,struct data *ref);

/*-------------------------------*/
/*---- asltest2D.c functions ----*/
/*-------------------------------*/
void asltest2D(struct data *d);

/*-----------------------------------*/
/*---- multiblock3D.c  functions ----*/
/*-----------------------------------*/
void multiblock3D(struct data *d);

/*------------------------------*/
/*---- dprocEPI.c functions ----*/
/*------------------------------*/
void addscaledEPIref(struct data *d,struct data *ref1,struct data *ref2);
void addEPIref(struct data *d,struct data *ref);
void ftnvEPI(struct data *d);
void nvfillEPI(struct data *d);
void ftnpEPI(struct data *d);
void zoomEPI(struct data *d);
void revreadEPI(struct data *d);
void prepEPIref(struct data *ref1,struct data *ref2);
void phaseEPIref(struct data *ref1,struct data *ref2,struct data *ref3);
void phaseEPI(struct data *d,struct data *ref);
void phaserampEPI(struct data *d,int mode);
void navcorrEPI(struct data *d);
void stripEPInav(struct data *d);
void weightnavs(struct data *d,double gf);
void setsegscale(struct data *d,struct segscale *scale);
void segscale(struct data *d,struct segscale *scale);
void analyseEPInav(struct data *d);
void getblockEPI(struct data *d,int volindex,int DCCflag);
void setblockEPI(struct data *d);
void setoutvolEPI(struct data *d);
int outvolEPI(struct data *d);
void setnvolsEPI(struct data *d);
void setblockSGE(struct data *d);

/*---------------------------------*/
/*---- prescanEPI.c functions -----*/
/*---------------------------------*/
void prescanEPI(struct data *d);
void settep(struct data *d, int mode);

/*---------------------------*/
/*---- dproc.c functions ----*/
/*---------------------------*/
int *sliceorder(struct data *d,int dim,char *par);
int *phaseorder(struct data *d,int views,int dim,char *par);
int *phaselist(struct data *d,char *par);
void getmax(struct data *d);
void getnoise(struct data *d,int mode);
void equalizenoise(struct data *d,int mode);
void dccorrect(struct data *d);
void scaledata(struct data *d,double factor);
void combine_channels(struct data *d);
void opti_comb(struct data *d, int offset);
double *unwrap1D(double *angles, int nangles);

/*-----------------------------*/
/*---- dproc1D.c functions ----*/
/*-----------------------------*/
void fft1D(struct data *d,int dataorder);
void shiftdata1D(struct data *d,int mode,int dataorder);
void shift1Ddata(struct data *d,int shft,int dataorder);
void weightdata1D(struct data *d,int mode,int dataorder);
void zerofill1D(struct data *d,int mode,int dataorder);
void phaseramp1D(struct data *d,int mode);
void phasedata1D(struct data *d,int mode,int dataorder);
void getmaxtrace1D(struct data *d,int trace);

/*-----------------------------*/
/*---- dproc2D.c functions ----*/
/*-----------------------------*/
void dimorder2D(struct data *d);
void fft2D(struct data *d,int mode);
void ifft2D(struct data *d,int mode);
void shiftdata2D(struct data *d,int mode);
void shift2Ddata(struct data *d,int npshft,int nvshft);
void shift2DCSIdata(struct data *d,int nvshft,int nv2shft);
void weightdata2D(struct data *d,int mode);
void zerofill2D(struct data *d,int mode);
void phaseramp2D(struct data *d,int mode);
void phasedata2D(struct data *d,int mode);
void zoomdata2D(struct data *d,int startdim1, int widthdim1, int startdim2, int widthdim2);
void revread(struct data *d);
void navcorr(struct data *d,struct data *nav);

/*-----------------------------*/
/*---- dproc3D.c functions ----*/
/*-----------------------------*/
void dimorder3D(struct data *d);
void getblock3D(struct data *d,int volindex,int DCCflag);
void getnavblock3D(struct data *d,int volindex,int DCCflag);
void shiftdatadim3(struct data *d,int dataorder,int mode);
void shiftdim3data(struct data *d,int dim3shft);
void shiftD3data(struct data *d,int dim3shft);
void weightdatadim3(struct data *d,int mode);
void zerofilldim3(struct data *d,int mode);
void fftdim3(struct data *d);
void phaserampdim3(struct data *d,int mode);
void phasedatadim3(struct data *d,int mode);

/*-----------------------------*/
/*---- noise2D.c functions ----*/
/*-----------------------------*/
void get2Dnoisematrix(struct data *d,int mode);
void zero2Dnoisematrix(struct data *d);
void print2Dnoisematrix(struct data *d);

/*-----------------------------*/
/*---- dmask2D.c functions ----*/
/*-----------------------------*/
void get2Dmask(struct data *d,int mode);
void fill2Dmask(struct data *d,int mode);

/*---------------------------*/
/*---- dread.c functions ----*/
/*---------------------------*/
void setblock(struct data *d,int dim);
void getblock(struct data *d,int volindex,int DCCflag);
int readfblock(struct data *d,int volindex,int DCCflag);
int readlblock(struct data *d,int volindex,int DCCflag);
int readsblock(struct data *d,int volindex,int DCCflag);
int setoffset(struct data *d,int volindex,int receiver,int dim3index,int dim2index);
void getnavblock(struct data *d,int volindex,int DCCflag);
int readnavfblock(struct data *d,int volindex,int DCCflag);
int readnavlblock(struct data *d,int volindex,int DCCflag);
int readnavsblock(struct data *d,int volindex,int DCCflag);
int setnavoffset(struct data *d,int volindex,int receiver,int dim3index,int seg,int nav);
int sliceindex(struct data *d,int slice);
int segindex(struct data *d,int phase);
int phaseindex(struct data *d,int phase);
int phase2index(struct data *d,int phase2);
int phase3index(struct data *d,int phase3);
void synctablesort2D(struct data *d, int d2, int d3, int *p2, int *p3);
int synctablesort3D(struct data *d, int d2, int d3, int d4, int *p2, int *p3, int *p4);

/*-----------------------------*/
/*---- dread1D.c functions ----*/
/*-----------------------------*/
void getblock1D(struct data *d,int volindex,int DCCflag);

/*-----------------------------*/
/*---- dread2D.c functions ----*/
/*-----------------------------*/
void getblock2D(struct data *d,int volindex,int DCCflag);
void getnavblock2D(struct data *d,int volindex,int DCCflag);
int process2Dblock(struct data *d,char *startpar,char *endpar);

/*---------------------------*/
/*---- dhead.c functions ----*/
/*---------------------------*/
void getdfh(struct data *d);
void getdbh(struct data *d,int blockindex);
void gethcbh(struct data *d,int blockindex);
void reversedfh(struct datafilehead *dfh);
void reversedbh(struct datablockhead *dbh);
void reversehcbh(struct hypercmplxbhead *hcbh);
void printfileheader(struct datafilehead *dfh);
void printfilestatus(struct datafilehead *dfh);
void printblockheader(struct datablockhead *dbh,int blockindex);
void printhcblockheader(struct hypercmplxbhead *hcbh,int blockindex);

/*----------------------------*/
/*---- dutils.c functions ----*/
/*----------------------------*/
void opendata(char *datafile,struct data *d);
void setnvols(struct data *d);
void setstartendvol(struct data *d);
void setdatapars(struct data *d);
void setseqmode(struct data *d);
void setdim(struct data *d);
int spar(struct data *d,char *par,char *str);
int im1D(struct data *d);
int im2D(struct data *d);
int im2DLL(struct data *d);
int im3D(struct data *d);
int im4D(struct data *d);
int imCSI(struct data *d);
double getelem(struct data *d,char *par,int image);
int wprocpar(struct data *d,char *filename);
void copydbh(struct datablockhead *dbh1,struct datablockhead *dbh2);
void setfile(struct file *datafile,char *datadir);
void setreffile(struct file *datafile,struct data *d,char *refpar);
void setfn(struct data *d1,struct data *d2,double multiplier);
void setfn1(struct data *d1,struct data *d2,double multiplier);
void setdimstatus(struct data *d,struct dimstatus *status,int block);
void copynblocks(struct data *d1,struct data *d2);
void copymaskpars(struct data *d1,struct data *d2);
void copysmappars(struct data *d1,struct data *d2);
void copysensepars(struct data *d1,struct data *d2);
int checkequal(struct data *d1,struct data *d2,char *par,char *comment);
void initdatafrom(struct data *d1, struct data *d2);
void initdata(struct data *d);
void nulldata(struct data *d);
void cleardimorder(struct data *d);
void clearstatus(struct data *d);
void zeromax(struct data *d);
void zeronoise(struct data *d);
void initnoise(struct data *d);
void nullnoise(struct data *d);
void closedata(struct data *d);

/*------------------------------*/
/*---- dutils1D.c functions ----*/
/*------------------------------*/
void clear1Dall(struct data *d);
void clear1Ddata(struct data *d);
void clearnoise1D(struct data *d);
void clear1Dmask(struct data *d);

/*------------------------------*/
/*---- dutils2D.c functions ----*/
/*------------------------------*/
int check2Dref(struct data *d,struct data *ref);
void copy2Ddata(struct data *d1,struct data *d2);
void clear2Dall(struct data *d);
void clear2Ddata(struct data *d);
void clearCSIdata(struct data *d);
void clearnoise2D(struct data *d);
void clear2Dmask(struct data *d);
void checkCrop(struct data *d);
void checkCrop3D(struct data *d);

/*------------------------------------*/
/*---- pars.c 'procpar' functions ----*/
/*------------------------------------*/
void getpars(char *procpar,struct data *d);
void getprocparpars(char *procpar,struct pars *p);
void skipprocparvals(FILE *fp,struct pars *p,int type,int nvals);
void procparerror(const char *function,char *par,char *str);
void setarray(struct pars *p,struct pars *a);
int list2array(char *par,struct pars *p1,struct pars *p2);
void setval(struct pars *p,char *par,double value);
void setsval(struct pars *p,char *par,char *svalue);
void cppar(char *par1,struct pars *p1,char *par2,struct pars *p2);
void copypar(char *par,struct pars *p1,struct pars *p2);
void copypars(struct pars *p1, struct pars *p2);
void copyvalues(struct pars *p1,int ix1,struct pars *p2,int ix2);
void mallocnpars(struct pars *p);
void freenpars(struct pars *p);
void printarray(struct pars *a,char *array);
void printprocpar(struct pars *p);
int ptype(char *par,struct pars *p);
int nvals(char *par,struct pars *p);
int *ival(char *par,struct pars *p);
double *val(char *par,struct pars *p);
char **sval(char *par,struct pars *p);
int parindex(char *par,struct pars *p);
int arraycheck(char *par,struct pars *a);
int getcycle(char *par,struct pars *a);
int cyclefaster(char *par1,char *par2,struct pars *a);
void clearpars(struct pars *p);
void cleararray(struct pars *a);
void nullpars(struct pars *p);

/*-------------------------------------*/
/*---- utils.c utilities functions ----*/
/*-------------------------------------*/
int check_CPU_type();
void complexmatrixmultiply(gsl_matrix_complex *A,gsl_matrix_complex *B,gsl_matrix_complex *C);
void complexmatrixconjugate(gsl_matrix_complex *A);
int round2int(double dval);
void reverse2ByteOrder(int nele,char *ptr);
void reverse4ByteOrder(int nele,char *ptr);
void reverse8ByteOrder(int nele,char *ptr);
int doublecmp(const void *double1,const void *double2);
int nomem(char *file,const char *function,int line);
int createdir(char *dirname);
int checkdir(char *dirname);
int checkfile(char *filename);

/*-----------------------------*/
/*---- write1D.c functions ----*/
/*-----------------------------*/
void w1Dtrace(struct data *d,int receiver,int trace,int fileid);
void w1Dblock(struct data *d,int receiver,int fileid);
void wdfh(struct data *d,int fileid);
void wdbh(struct data *d,int fileid);
void openfpw(struct data *d,int fileid);
void closefp(struct data *d,int fileid);

/*--------------------------------*/
/*---- fdfwrite2D.c functions ----*/
/*--------------------------------*/
void w2Dfdfs(struct data *d,int type,int scalevalue,int volindex);
int magnitude2Dfdfs(struct data *d,char *par,char *parIR,char *outdir,int type,int precision,int volindex);
void other2Dfdfs(struct data *d,char *par,char *outdir,int type,int precision,int volindex);
void gen2Dfdfs(struct data *d,int mode,char *outdir,int type,int precision,int volindex);
void w2Dfdf(char *filename,struct data *d,int image,int slice,int sliceindex,int echo,int receiver,int type,int precision);
void wcomb2Dfdf(char *filename,struct data *d,int image,int slice,int sliceindex,int echo,int type,int precision);
void gen2Dfdfhdr(struct data *d,int image,int slice,int echo,int receiver,int type,int precision);
int validtype(int type);
int addpar2hdr(struct pars *p,int ix);
double sliceposition(struct data *d,int slice);

/*--------------------------------*/
/*---- fdfwrite3D.c functions ----*/
/*--------------------------------*/
void w3Dfdfs(struct data *d,int type,int scalevalue,int volindex);
int magnitude3Dfdfs(struct data *d,char *par,char *parIR,char *outdir,int type,int precision,int volindex);
void other3Dfdfs(struct data *d,char *par,char *outdir,int type,int precision,int volindex);
void gen3Dfdfs(struct data *d,int mode,char *outdir,int type,int precision,int volindex);
void w3Dfdf(char *filename,struct data *d,int image,int slab,int echo,int receiver,int type,int precision);
void wcomb3Dfdf(char *filename,struct data *d,int image,int slab,int echo,int type,int precision);
void gen3Dfdfhdr(struct data *d,int image,int slice,int echo,int receiver,int type,int precision);

/*-----------------------------*/
/*---- rawIO3D.c functions ----*/
/*-----------------------------*/
void wrawbin3D(struct data *d,int dataorder,int type,int precision);
void gen3Draw(struct data *d,int output,char *outdir,int dataorder,int type,int precision);
void w3Draw(char *filename,struct data *d,int receiver,int type,int precision);
void w3DrawD3(char *filename,struct data *d,int receiver);
void wcomb3Draw(char *filename,struct data *d,int type,int precision);
void cleardim3data(struct data *d);
void rrawbin3D(struct data *d,int dataorder,int type,int precision);
void get3Draw(struct data *d,int mode,char *indir,int dataorder);
void r3Draw(char *filename,struct data *d,int receiver);
void r3DrawD3(char *filename,struct data *d,int receiver);
void delrawbin3D(struct data *d,int dataorder,int type);
void del3Draw(struct data *d,char *dir);

/*--------------------------------*/
/*---- niftiwrite.c functions ----*/
/*--------------------------------*/
void wnifti(struct data *d,int type,int precision,int volindex);
void writenifti(struct data *d,int type,int precision,int format,int volindex);
int magnifti(struct data *d,char *par,char *parIR,char *outdir,int type,int precision,int format,int volindex);
void othernifti(struct data *d,char *par,char *outdir,int type,int precision,int format,int volindex);
void gennifti(struct data *d,int mode,char *outdir,int type,int precision,int format,int volindex);
void wniftihdr(char *filebase,struct data *d,int precision,int format,int image);
void wniftidata(char *filebase,struct data *d,int receiver,int type,int precision,int format,int newfile);
void wcombniftidata(char *filebase,struct data *d,int type,int precision,int format,int newfile);

/*--------------------------------*/
/*---- tifwrite2D.c functions ----*/
/*--------------------------------*/
int wtifs2D(struct data *d,int type,int scalevalue,int volindex);
int wtif2D(char *filename,struct data *d,int receiver,int slice,int type,double scale);

/*--------------------------------*/
/*---- rawwrite2D.c functions ----*/
/*--------------------------------*/
int wrawbin2D(struct data *d,int type,int precision,int volindex);
int wraw2D(char *filename,struct data *d,int receiver,int slice,int type,int precision);

/*-----------------------------*/
/*---- options.c functions ----*/
/*-----------------------------*/
void getoptions(struct file *f,int argc,char *argv[]);

/*---------------------------*/
/*---- For byte swapping ----*/
/*---------------------------*/
EXTERN int reverse_byte_order;

/*-------------------------*/
/*---- For VnmrJ recon ----*/
/*-------------------------*/
EXTERN int vnmrj_recon;
EXTERN char vnmrj_path[MAXPATHLEN];

/*-------------------------------------------------------------*/
/*---- If C++ compiler is used it needs to know above is C ----*/
/*-------------------------------------------------------------*/
//#ifdef __cplusplus
//}
//#endif

/*----------------------------*/
/*---- End of header wrap ----*/
/*----------------------------*/
#endif

max_stars_repo_name: DanIverson/OpenVnmrJ

max_issues_repo_name: DanIverson/OpenVnmrJ

max_forks_repo_name: DanIverson/OpenVnmrJ