text:
/* integration/rational.c
 *
 * Copyright (C) 2017 Konrad Griessinger, Patrick Alken
 *
 * 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.
 */

/*
 * The code in this module is based on IQPACK, specifically the LGPL
 * implementation found in HERMITE_RULE:
 * https://people.sc.fsu.edu/~jburkardt/c_src/hermite_rule/hermite_rule.html
 */

#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_sf_gamma.h>

static int
rational_check(const size_t n, const gsl_integration_fixed_params * params)
{
  if (fabs(params->b - params->a) <= GSL_DBL_EPSILON)
    {
      GSL_ERROR("|b - a| too small", GSL_EDOM);
    }
  else if (params->alpha <= -1.0)
    {
      GSL_ERROR("alpha must be > -1", GSL_EDOM);
    }
  else if (params->beta >= 0.0 || params->alpha+params->beta+2*n >= 0.0 || 0.0 >= params->alpha+2*n)
    {
      GSL_ERROR("beta < alpha + beta + 2n < 0 is required", GSL_EDOM);
    }
  else if (params->a + params->b <= 0.0)
    {
      GSL_ERROR("a + b <= 0 is not allowed", GSL_EDOM);
    }
  else
    {
      return GSL_SUCCESS;
    }
}

static int
rational_init(const size_t n, double * diag, double * subdiag, gsl_integration_fixed_params * params)
{
  const double absum = params->beta + params->alpha;
  const double a1 = params->alpha + 1.0;
  const double aba1 = absum*a1;
  double ab2i = absum + 2.0;
  size_t i;

  /* construct the diagonal and subdiagonal elements of Jacobi matrix */

  diag[0] = -a1/(absum + 2.0);
  subdiag[0] = sqrt( -diag[0] * ( params->beta + 1.0 ) / ( (absum + 2.0)*(absum + 3.0) ) );

  for (i = 1; i < n-1; i++)
    {
      ab2i += 2.0;
      diag[i] = ( -aba1 - 2.0 * i * ( absum + i + 1.0 ) ) / ( ab2i * ( ab2i - 2.0 ) );
      subdiag[i] = sqrt( (i+1.0) * ( params->alpha + i + 1.0 ) / ( ab2i - 1.0 ) * ( params->beta + i + 1.0 ) / ( ab2i * ab2i ) * ( absum + i + 1.0 ) / ( ab2i + 1.0 ) );
    }

  diag[n-1] = ( -aba1 - 2.0 * (n-1.0) * ( absum + n ) ) / ( (absum + 2.0*n) * ( absum + 2.0*n - 2.0 ) );
  subdiag[n-1] = 0.0;

  params->zemu = gsl_sf_gamma(params->alpha + 1.0) * gsl_sf_gamma(-absum - 1.0) / gsl_sf_gamma(-params->beta);
  params->shft = params->a;
  params->slp = params->b + params->a;
  params->al = params->alpha;
  params->be = params->beta;

  return GSL_SUCCESS;
}

static const gsl_integration_fixed_type rational_type =
{
  rational_check,
  rational_init
};

const gsl_integration_fixed_type *gsl_integration_fixed_rational = &rational_type;

ext: c

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#include <jni.h>
#include <assert.h>
#include <lapacke.h>

/* LU */

JNIEXPORT jint Java_JAMAJni_LUDecomposition_dgetrf (JNIEnv *env, jclass klass, jint matrix_layout, jint m, jint n, jdoubleArray a, jint lda, jintArray ipiv){

    double *aElems;
    lapack_int *ipivElems;
    int info;

    aElems = (*env)-> GetDoubleArrayElements (env, a, NULL);
    ipivElems = (*env)-> GetIntArrayElements (env, ipiv, NULL);

    assert(aElems && ipivElems);

    info = LAPACKE_dgetrf ((int) matrix_layout, (lapack_int) m, (lapack_int) n, aElems, (lapack_int) lda, ipivElems);

    (*env)-> ReleaseDoubleArrayElements (env, a, aElems, 0);
    (*env)-> ReleaseIntArrayElements (env, ipiv, ipivElems, 0);

    return info;
}

JNIEXPORT jint Java_JAMAJni_LUDecomposition_dgetrs (JNIEnv *env, jclass klass, jint matrix_layout, jchar trans, jint n, jint nrhs, jdoubleArray a, jint lda, jintArray ipiv, jdoubleArray b, jint ldb){

    double *bElems;
    double *aElems;
    lapack_int *ipivElems;
    int info;

    aElems = (*env)-> GetDoubleArrayElements (env, a, NULL);
    bElems = (*env)-> GetDoubleArrayElements (env, b, NULL);
    ipivElems = (*env)-> GetIntArrayElements (env, ipiv, NULL);

    assert(aElems && bElems && ipivElems);

    info = LAPACKE_dgetrs ((int) matrix_layout, (char) trans, (lapack_int) n, (lapack_int) nrhs, aElems, (lapack_int) lda, ipivElems, bElems, (lapack_int) ldb);

    (*env)-> ReleaseDoubleArrayElements (env, a, aElems, JNI_ABORT);
    (*env)-> ReleaseDoubleArrayElements (env, b, bElems, 0);
    (*env)-> ReleaseIntArrayElements (env, ipiv, ipivElems, JNI_ABORT);

    return info;
}


text:
#include <stdio.h>
#include <errno.h>
#include <string.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_permutation.h>
#include "matrixUtils.h"
#include "../mcmc/mvnormMCMC/mvrandist.h"

/*A collection of matrix functions*/
/*Alexander G. Shanku*/
/*Wed 12/05/2012 12:44:41 EST*/

//// Makes a matrix filled with random U(0,1) elements ////
gsl_matrix fill_rand_matrix(gsl_matrix *X, gsl_rng *r_num){
	int i,j;
	/*int step = 0;*/
	int n = X->size1;
	int m = X->size2;
	for (i = 0; i < m; ++i){
		for (j = 0; j < n; ++j){
			gsl_matrix_set(X, i, j, gsl_rng_uniform(r_num));
		}
	}
	return(*X);
}

//// Makes a matrix filled with sequential elements, starting with 0 ////
gsl_matrix fill_matrix(gsl_matrix *X){
	int i,j;
	int n = X->size1;
	int m = X->size2;
	for (i = 0; i < m; ++i){
		for (j = 0; j < n; ++j){
			gsl_matrix_set(X, i, j, i+(j*m));
		}
	}
	return(*X);
}

//// Do LU decomp, get determinant and keep input matrix untouched ////
double matrix_determ(gsl_matrix *X){
	int s;
	int n = X->size1;
	int m = X->size2;
	gsl_matrix *a_copy = gsl_matrix_alloc(m, n);
	gsl_matrix_memcpy(a_copy, X );
	gsl_permutation *P = gsl_permutation_alloc(n);
	gsl_linalg_LU_decomp(a_copy, P, &s);
	double my_det = gsl_linalg_LU_det (a_copy, s);
	return(my_det);
}

//// Get trace of given matrix ////
double matrix_trace(gsl_matrix *X){
	int i, m;
	m = X->size1;
	double trace = 0.0;
	for(i=0;i<m;i++){
		trace += gsl_matrix_get(X,i,i);
	}
	return(trace);
}

//// Multivariate Gamma function ////
double mv_gamma(double a, double d){
	double val = 1.0;
	int i;
		for(i = 1; i <= d; i++){
		 val *= gsl_sf_gamma(a - (0.5 * (i - 1)));
		}
	val *=  pow(M_PI, (d * (d - 1) / 4.0));
	return(val);
}

//// Get Inverse Wishart PDF ////
double iwishpdf(gsl_matrix *X, gsl_matrix *Scale, gsl_matrix *inv, double dof){
	// From http://en.wikipedia.org/wiki/Inverse-Wishart_distribution
	double X_det, scale_det, denom, pdf, trace, numerator;
	int m = X->size1;
	int n = X->size2;

	// Allocate matrix for inv(X)
	gsl_matrix *for_mult = gsl_matrix_alloc(m, n);

	// Get determinant of X and Scale matrix
	X_det = matrix_determ(X);
	scale_det = matrix_determ(Scale);

	// Invert X
	inv_matrix(X,inv);

	// Multiple Scale * inv(X)
	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, Scale, inv, 1.0, for_mult);

	// Get trace of above.
	trace = matrix_trace(for_mult);

	numerator =  pow(scale_det, dof / 2.0) * pow(X_det, (-dof-m-1)/ 2.0) * exp(-0.5 * trace);
	denom = pow(2,dof * m / 2) * mv_gamma(dof/2, m);

	pdf = (numerator/denom);

	return(pdf);
}

//// Invert a matrix and keep input matrix untouched ////
gsl_matrix inv_matrix(gsl_matrix *X, gsl_matrix *inv){
	int s;
	int n = X->size1;
	int m = X->size2;
	gsl_matrix *a_copy = gsl_matrix_alloc(m, n);
	gsl_matrix_memcpy( a_copy, X );
	gsl_permutation *P = gsl_permutation_alloc(n);
	gsl_linalg_LU_decomp(a_copy, P, &s);
	gsl_linalg_LU_invert(a_copy, P, inv);
	return(*inv);
}

//// Print a matrix to stdout ////
void print_matrix(gsl_matrix *X){
	int i, j;
	int n = X->size1;
	int m = X->size2;
	printf("\n");
	for(i=0; i<n; i++){
		for(j=0; j<m; j++){
			printf("%05.2f ", gsl_matrix_get(X, i, j));
		}
		printf("\n");
	}
}

///////// Added Sat 01/19/2013 15:24:41 EST /////////

gsl_vector array_to_gsl_vec(gsl_vector *dest, double src[3]){
	int i;
	int n = dest->size;
		for (i = 0; i < n; ++i){

				gsl_vector_set(dest, i, src[i]);
			}
	return(*dest);
}

///// Transform user supplied matrix into a gsl matrix ////
gsl_matrix mdarray_to_gsl_mat(gsl_matrix *dest, double src[3][3]){
	int i, j;
	int n = dest->size1;
	int m = dest->size2;
		for (i = 0; i < m; ++i){
			for (j = 0; j < n; ++j){
				gsl_matrix_set(dest, i, j, src[i][j]);
			}
		}
	return(*dest);
}

//// Print a vector to stdout ////
void print_vector(gsl_vector *src){
	int i;
	int n = src->size;
	printf("\n");
	for (i = 0; i < n; ++i){
		printf("%05.4f ", gsl_vector_get(src,i));
	}
	printf("\n");
}

//// Do cholesky decomp to matrix ////
gsl_matrix get_chol(gsl_matrix *src, gsl_matrix *dest){
	gsl_matrix_memcpy(dest, src);
	gsl_linalg_cholesky_decomp(dest);
	return(*dest);
}

//// Add gaussian noise to chol decomp matrix ////
gsl_matrix mat_noise(gsl_matrix *src, gsl_rng *r_num){
	int i,j;
	int n = src->size1;
	int m = src->size2;
	double tmp;

	for(i = 0; i < n; i++){
		for(j = 0; j < m; j++){
			if(i >= j){
				tmp = gsl_matrix_get(src,i,j) + gsl_ran_gaussian(r_num,0.05);
				gsl_matrix_set(src, i, j, tmp);
			}

		}
	}
	return(*src);
}

//// Transform cholesky decomp matrix into vector ////
void fill_chol_vec(gsl_matrix *src, gsl_vector *dest){
	int i, j, count = 0;

	for(i = 0; i < src->size1; i++){
		for(j=0;j<src->size2;j++){
			if(i<=j){
				 gsl_vector_set(dest,count++, gsl_matrix_get(src,i,j));
			}

		}
	}
}

//// Add gaussian noise to each element in vector ////
gsl_vector vec_noise(gsl_rng *r_num, gsl_vector *src, gsl_vector *dest){
	int i;
	int n = src->size;
	double tmp;

	for(i=0;i<n;i++){
			tmp = gsl_vector_get(src,i) + gsl_ran_gaussian(r_num, 0.05);
			gsl_vector_set(dest,i,tmp);
		}
		return(*dest);
}

//// Multiply cholesky decomp matrix by its transpose ////
gsl_matrix chol_mult(gsl_matrix *src, gsl_matrix *dest){
	gsl_blas_dgemm(CblasNoTrans, CblasTrans,1.0, src, src,1.0, dest);
	return(*dest);
}

//// Zeros out everything but lower tri-- in place ////
//// From adk's adkGSL.c ////
gsl_matrix gsl_matrix_lower_tri(gsl_matrix *src){
	int i, j;
	for(i = 0; i < src->size1; i++){
		for(j = 0; j < src->size2; j++){
			if(i<j)	gsl_matrix_set(src,i,j,0.0);
		}
	}
	return(*src);
}

// Random Wishart matrix.  Taken from Andy Kern's mvnorm_gsl.c
/*void rwishartGelman(gsl_matrix *Scale, double dof){
*//* Wishart distribution random number generator */
/*
*	n	 gives the dimension of the random matrix
*	dof	 degrees of freedom
*	scale	 scale matrix of dimension n x n
*	result	 output variable with a single random matrix Wishart distribution generation
*/
/*	int i;
	gsl_vector gsl_vector_alloc(ndim)

	gsl_vector_set_all(mean,0.0);
	gsl_matrix_set_all(output,0.0);
	for(i = 0; i < dof; i++){
		gsl_vector_set_all(xm,0);
		rmvnorm_prealloc(r, n, mean, scale, work, xm);
		gsl_vector_outer_product(xm,xm,work2);
		gsl_matrix_add(output,work2);
	}
}
*/

/*void rwishartGelman(const gsl_rng *r, const int n, const int dof, const gsl_matrix *scale,
	gsl_matrix *work,gsl_matrix *work2, gsl_vector *mean,gsl_vector *xm, gsl_matrix *output)

rwishartGelman(r, ndim, nsnps+ndim, p2->winv,
			p2->work,p2->work2, p2->means, p2->xm, p2->covMat);

*/

max_stars_repo_name: andrewkern/pgLib

max_issues_repo_name: andrewkern/pgLib

max_forks_repo_name: andrewkern/pgLib


text:
#include <cblas.h>

void openblas_mul_double(double *a, double *b, double *c, int m, int k, int n)
{
    cblas_dgemm(
        CblasRowMajor, CblasNoTrans, CblasNoTrans,
        m, n, k,
        1.,
        a, k,
        b, n,
        0.,
        c, n
    );
}


text:
/* Useful general-purpose functions.
 *
 * Krzysztof Chalupka, 2017.
 */

#include <stdio.h>
#include <gsl/gsl_blas.h>
#include "utils.h"

void crit_err(char *msg){
  fprintf(stderr, "%s, exiting.\n", msg);
  exit(1);
}

int get_nlines(char *fname){
  /* Check how many lines there are in a file. */
  int nlines = 0;
  FILE *f = fopen(fname, "r");
  while(!feof(f))
    if (fgetc(f) == '\n')
      nlines++;
  fclose(f);
  printf("%d lines read from %s.\n\n", nlines, fname);
  return nlines;
}

void print_intarray(int *arr, int n){
  int i;
  printf("[");
  for(i = 0; i < n; i++)
    printf("%d ", arr[i]);
  printf("]");
}

void print_gsl_vector(gsl_vector *vec){
  int i;
  for (i = 0; i < vec->size; i++)
    printf("%.4f ", vec->data[i]);
}

gsl_vector *vector_sub(gsl_vector *v1, gsl_vector *v2){
  /* Create a new vector equal to v1 - v2. */
  gsl_vector *res = gsl_vector_calloc(v1->size);
  gsl_vector_add(res, v1);
  gsl_vector_sub(res, v2);
  return res;
}

double vector_dist(gsl_vector *v1, gsl_vector *v2){
  /* Compute the Euclidean distance between vectors. */
  double dist;
  gsl_vector *tmp =  vector_sub(v1, v2);
  dist = gsl_blas_dnrm2(tmp);
  free(tmp);
  return dist;
}

gsl_vector **load_vectors_from_file(char *fname, int nvecs)
{
  int vec_id;
  FILE *f;
  gsl_vector **vecs;
  double vals[2] = {0., 0.};

  /* Allocate the vector array. */
  vecs = (gsl_vector **) calloc(nvecs, sizeof(gsl_vector *));

  /* Load the vectors from a file into an array. */
  f = fopen(fname, "r");
  vec_id = 0;
  while(fscanf(f,"%lf %lf", &vals[0], &vals[1]) != EOF){
    vecs[vec_id] = gsl_vector_alloc(2);
    gsl_vector_set(vecs[vec_id], 0, vals[0]);
    gsl_vector_set(vecs[vec_id], 1, vals[1]);
    vec_id++;
  }

  /* Clean up. */
  fclose(f);
  printf("Data: \n");
  for (vec_id = 0; vec_id < nvecs; vec_id++)
    printf("[%f, %f]\n",
	   vecs[vec_id]->data[0],
	   vecs[vec_id]->data[1]);
  printf("\n");
  return vecs;
}


text:
/*
# Program to run a Monte Carlo radiation transfer through the 2D
# simulations of GRB jets.
#
# Python code written by D. Lazzati at Oregonstate, C code written by Tyler Parsotan @ Oregon State
# ver 0.1 July 8, 2015
# ver 1.1 July 20, 2015: added record of number of scatterings, included
# 	all terms in weight. Should now give correct light curves.
# ver 1.2 July 21, 2015: added parameter file to keep track of input
# 	params of each simulation

# ver 2.0 July 22, 2015: corrected the problem that arises when there is
# 	no scattering in the time span of one frame. Fixed output arrays dimension.

# ver 2.1 July 25, 2015: fixed bug that did not make the number of
# 	scattering grow with the number of photons.

# ver 3.0 July 28, 2015: using scipy nearest neighbor interpolation to
# 	speed things up. Gained about factor 2

# ver 3.1 July 29, 2015: added radial spread of photon injection points
# ver 3.2 July 31, 2015: added Gamma to the weight of photons!!!

# ver 4.0 Aug 5, 2015: try to speed up by inverting cycle
# ver 4.1 Aug 8, 2015: add spherical test as an option
# ver 4.2 Aug 9, 2015: saving files appending rather than re-writing
# ver 4.3 Aug 11, 2015: corrected error in the calculation of the local temperature
# ver 4.4 Aug 13, 2015: added cylindrical test
# ver 4.5 Aug 18, 2015: fixd various problems pointed by the cylindrical test
# ver 4.6 Aug 21, 2015: corrected mean free path for large radii

# ver 5.0 Aug 25, 2015: corrected problem with high-T electrons and excess scatterings
# ver 5.1 Aug 25, 2015: cleaned-up coding
# ver 5.2 Sept 3, 2015: fixed problem with number of scatterings for multiple injections
 *
 * ver 6.0 Dec 28, 2016: rewrote the code in C, added checkpoint file so if the code is interrupted all the progress wont be lost, made the code only need to be compiled once for a given MC_XXX directory path
                                            so you just need to supply the sub directory of MC_XXX as a command line argument
*/

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <unistd.h>
#include <dirent.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include "mclib.h"

#define THISRUN "Spherical"
#define FILEPATH "/Volumes/parsotat/Documents/16TI/"
#define FILEROOT "rhd_jet_big_13_hdf5_plt_cnt_"
#define MC_PATH "CMC_16TI_SPHERICAL/"
#define MCPAR "mc.par"

int main(int argc, char **argv)
{
    //compile each time a macro is changed, have to supply the subfolder within the MC_PATH directory as a command line argument to the C program eg. MCRAT 1/

	// Define variables
	char flash_prefix[200]="";
	char flash_file[200]="";
	char mc_dir[200]="" ;
	char mc_file[200]="" ;
    char mc_filename[200]="";
    char mc_operation[200]="";
    char this_run[200]=THISRUN;
    char *cyl="Cylindrical";
    char *sph="Spherical";
    int file_count = 0;
    DIR * dirp;
    struct dirent * entry;
	char spect;//type of spectrum
    char restrt;//restart or not
	double fps, theta_jmin, theta_jmax ;//frames per second of sim, min opening angle of jet, max opening angle of jet in radians
	double inj_radius, ph_weight ;//radius at chich photons are injected into sim
	int frm0,last_frm, frm2, photon_num ;//frame starting from, last frame of sim, frame of last injection, number of photons (not necessary any more)
	double *xPtr=NULL,  *yPtr=NULL,  *rPtr=NULL,  *thetaPtr=NULL,  *velxPtr=NULL,  *velyPtr=NULL,  *densPtr=NULL,  *presPtr=NULL,  *gammaPtr=NULL,  *dens_labPtr=NULL;
    double *szxPtr=NULL,*szyPtr=NULL, *tempPtr=NULL; //pointers to hold data from FLASH files
    int num_ph=0, array_num=0, ph_scatt_index=0, max_scatt=0, min_scatt=0,i=0; //number of photons produced in injection algorithm, number of array elleemnts from reading FLASH file, index of photon whch does scattering, generic counter
    double dt_max=0, thescatt=0, accum_time=0;
    double  gamma_infinity=0;
    double time_now=0, time_step=0, avg_scatt=0;
    double ph_dens_labPtr=0, ph_vxPtr=0, ph_vyPtr=0, ph_tempPtr=0;// *ph_cosanglePtr=NULL ;
    int frame=0, scatt_frame=0, frame_scatt_cnt=0, scatt_framestart=0, framestart=0;
    struct photon *phPtr=NULL; //pointer to array of photons

    long seed;
    const gsl_rng_type * T;
    gsl_rng * rand;

    gsl_rng_env_setup();
    T = gsl_rng_ranlxs0;
    rand = gsl_rng_alloc (T);

    //make strings of proper directories etc.
	snprintf(flash_prefix,sizeof(flash_prefix),"%s%s",FILEPATH,FILEROOT );
	snprintf(mc_dir,sizeof(flash_prefix),"%s%s%s",FILEPATH,MC_PATH, argv[1]);
	snprintf(mc_file,sizeof(flash_prefix),"%s%s",mc_dir,MCPAR);

    printf(">> mc.py: I am working on path: %s \n",mc_dir);

    printf(">> mc.py:  Reading mc.par\n");

    readMcPar(mc_file, &fps, &theta_jmin, &theta_jmax,&inj_radius,&frm0,&last_frm,&frm2,&photon_num, &ph_weight, &spect, &restrt);

    if (restrt=='c')
    {
        printf(">> mc.py:  Reading checkpoint\n");
        readCheckpoint(mc_dir, &phPtr, &framestart, &scatt_framestart, &num_ph, &restrt, &time_now);
        for (i=0;i<num_ph;i++)
        {
            printf("%e,%e,%e, %e,%e,%e, %e, %e\n",(phPtr+i)->p0, (phPtr+i)->p1, (phPtr+i)->p2, (phPtr+i)->p3, (phPtr+i)->r0, (phPtr+i)->r1, (phPtr+i)->r2, (phPtr+i)->num_scatt );
        }

        if (restrt=='c')
        {
            printf("Starting from photons injected at frame: %d out of %d\n", framestart, frm2);
            printf("Continuing scattering %d photons from frame: %d\n", num_ph, scatt_framestart);
            printf("The time now is: %e\n", time_now);
        }
        else
        {
            printf("Continuing simulation by injecting photons at frame: %d out of %d\n", framestart, frm2); //starting with new photon injection is same as restarting sim
        }
    }
    else
    {
        //remove everything from MC directory to ensure no corruption of data if theres other files there besides the mc.par file
        //for a checkpoint implementation, need to find the latest checkpoint file and read it and not delete the files
        printf(">> mc.py:  Cleaning directory\n");
        dirp = opendir(mc_dir);
        while ((entry = readdir(dirp)) != NULL)
        {
            if (entry->d_type == DT_REG) { /* If the entry is a regular file */
                 file_count++; //count how many files are in dorectory
            }
        }
        //printf("File count %d\n", file_count);

        if (file_count>2)
        {
            for (i=0;i<=last_frm;i++)
            {
                snprintf(mc_filename,sizeof(flash_prefix),"%s%s%d%s", mc_dir,"mcdata_",i,"_P0.dat");
                if( access( mc_filename, F_OK ) != -1 )
                {
                    snprintf(mc_operation,sizeof(flash_prefix),"%s%s%s%d%s","exec rm ", mc_dir,"mcdata_",i,"_*.dat"); //prepares string to remove *.dat in mc_dir
                    //printf("%s",mc_operation);
                    system(mc_operation);
                }
            }
        }
        framestart=frm0; //if restarting then start from parameters given in mc.par file
        scatt_framestart=framestart;
    }

    dt_max=1.0/fps;

    //loop over frames
    //for a checkpoint implementation, start from the last saved "frame" value and go to the saved "frm2" value
    for (frame=framestart;frame<=frm2;frame++)
    {
         if (restrt=='r')
         {
            time_now=frame/fps; //for a checkpoint implmentation, load the saved "time_now" value when reading the ckeckpoint file otherwise calculate it normally
         }

        printf(">> mc.py: Working on Frame %d\n", frame);

        if (restrt=='r')
        {
            //put proper number at the end of the flash file
            if (frame<10)
            {
                snprintf(flash_file,sizeof(flash_prefix), "%s%.3d%d",flash_prefix,000,frame);
            }
            else if (frame<100)
            {
                snprintf(flash_file,sizeof(flash_prefix), "%s%.2d%d",flash_prefix,00,frame);
            }
            else if (frame<1000)
            {
                snprintf(flash_file,sizeof(flash_prefix), "%s%d%d",flash_prefix,0,frame);
            }
            else
            {
                snprintf(flash_file,sizeof(flash_prefix), "%s%d",flash_prefix,frame);
            }

            printf(">> mc.py: Opening FLASH file %s\n",flash_file);
            //read in FLASH file
            //for a checkpoint implmentation, dont need to read the file yet
            readAndDecimate(flash_file, inj_radius, &xPtr,  &yPtr,  &szxPtr, &szyPtr, &rPtr,\
                &thetaPtr, &velxPtr,  &velyPtr,  &densPtr,  &presPtr,  &gammaPtr,  &dens_labPtr, &tempPtr, &array_num);

            //check for run type
            if(strcmp(cyl, this_run)==0)
            {
                //printf("In cylindrical prep\n");
                cylindricalPrep(gammaPtr, velxPtr, velyPtr, densPtr, dens_labPtr, presPtr, tempPtr, array_num);
            }
            else if (strcmp(sph, this_run)==0)
            {
                sphericalPrep(rPtr, xPtr, yPtr,gammaPtr, velxPtr, velyPtr, densPtr, dens_labPtr, presPtr, tempPtr, array_num );
            }

            //determine where to place photons and how many should go in a given place
            //for a checkpoint implmentation, dont need to inject photons, need to load photons' last saved data
            printf(">> mc.py: Injecting photons\n");
            photonInjection(&phPtr, &num_ph, inj_radius, ph_weight, spect, array_num, fps, theta_jmin, theta_jmax, xPtr, yPtr, szxPtr, szyPtr,rPtr,thetaPtr, tempPtr, velxPtr, velyPtr,rand );
            printf("%d\n",num_ph); //num_ph is one more photon than i actually have
            /*
            for (i=0;i<num_ph;i++)
                printf("%e,%e,%e \n",(phPtr+i)->r0, (phPtr+i)->r1, (phPtr+i)->r2 );
            */
        }

        //scatter photons all the way thoughout the jet
        //for a checkpoint implmentation, start from the last saved "scatt_frame" value eh start_frame=frame or start_frame=cont_frame
        if (restrt=='r')
        {
            scatt_framestart=frame; //have to make sure that once the inner loop is done and the outer loop is incrememnted by one the inner loop starts at that new value and not the one read by readCheckpoint()
        }

        for (scatt_frame=scatt_framestart;scatt_frame<=last_frm;scatt_frame++)
        {
            printf(">>\n");
            printf(">> mc.py: Working on photons injected at frame: %d out of %d\n", frame, frm2);
            printf(">> mc.py: %s: Working on frame %d\n", THISRUN, scatt_frame);
            printf(">> mc.py: Opening file...\n");

            //put proper number at the end of the flash file
            if (scatt_frame<10)
            {
                snprintf(flash_file,sizeof(flash_prefix), "%s%.3d%d",flash_prefix,000,scatt_frame);
            }
            else if (scatt_frame<100)
            {
                snprintf(flash_file,sizeof(flash_prefix), "%s%.2d%d",flash_prefix,00,scatt_frame);
            }
            else if (scatt_frame<1000)
            {
                snprintf(flash_file,sizeof(flash_prefix), "%s%d%d",flash_prefix,0,scatt_frame);
            }
            else
            {
                snprintf(flash_file,sizeof(flash_prefix), "%s%d",flash_prefix,scatt_frame);
            }

            readAndDecimate(flash_file, inj_radius, &xPtr,  &yPtr,  &szxPtr, &szyPtr, &rPtr,\
                &thetaPtr, &velxPtr,  &velyPtr,  &densPtr,  &presPtr,  &gammaPtr,  &dens_labPtr, &tempPtr, &array_num);
                //printf("The result of read and decimate are arrays with %d elements\n", array_num);

            //check for run type
            if(strcmp(cyl, this_run)==0)
            {
                printf("In cylindrical prep\n");
                cylindricalPrep(gammaPtr, velxPtr, velyPtr, densPtr, dens_labPtr, presPtr, tempPtr, array_num);
            }
            else if (strcmp(sph, this_run)==0)
            {
                sphericalPrep(rPtr, xPtr, yPtr,gammaPtr, velxPtr, velyPtr, densPtr, dens_labPtr, presPtr, tempPtr, array_num );
            }

            printf(">> mc.py: propagating and scattering %d photons\n", num_ph);

            frame_scatt_cnt=0;
            while (time_now<((scatt_frame+1)/fps))
            {
                //if simulation time is less than the simulation time of the next frame, keep scattering in this frame
                //go through each photon and find blocks closest to each photon and properties of those blocks to calulate mean free path
                //and choose the photon with the smallest mfp and calculate the timestep

                ph_scatt_index=findNearestPropertiesAndMinMFP(phPtr, num_ph, array_num, &time_step, xPtr,  yPtr, velxPtr,  velyPtr,  dens_labPtr, tempPtr,\
                    &ph_dens_labPtr, &ph_vxPtr, &ph_vyPtr, &ph_tempPtr, rand);

                //printf("In main: %e, %d, %e, %e\n", *(ph_num_scatt+ph_scatt_index), ph_scatt_index, time_step, time_now);
                //printf("In main: %e, %d, %e, %e\n",((phPtr+ph_scatt_index)->num_scatt), ph_scatt_index, time_step, time_now);

                 if (time_step<dt_max)
                {

                    //update number of scatterings and time
                    //(*(ph_num_scatt+ph_scatt_index))+=1;
                    ((phPtr+ph_scatt_index)->num_scatt)+=1;
                    frame_scatt_cnt+=1;
                    time_now+=time_step;

                    updatePhotonPosition(phPtr, num_ph, time_step);

                    //scatter the photon
                    //printf("Passed Parameters: %e, %e, %e\n", (ph_vxPtr), (ph_vyPtr), (ph_tempPtr));

                    photonScatter( (phPtr+ph_scatt_index), (ph_vxPtr), (ph_vyPtr), (ph_tempPtr), rand );

                    if (frame_scatt_cnt%1000 == 0)
                    {
                        printf("Scattering Number: %d\n", frame_scatt_cnt);
                        printf("The local temp is: %e, Local energy: %e\n", (ph_tempPtr),  (ph_tempPtr)*(1.380658e-16/1.6e-9));
                        printf("Average photon energy is: %e keV\n", averagePhotonEnergy(phPtr, num_ph)/1.6e-9); //write function to average over the photons p0*c and then do (/1.6e-9)
                    }

                }
                else
                {
                    time_now+=dt_max;

                    //for each photon update its position based on its momentum

                    updatePhotonPosition(phPtr, num_ph, dt_max);
                }

                //printf("In main 2: %e, %d, %e, %e\n", ((phPtr+ph_scatt_index)->num_scatt), ph_scatt_index, time_step, time_now);

            }

        //get scattering statistics
        phScattStats(phPtr, num_ph, &max_scatt, &min_scatt, &avg_scatt);

        printf("The number of scatterings in this frame is: %d\n", frame_scatt_cnt);
        printf("The last time step was: %lf.\nThe time now is: %lf\n", time_step,time_now);
        printf("The maximum number of scatterings for a photon is: %d\nThe minimum number of scattering for a photon is: %d\n", max_scatt, min_scatt);
        printf("The average number of scatterings thus far is: %lf\n", avg_scatt);

        printPhotons(phPtr, num_ph,  scatt_frame , mc_dir);

        //for a checkpoint implmentation,save the checkpoint file here after every 5 frames or something
        //save the photons data, the scattering number data, the scatt_frame value, and the frame value
        //WHAT IF THE PROGRAM STOPS AFTER THE LAST SCATT_FRAME, DURING THE FIRST SCATT_FRAME OF NEW FRAME VARIABLE - save restrt variable as 'r'
        printf(">> mc.py: Making checkpoint file\n");
        saveCheckpoint(mc_dir, frame, scatt_frame, num_ph, time_now, phPtr, last_frm);

        free(xPtr);free(yPtr);free(szxPtr);free(szyPtr);free(rPtr);free(thetaPtr);free(velxPtr);free(velyPtr);free(densPtr);free(presPtr);
        free(gammaPtr);free(dens_labPtr);free(tempPtr);
        xPtr=NULL; yPtr=NULL;  rPtr=NULL;thetaPtr=NULL;velxPtr=NULL;velyPtr=NULL;densPtr=NULL;presPtr=NULL;gammaPtr=NULL;dens_labPtr=NULL;
        szxPtr=NULL; szyPtr=NULL; tempPtr=NULL;
        }
        restrt='r';//set this to make sure that the next iteration of propogating photons doesnt use the values from the last reading of the checkpoint file
        free(phPtr);
        phPtr=NULL;

    }

    gsl_rng_free (rand);
	return 0;

}

max_stars_repo_name: outflows/MCRaT

max_issues_repo_name: outflows/MCRaT

max_forks_repo_name: outflows/MCRaT


text:
/*
 * Copyright (c) 2016-2021 lymastee, All rights reserved.
 * Contact: lymastee@hotmail.com
 *
 * This file is part of the gslib project.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */

#pragma once

#ifndef error_e0b666c0_4caa_4710_bb62_0f8a9d0e579e_h
#define error_e0b666c0_4caa_4710_bb62_0f8a9d0e579e_h

#include <gslib/type.h>

__gslib_begin__

struct errinfo
{
    gchar*      desc;
    gchar*      file;
    int         line;
    void*       user;
};

typedef void* (*error_dump)(void*);
gs_export extern void* error_dump_callstack(void*);

gs_export extern void _set_error(const gchar* desc, error_dump dump, void* user, const gchar* file, int line);
gs_export extern void _set_last_error(const gchar* desc, error_dump dump, void* user, const gchar* file, int line);
gs_export extern void pop_error();
gs_export extern errinfo* get_last_error();
gs_export extern void reset_error();

#ifndef swsc_char
#define swsc_char(x) _t(x)
#endif

#ifdef set_error
#undef set_error
#endif
#define set_error(desc, ...) do { \
    assert(!desc); \
    string str; \
    str.format(desc, __VA_ARGS__); \
    _set_error(str.c_str(), 0, 0, swsc_char(__FILE__), __LINE__); \
} while(0)

#ifdef set_last_error
#undef set_last_error
#endif
#define set_last_error(desc, ...) do { \
    assert(!desc); \
    string str; \
    str.format(desc, _VA_ARGS__); \
    _set_last_error(str.c_str(), 0, 0, swsc_char(__FILE__), __LINE__); \
} while(0)

//gs_export extern void dumperr_file(void);
gs_export extern void trace_hold(const gchar* fmt, ...);
gs_export extern void _trace_to_clipboard();

#ifdef _GS_TRACE_TO_CLIPBOARD
#define trace(fmt, ...) do { trace_hold(fmt, __VA_ARGS__); } while(0)
#define trace_to_clipboard _trace_to_clipboard
#elif defined (DEBUG) || defined (_DEBUG)
gs_export extern void trace(const gchar* fmt, ...);
gs_export extern void trace_all(const gchar* str);
#define trace_to_clipboard __noop
#else
#define trace(fmt, ...) __noop
#define trace_all __noop;
#define trace_to_clipboard __noop
#endif

/* beep an alarm */
gs_export extern void sound_alarm();

__gslib_end__

#endif

ext: h

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


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

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

void line_convolution()
{

  FILE *fp;
  int info, i, j, nq;
  double flux, err, tline, taup, tcon, fcon, fcon_err, dtau;
  double *yrec_line;
//  double sigma, taud;
  double *Larr, *ybuf, *Cq, *ICq, *ave, *yave;

  gsl_interp_accel *gsl_acc_conv, *gsl_acc_error_conv;
  gsl_interp  *gsl_linear_conv, *gsl_linear_error_conv;

  gsl_acc_conv = gsl_interp_accel_alloc();
  gsl_acc_error_conv = gsl_interp_accel_alloc();

  gsl_linear_conv = gsl_interp_alloc(gsl_interp_linear, ncon);
  gsl_linear_error_conv = gsl_interp_alloc(gsl_interp_linear, ncon);

  gsl_interp_init(gsl_linear_conv, Tcon, Fcon, ncon);
  gsl_interp_init(gsl_linear_error_conv, Tcon, Fcerrs, ncon);

  yrec_line = array_malloc(nline_data);

// ave
//  taud = exp(theta_best[1]);
//  sigma = exp(theta_best[0]) * sqrt(taud/2.0);

  nq = 2*(flag_detrend + 1);
  Larr = workspace;
  ybuf = Larr + nq*nall_data;
  Cq = ybuf + nall_data;
  ICq = Cq + nq*nq;
  ave = ICq + nq*nq;
  yave = ave + nall_data;

  set_covar_mat(theta_best);

  for(i=0; i<nall_data*nall_data; i++)
  {
    Cmat[i] = Smat[i] + Nmat[i];
  }

  for(i=0; i<ncon_data; i++)
  {
    if(flag_detrend==0)
    {
      Larr[i*nq+0] = 1.0;
      Larr[i*nq+1] = 0.0;
    }
    else
    {
      Larr[i*nq + 0] = 1.0;
      Larr[i*nq + 1] = Tcon_data[i];
      Larr[i*nq + 2] = 0.0;
      Larr[i*nq + 3] = 0.0;
    }
  }
  for(i=0; i<nline_data; i++)
  {
    if(flag_detrend==0)
    {
      Larr[(i+ncon_data)*nq+0] = 0.0;
      Larr[(i+ncon_data)*nq+1] = 1.0;
    }
    else
    {
      Larr[(i+ncon_data)*nq + 0] = 0.0;
      Larr[(i+ncon_data)*nq + 1] = 0.0;
      Larr[(i+ncon_data)*nq + 2] = 1.0;
      Larr[(i+ncon_data)*nq + 3] = Tline_data[i];
    }
  }

  memcpy(ICmat, Cmat, nall_data*nall_data*sizeof(double));
  inverse_mat(ICmat, nall_data, &info);

  multiply_mat_MN(ICmat, Larr, Tmat1, nall_data, nq, nall_data);
  multiply_mat_MN_transposeA(Larr, Tmat1, ICq, nq, nq, nall_data);

  memcpy(Cq, ICq, nq*nq*sizeof(double));
  inverse_mat(Cq, nq, &info);

  multiply_mat_MN_transposeA(Larr, ICmat, Tmat1, nq, nall_data, nall_data);
  multiply_mat_MN(Cq, Tmat1, Tmat2, nq, nall_data, nq);
  multiply_mat_MN(Tmat2, Fall_data, ave, nq, 1, nall_data);
  multiply_mat_MN(Larr, ave, yave, nall_data, 1, nq);
  //printf("%f %f\n", ave[0], ave[1]);

  dtau = TF_tau[1] - TF_tau[0];
  fp = fopen("data/line_conv.txt", "w");

  for(i=0; i<nline_data; i++)
  {
    flux = 0.0;
    err = 0.0;
    tline = Tline_data[i];
    for(j=0; j<ntau; j++)
    {
      taup = TF_tau[j];
      tcon = tline - taup;
      if(tcon >=Tcon[0] && tcon <= Tcon[ncon-1])
      {
        fcon = gsl_interp_eval(gsl_linear_conv, Tcon, Fcon, tcon, gsl_acc_conv);
        fcon_err = gsl_interp_eval(gsl_linear_error_conv, Tcon, Fcerrs, tcon, gsl_acc_error_conv);
        flux += fcon * TF[j];
        err += pow(fcon_err * TF[j] * dtau, 2);
      }
    }
    flux *= dtau;
    err = sqrt(err);
    yrec_line[i] = flux + yave[i+ncon_data];
    fprintf(fp, "%f  %e %e\n", tline, yrec_line[i]*scale_line, err*scale_line);
  }
  fclose(fp);

  gsl_interp_free(gsl_linear_conv);
  gsl_interp_free(gsl_linear_error_conv);
  gsl_interp_accel_free(gsl_acc_conv);
  gsl_interp_accel_free(gsl_acc_error_conv);
}

max_stars_repo_name: LiyrAstroph/MICA

max_issues_repo_name: LiyrAstroph/MICA

max_forks_repo_name: LiyrAstroph/MICA


text:
#ifndef TFDH_GSL_WRAPPERS_H
#define TFDH_GSL_WRAPPERS_H

#include <gsl/gsl_spline.h>
#include <vector>

namespace GSL {

  // simple functor class to help interfacing with gsl_function.
  // derived classes representing a function f(x) should hold any needed
  // parameters as member variables, and operator()(double x) should evaluate
  // the function f(x).
  //
  // NOTE: the GSL wrappers below could in principle take C++11 lambdas instead
  //       of this custom function object, however i did not find an elegant
  //       way to call these lambdas using GSL's C interfaces...
  class FunctionObject {
    public:
      virtual ~FunctionObject() = default;
      virtual double operator()(double x) const = 0;
  };

  // simple wrapper class around GSL splines
  class Spline {
    private:
      gsl_interp_accel* acc;
      gsl_spline* spline;
    public:
      Spline(const std::vector<double>& x, const std::vector<double>& f);
      ~Spline();
      double eval(double r) const;
  };

  // find a root of func in the interval x1 to x2
  // one of (x1,x2) or (x2,x1) MUST bracket a root
  // implemented using GSL's Brent rootfinder method
  double findRoot(const GSL::FunctionObject& func, double x1, double x2,
      double eps_abs, double eps_rel);

  // definite integral of func from x1 to x2
  // implemented using GSL's adaptive Gauss quadrature
  double integrate(const GSL::FunctionObject& func, double x1, double x2,
      double eps_abs, double eps_rel);

}

#endif // TFDH_GSL_WRAPPERS_H

max_stars_repo_name: fmahebert/tfdh

max_issues_repo_name: fmahebert/tfdh

max_forks_repo_name: fmahebert/tfdh


text:
#include <stdio.h>
#include <math.h>
#include <sys/types.h>
#include <time.h>
#include <complex.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_sf_gamma.h>

#define FALSE 0
#define TRUE 1

#ifndef M_PI
#define M_PI            3.141592653589793238462643383280        /* pi */
#endif

#define LN2PI 1.83787706640934548356065947

extern char parameterfilename[50];
extern char parameterlogfilename[50];

#define MBIG 1000000000
#define MSEED 161803398
#define MZ 0
#define FAC (1.0/MBIG)
#define MIN(A,B) ((A) < (B) ? (A) : (B))
#define MAX(A,B) ((A) > (B) ? (A) : (B))
enum {NOT_AVAILABLE=0, ABOVE_DETECTION=1, BELOW_DETECTION=2, MEASURED=3};

typedef struct {
  int n;
  int m;
  int **a;
} Matrix;

void ReadParam( char*, void*, char* ) ;
void ExperimentNumber( char*, char* );
void Error( char* );
void FileError( char*, char* );
double ***doublematrix3( int, int, int );
double **doublematrix( int, int );
int **integermatrix( int, int );
uint **uintegermatrix( int, int );
uint *uintegervector( int );
char **charmatrix( int, int );
double *doublevector( int );
int *integervector( int );
char *charvector( int );

void Int_HeapSort(int*, int);
void HeapSort( double*, int );
void circmeanvar( double*, int, double*, double* );
void meanvar( double*, int, double*, double* );
void meanvarl( long double*, int, long double*, long double* );
void wmeanvar( long double *x, int n, long double *m, long double *s2 );
double mean( double*, int, int );
int choldc_block(double**, long long);

uint   binomialRNG(uint, double, const gsl_rng*);
void   multinomialRNG(uint, const double*, size_t, unsigned int*,
		      const gsl_rng*);
int    poissonRNG(double, const gsl_rng*);
int    uniform_intRNG(int, int, const gsl_rng*);
double normalRNG(double, double, const gsl_rng*);
double betaRNG(double, double, const gsl_rng*);
double chisqRNG(double, const gsl_rng*);
double scaleinvchisqRNG(double, double, const gsl_rng*);
double invgammaRNG(double, double, const gsl_rng*);
double gammaRNG(double, double, const gsl_rng*);
double uniformRNG(double, double, const gsl_rng*);
double paretoRNG(double, double, const gsl_rng*);
void   dirichlet_multinomialRNG(uint, const double*, size_t, uint*,
				const gsl_rng*);
void   dirichletRNG(const double*, size_t, double*, const gsl_rng*);
uint   beta_binomialRNG(uint, double, double, const gsl_rng*);
uint   negative_binomialRNG(double, double, const gsl_rng*);
uint   beta_negative_binomialRNG(uint, double, double, const gsl_rng*);
/* Matrix *wishartRNG(int, int, Matrix*, Matrix*, const gsl_rng*); */
/* void   MVNRNG(int, Matrix*, Matrix*, Matrix*, const gsl_rng*); */

double poissonCMF(const double, const double, int);
double gamma_scalePDF(const double, const double, const double);
double betaPDF(const double, const double, const double);
double exponentialPDF(const double, const double);
double poissonPMF(const double, const double);
double geometricPDF(const double, const double);
double trnormalPDF(const double, const double, const double);
double normalCDF(const double, const double, const double, int);
double normalPDF(const double, const double, const double);
double binomialPMF(const double, const double, const double);
double multinomialPMF(const unsigned int, const double*, const double*);
double uniformPDF(const double, const double, const double);
double logisticPDF(const double, const double, const double);
double betabinomialPMF(const double, const double, const double, const double);
/* double MVNPDF(Matrix*, Matrix*, Matrix*, const double, int); */

double identity(double);
double logit(double);
double invlogit(double);

Matrix *MatrixInit(Matrix*, int, int, int*);
Matrix *MatrixDestroy(Matrix*);
/* Matrix *Matrix_Init(Matrix*, char*, int, int); */
/* Matrix *Matrix_Fill(Matrix*, int, int, ...); */
/* Matrix *Matrix_Destroy(Matrix*); */
/* Matrix *Matrix_Copy(Matrix*, Matrix*); */
/* Matrix *Matrix_Sub_Copy(Matrix*, Matrix*, int); */
/* Matrix *Matrix_Scalar_Multiply(Matrix*, double, Matrix*); */
/* Matrix *Matrix_Multiply(Matrix*, Matrix*, Matrix*); */
/* Matrix *Matrix_Multiply(Matrix*, Matrix*, Matrix*); */
/* Matrix *Matrix_Cholesky(Matrix*, Matrix*); */
/* Matrix *Matrix_Invert(Matrix*, Matrix*); */
/* Matrix *Matrix_Add(Matrix*, Matrix*, Matrix*); */
/* Matrix *Matrix_Add_Diagonal(Matrix*, Matrix*); */
/* Matrix *Matrix_Subtract(Matrix*, Matrix*, Matrix*); */
/* Matrix *Matrix_Transpose(Matrix*, Matrix*); */
/* double Matrix_Log_Determinant(Matrix*, int); */
/* void   Matrix_Print(Matrix*); */

double rtnewt(void(*)(double, double*, double*, double*),
	      double, double, double, double*);
double rtbis( double (*)(double, double*), double, double, double, int*,
	      double* );

size_t strlen(const char*);
char *strncpy(char*, const char*, size_t);
char *strcpy(char*, const char*);
char *strcat(char*, const char*);
char *strchr ( const char*, int );
int  strcmp (const char*, const char*);
void *calloc(size_t, size_t);
void *malloc(size_t);
void free(void*);
void exit(int);

double von_mises_cdf ( double x, double a, double b );
double von_mises_cdf_inv ( double cdf, double a, double b );


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

/* These tests are based on the NIST Statistical Reference Datasets
   See http://www.nist.gov/itl/div898/strd/index.html for more
   information. */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_multifit_nlin.h>
#include <gsl/gsl_blas.h>

#include <gsl/gsl_ieee_utils.h>

#include "test_longley.c"
#include "test_filip.c"
#include "test_pontius.c"
#include "test_brown.c"
#include "test_enso.c"
#include "test_kirby2.c"
#include "test_hahn1.c"
#include "test_nelson.c"
#include "test_fn.c"
#include "test_estimator.c"

void
test_lmder (gsl_multifit_function_fdf * f, double x0[],
            double * X, double F[], double * cov);

void
test_fdf (const char * name, gsl_multifit_function_fdf * f,
          double x0[], double x[], double sumsq,
          double sigma[]);

int
main (void)
{
  gsl_ieee_env_setup();

  test_longley();
  test_filip();
  test_pontius();
  test_estimator();

  {
    gsl_multifit_function_fdf f = make_fdf (&brown_f, &brown_df, &brown_fdf,
                                            brown_N, brown_P, 0);

    test_lmder(&f, brown_x0, &brown_X[0][0], brown_F, &brown_cov[0][0]);
  }

  {
    gsl_multifit_function_fdf f = make_fdf (&enso_f, &enso_df, &enso_fdf,
                                            enso_N, enso_P, 0);

    test_fdf("nist-ENSO", &f, enso_x0, enso_x, enso_sumsq, enso_sigma);
  }

  {
    gsl_multifit_function_fdf f = make_fdf (&kirby2_f, &kirby2_df, &kirby2_fdf,
                                            kirby2_N, kirby2_P, 0);

    test_fdf("nist-kirby2", &f, kirby2_x0, kirby2_x, kirby2_sumsq, kirby2_sigma);
  }

  {
    gsl_multifit_function_fdf f = make_fdf (&hahn1_f, &hahn1_df, &hahn1_fdf,
                                            hahn1_N, hahn1_P, 0);

    test_fdf("nist-hahn1", &f, hahn1_x0, hahn1_x, hahn1_sumsq, hahn1_sigma);
  }

#ifdef JUNK
  {
    gsl_multifit_function_fdf f = make_fdf (&nelson_f, &nelson_df, &nelson_fdf,
                                            nelson_N, nelson_P, 0);

    test_fdf("nist-nelson", &f, nelson_x0, nelson_x, nelson_sumsq, nelson_sigma);
  }
#endif

  /* now summarize the results */

  exit (gsl_test_summary ());
}

void
test_lmder (gsl_multifit_function_fdf * f, double x0[],
            double * X, double F[], double * cov)
{
  const gsl_multifit_fdfsolver_type *T;
  gsl_multifit_fdfsolver *s;

  const size_t n = f->n;
  const size_t p = f->p;

  int status;
  size_t iter = 0, i;

  gsl_vector_view x = gsl_vector_view_array (x0, p);

  T = gsl_multifit_fdfsolver_lmsder;
  s = gsl_multifit_fdfsolver_alloc (T, n, p);
  gsl_multifit_fdfsolver_set (s, f, &x.vector);

  do
    {
      status = gsl_multifit_fdfsolver_iterate (s);

      for (i = 0 ; i < p; i++)
        {
          gsl_test_rel (gsl_vector_get (s->x, i), X[p*iter+i], 1e-5,
                        "lmsder, iter=%u, x%u", iter, i);
        }

      gsl_test_rel (gsl_blas_dnrm2 (s->f), F[iter], 1e-5,
                    "lmsder, iter=%u, f", iter);

      iter++;
    }
  while (iter < 20);

  {
    size_t i, j;
    gsl_matrix * covar = gsl_matrix_alloc (4, 4);
    gsl_multifit_covar (s->J, 0.0, covar);

    for (i = 0; i < 4; i++)
      {
        for (j = 0; j < 4; j++)
          {
            gsl_test_rel (gsl_matrix_get(covar,i,j), cov[i*p + j], 1e-7,
                          "gsl_multifit_covar cov(%d,%d)", i, j) ;
          }
      }

    gsl_matrix_free (covar);
  }

  gsl_multifit_fdfsolver_free (s);

}

void
test_fdf (const char * name, gsl_multifit_function_fdf * f,
          double x0[], double x_final[],
          double f_sumsq, double sigma[])
{
  const gsl_multifit_fdfsolver_type *T;
  gsl_multifit_fdfsolver *s;

  const size_t n = f->n;
  const size_t p = f->p;

  int status;
  size_t iter = 0;

  gsl_vector_view x = gsl_vector_view_array (x0, p);

  T = gsl_multifit_fdfsolver_lmsder;
  s = gsl_multifit_fdfsolver_alloc (T, n, p);
  gsl_multifit_fdfsolver_set (s, f, &x.vector);

  do
    {
      status = gsl_multifit_fdfsolver_iterate (s);

#ifdef DEBUG
       printf("iter = %d  status = %d  |f| = %.18e x = \n",
         iter, status, gsl_blas_dnrm2 (s->f));

         gsl_vector_fprintf(stdout, s->x, "%.8e");
#endif
      status = gsl_multifit_test_delta (s->dx, s->x, 0.0, 1e-7);

      iter++;
    }
  while (status == GSL_CONTINUE && iter < 1000);

  {
    size_t i;
    gsl_matrix * covar = gsl_matrix_alloc (p, p);
    gsl_multifit_covar (s->J, 0.0, covar);

    for (i = 0 ; i < p; i++)
      {
        gsl_test_rel (gsl_vector_get (s->x, i), x_final[i], 1e-5,
                      "%s, lmsder, x%u", name, i);
      }

    {
      double s2 = pow(gsl_blas_dnrm2 (s->f), 2.0);

      gsl_test_rel (s2, f_sumsq, 1e-5, "%s, lmsder, |f|^2", name);

      for (i = 0; i < p; i++)
        {
          double ei = sqrt(s2/(n-p))*sqrt(gsl_matrix_get(covar,i,i));
          gsl_test_rel (ei, sigma[i], 1e-4,
                        "%s, sigma(%d)", name, i) ;
        }
    }

    gsl_matrix_free (covar);
  }

  /* Check that there is no hidden state, restarting should
     produce identical results. */

  {
    int status0, status1;
    size_t i;
    gsl_multifit_fdfsolver *t = gsl_multifit_fdfsolver_alloc (T, n, p);
    gsl_multifit_fdfsolver_set (t, f, &x.vector);

    /* do a few extra iterations to stir things up */

    gsl_multifit_fdfsolver_set (s, f, &x.vector);

    for (i = 0; i < 3; i++)
      {
        gsl_multifit_fdfsolver_iterate (s);
      }

    gsl_multifit_fdfsolver_set (s, f, &x.vector);

    do
      {
        status0 = gsl_multifit_fdfsolver_iterate (s);
        status1 = gsl_multifit_fdfsolver_iterate (t);

        gsl_test_int(status0, status1, "%s, lmsder status after set iter=%u", name, iter);

        for (i = 0; i < p; i++) {
          double sxi = gsl_vector_get(s->x,i);
          double txi = gsl_vector_get(t->x,i);
#ifdef DEBUG
          printf("%d %g %g\n", i, sxi, txi);
#endif
          gsl_test_rel(sxi, txi, 1e-15, "%s, lmsder after set, %u/%u", name, iter, i);
        }

#ifdef DEBUG
        printf("iter = %d  status = %d  |f| = %.18e x = \n",
               iter, status, gsl_blas_dnrm2 (s->f));

        gsl_vector_fprintf(stdout, s->x, "%.8e");
#endif
        status0 = gsl_multifit_test_delta (s->dx, s->x, 0.0, 1e-7);
        status1 = gsl_multifit_test_delta (t->dx, s->x, 0.0, 1e-7);

        gsl_test_int(status0, status1, "%s, lmsder test delta status after set iter=%u", name, iter);

        iter++;
      }
    while (status1 == GSL_CONTINUE && iter < 1000);

    gsl_multifit_fdfsolver_free (t);
  }

  gsl_multifit_fdfsolver_free (s);
}


text:
/* specfunc/legendre_H3d.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 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.
 */

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_exp.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_trig.h>
#include <gsl/gsl_sf_legendre.h>

#include "error.h"

#include "legendre.h"

/* See [Abbott+Schaefer, Ap.J. 308, 546 (1986)] for
 * enough details to follow what is happening here.
 */

/* Logarithm of normalization factor, Log[N(ell,lambda)].
 * N(ell,lambda) = Product[ lambda^2 + n^2, {n,0,ell} ]
 *               = |Gamma(ell + 1 + I lambda)|^2  lambda sinh(Pi lambda) / Pi
 * Assumes ell >= 0.
 */
static
int
legendre_H3d_lnnorm(const int ell, const double lambda, double * result)
{
  double abs_lam = fabs(lambda);

  if(abs_lam == 0.0) {
    *result = 0.0;
    GSL_ERROR ("error", GSL_EDOM);
  }
  else if(lambda > (ell + 1.0)/GSL_ROOT3_DBL_EPSILON) {
    /* There is a cancellation between the sinh(Pi lambda)
     * term and the log(gamma(ell + 1 + i lambda) in the
     * result below, so we show some care and save some digits.
     * Note that the above guarantees that lambda is large,
     * since ell >= 0. We use Stirling and a simple expansion
     * of sinh.
     */
    double rat = (ell+1.0)/lambda;
    double ln_lam2ell2  = 2.0*log(lambda) + log(1.0 + rat*rat);
    double lg_corrected = -2.0*(ell+1.0) + M_LNPI + (ell+0.5)*ln_lam2ell2 + 1.0/(288.0*lambda*lambda);
    double angle_terms  = lambda * 2.0 * rat * (1.0 - rat*rat/3.0);
    *result = log(abs_lam) + lg_corrected + angle_terms - M_LNPI;
    return GSL_SUCCESS;
  }
  else {
    gsl_sf_result lg_r;
    gsl_sf_result lg_theta;
    gsl_sf_result ln_sinh;
    gsl_sf_lngamma_complex_e(ell+1.0, lambda, &lg_r, &lg_theta);
    gsl_sf_lnsinh_e(M_PI * abs_lam, &ln_sinh);
    *result = log(abs_lam) + ln_sinh.val + 2.0*lg_r.val - M_LNPI;
    return GSL_SUCCESS;
  }
}

/* Calculate series for small eta*lambda.
 * Assumes eta > 0, lambda != 0.
 *
 * This is just the defining hypergeometric for the Legendre function.
 *
 * P^{mu}_{-1/2 + I lam}(z) = 1/Gamma(l+3/2) ((z+1)/(z-1)^(mu/2)
 *                            2F1(1/2 - I lam, 1/2 + I lam; l+3/2; (1-z)/2)
 * We use
 *       z = cosh(eta)
 * (z-1)/2 = sinh^2(eta/2)
 *
 * And recall
 * H3d = sqrt(Pi Norm /(2 lam^2 sinh(eta))) P^{-l-1/2}_{-1/2 + I lam}(cosh(eta))
 */
static
int
legendre_H3d_series(const int ell, const double lambda, const double eta,
                    gsl_sf_result * result)
{
  const int nmax = 5000;
  const double shheta = sinh(0.5*eta);
  const double ln_zp1 = M_LN2 + log(1.0 + shheta*shheta);
  const double ln_zm1 = M_LN2 + 2.0*log(shheta);
  const double zeta = -shheta*shheta;
  gsl_sf_result lg_lp32;
  double term = 1.0;
  double sum  = 1.0;
  double sum_err = 0.0;
  gsl_sf_result lnsheta;
  double lnN;
  double lnpre_val, lnpre_err, lnprepow;
  int stat_e;
  int n;

  gsl_sf_lngamma_e(ell + 3.0/2.0, &lg_lp32);
  gsl_sf_lnsinh_e(eta, &lnsheta);
  legendre_H3d_lnnorm(ell, lambda, &lnN);
  lnprepow = 0.5*(ell + 0.5) * (ln_zm1 - ln_zp1);
  lnpre_val  = lnprepow + 0.5*(lnN + M_LNPI - M_LN2 - lnsheta.val) - lg_lp32.val - log(fabs(lambda));
  lnpre_err  = lnsheta.err + lg_lp32.err + GSL_DBL_EPSILON * fabs(lnpre_val);
  lnpre_err += 2.0*GSL_DBL_EPSILON * (fabs(lnN) + M_LNPI + M_LN2);
  lnpre_err += 2.0*GSL_DBL_EPSILON * (0.5*(ell + 0.5) * (fabs(ln_zm1) + fabs(ln_zp1)));
  for(n=1; n<nmax; n++) {
    double aR = n - 0.5;
    term *= (aR*aR + lambda*lambda)*zeta/(ell + n + 0.5)/n;
    sum  += term;
    sum_err += 2.0*GSL_DBL_EPSILON*fabs(term);
    if(fabs(term/sum) < 2.0 * GSL_DBL_EPSILON) break;
  }

  stat_e = gsl_sf_exp_mult_err_e(lnpre_val, lnpre_err, sum, fabs(term)+sum_err, result);
  return GSL_ERROR_SELECT_2(stat_e, (n==nmax ? GSL_EMAXITER : GSL_SUCCESS));
}

/* Evaluate legendre_H3d(ell+1)/legendre_H3d(ell)
 * by continued fraction.
 */
#if 0
static
int
legendre_H3d_CF1(const int ell, const double lambda, const double coth_eta,
                 gsl_sf_result * result)
{
  const double RECUR_BIG = GSL_SQRT_DBL_MAX;
  const int maxiter = 5000;
  int n = 1;
  double Anm2 = 1.0;
  double Bnm2 = 0.0;
  double Anm1 = 0.0;
  double Bnm1 = 1.0;
  double a1 = hypot(lambda, ell+1.0);
  double b1 = (2.0*ell + 3.0) * coth_eta;
  double An = b1*Anm1 + a1*Anm2;
  double Bn = b1*Bnm1 + a1*Bnm2;
  double an, bn;
  double fn = An/Bn;

  while(n < maxiter) {
    double old_fn;
    double del;
    n++;
    Anm2 = Anm1;
    Bnm2 = Bnm1;
    Anm1 = An;
    Bnm1 = Bn;
    an = -(lambda*lambda + ((double)ell + n)*((double)ell + n));
    bn = (2.0*ell + 2.0*n + 1.0) * coth_eta;
    An = bn*Anm1 + an*Anm2;
    Bn = bn*Bnm1 + an*Bnm2;

    if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
      An /= RECUR_BIG;
      Bn /= RECUR_BIG;
      Anm1 /= RECUR_BIG;
      Bnm1 /= RECUR_BIG;
      Anm2 /= RECUR_BIG;
      Bnm2 /= RECUR_BIG;
    }

    old_fn = fn;
    fn = An/Bn;
    del = old_fn/fn;

    if(fabs(del - 1.0) < 4.0*GSL_DBL_EPSILON) break;
  }

  result->val = fn;
  result->err = 2.0 * GSL_DBL_EPSILON * (sqrt(n)+1.0) * fabs(fn);

  if(n >= maxiter)
    GSL_ERROR ("error", GSL_EMAXITER);
  else
    return GSL_SUCCESS;
}
#endif /* 0 */

/* Evaluate legendre_H3d(ell+1)/legendre_H3d(ell)
 * by continued fraction. Use the Gautschi (Euler)
 * equivalent series.
 */
 /* FIXME: Maybe we have to worry about this. The a_k are
  * not positive and there can be a blow-up. It happened
  * for J_nu once or twice. Then we should probably use
  * the method above.
  */
static
int
legendre_H3d_CF1_ser(const int ell, const double lambda, const double coth_eta,
                     gsl_sf_result * result)
{
  const double pre = hypot(lambda, ell+1.0)/((2.0*ell+3)*coth_eta);
  const int maxk = 20000;
  double tk   = 1.0;
  double sum  = 1.0;
  double rhok = 0.0;
  double sum_err = 0.0;
  int k;

  for(k=1; k<maxk; k++) {
    double tlk = (2.0*ell + 1.0 + 2.0*k);
    double l1k = (ell + 1.0 + k);
    double ak = -(lambda*lambda + l1k*l1k)/(tlk*(tlk+2.0)*coth_eta*coth_eta);
    rhok = -ak*(1.0 + rhok)/(1.0 + ak*(1.0 + rhok));
    tk  *= rhok;
    sum += tk;
    sum_err += 2.0 * GSL_DBL_EPSILON * k * fabs(tk);
    if(fabs(tk/sum) < GSL_DBL_EPSILON) break;
  }

  result->val  = pre * sum;
  result->err  = fabs(pre * tk);
  result->err += fabs(pre * sum_err);
  result->err += 4.0 * GSL_DBL_EPSILON * fabs(result->val);

  if(k >= maxk)
    GSL_ERROR ("error", GSL_EMAXITER);
  else
    return GSL_SUCCESS;
}

/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/

int
gsl_sf_legendre_H3d_0_e(const double lambda, const double eta, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(eta < 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(eta == 0.0 || lambda == 0.0) {
    result->val = 1.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else {
    const double lam_eta = lambda * eta;
    gsl_sf_result s;
    gsl_sf_sin_err_e(lam_eta, 2.0*GSL_DBL_EPSILON * fabs(lam_eta), &s);
    if(eta > -0.5*GSL_LOG_DBL_EPSILON) {
      double f = 2.0 / lambda * exp(-eta);
      result->val  = f * s.val;
      result->err  = fabs(f * s.val) * (fabs(eta) + 1.0) * GSL_DBL_EPSILON;
      result->err += fabs(f) * s.err;
      result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    }
    else {
      double f = 1.0/(lambda*sinh(eta));
      result->val  = f * s.val;
      result->err  = fabs(f * s.val) * (fabs(eta) + 1.0) * GSL_DBL_EPSILON;
      result->err += fabs(f) * s.err;
      result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    }
    return GSL_SUCCESS;
  }
}

int
gsl_sf_legendre_H3d_1_e(const double lambda, const double eta, gsl_sf_result * result)
{
  const double xi    = fabs(eta*lambda);
  const double lsq   = lambda*lambda;
  const double lsqp1 = lsq + 1.0;

  /* CHECK_POINTER(result) */

  if(eta < 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(eta == 0.0 || lambda == 0.0) {
    result->val = 0.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(xi < GSL_ROOT5_DBL_EPSILON && eta < GSL_ROOT5_DBL_EPSILON) {
    double etasq = eta*eta;
    double xisq  = xi*xi;
    double term1 = (etasq + xisq)/3.0;
    double term2 = -(2.0*etasq*etasq + 5.0*etasq*xisq + 3.0*xisq*xisq)/90.0;
    double sinh_term = 1.0 - eta*eta/6.0 * (1.0 - 7.0/60.0*eta*eta);
    double pre = sinh_term/sqrt(lsqp1) / eta;
    result->val  = pre * (term1 + term2);
    result->err  = pre * GSL_DBL_EPSILON * (fabs(term1) + fabs(term2));
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    double sin_term;     /*  Sin(xi)/xi     */
    double cos_term;     /*  Cos(xi)        */
    double coth_term;    /*  eta/Tanh(eta)  */
    double sinh_term;    /*  eta/Sinh(eta)  */
    double sin_term_err;
    double cos_term_err;
    double t1;
    double pre_val;
    double pre_err;
    double term1;
    double term2;
    if(xi < GSL_ROOT5_DBL_EPSILON) {
      sin_term = 1.0 - xi*xi/6.0 * (1.0 - xi*xi/20.0);
      cos_term = 1.0 - 0.5*xi*xi * (1.0 - xi*xi/12.0);
      sin_term_err = GSL_DBL_EPSILON;
      cos_term_err = GSL_DBL_EPSILON;
    }
    else {
      gsl_sf_result sin_xi_result;
      gsl_sf_result cos_xi_result;
      gsl_sf_sin_e(xi, &sin_xi_result);
      gsl_sf_cos_e(xi, &cos_xi_result);
      sin_term = sin_xi_result.val/xi;
      cos_term = cos_xi_result.val;
      sin_term_err = sin_xi_result.err/fabs(xi);
      cos_term_err = cos_xi_result.err;
    }
    if(eta < GSL_ROOT5_DBL_EPSILON) {
      coth_term = 1.0 + eta*eta/3.0 * (1.0 - eta*eta/15.0);
      sinh_term = 1.0 - eta*eta/6.0 * (1.0 - 7.0/60.0*eta*eta);
    }
    else {
      coth_term = eta/tanh(eta);
      sinh_term = eta/sinh(eta);
    }
    t1 = sqrt(lsqp1) * eta;
    pre_val = sinh_term/t1;
    pre_err = 2.0 * GSL_DBL_EPSILON * fabs(pre_val);
    term1 = sin_term*coth_term;
    term2 = cos_term;
    result->val  = pre_val * (term1 - term2);
    result->err  = pre_err * fabs(term1 - term2);
    result->err += pre_val * (sin_term_err * coth_term + cos_term_err);
    result->err += pre_val * fabs(term1-term2) * (fabs(eta) + 1.0) * GSL_DBL_EPSILON;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
}

int
gsl_sf_legendre_H3d_e(const int ell, const double lambda, const double eta,
                         gsl_sf_result * result)
{
  const double abs_lam = fabs(lambda);
  const double lsq     = abs_lam*abs_lam;
  const double xi      = abs_lam * eta;
  const double cosh_eta = cosh(eta);

  /* CHECK_POINTER(result) */

  if(eta < 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(eta > GSL_LOG_DBL_MAX) {
    /* cosh(eta) is too big. */
    OVERFLOW_ERROR(result);
  }
  else if(ell == 0) {
    return gsl_sf_legendre_H3d_0_e(lambda, eta, result);
  }
  else if(ell == 1) {
    return gsl_sf_legendre_H3d_1_e(lambda, eta, result);
  }
  else if(eta == 0.0) {
    result->val = 0.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(xi < 1.0) {
    return legendre_H3d_series(ell, lambda, eta, result);
  }
  else if((ell*ell+lsq)/sqrt(1.0+lsq)/(cosh_eta*cosh_eta) < 5.0*GSL_ROOT3_DBL_EPSILON) {
    /* Large argument.
     */
    gsl_sf_result P;
    double lm;
    int stat_P = gsl_sf_conicalP_large_x_e(-ell-0.5, lambda, cosh_eta, &P, &lm);
    if(P.val == 0.0) {
      result->val = 0.0;
      result->err = 0.0;
      return stat_P;
    }
    else {
      double lnN;
      gsl_sf_result lnsh;
      double ln_abslam;
      double lnpre_val, lnpre_err;
      int stat_e;
      gsl_sf_lnsinh_e(eta, &lnsh);
      legendre_H3d_lnnorm(ell, lambda, &lnN);
      ln_abslam = log(abs_lam);
      lnpre_val  = 0.5*(M_LNPI + lnN - M_LN2 - lnsh.val) - ln_abslam;
      lnpre_err  = lnsh.err;
      lnpre_err += 2.0 * GSL_DBL_EPSILON * (0.5*(M_LNPI + M_LN2 + fabs(lnN)) + fabs(ln_abslam));
      lnpre_err += 2.0 * GSL_DBL_EPSILON * fabs(lnpre_val);
      stat_e = gsl_sf_exp_mult_err_e(lnpre_val + lm, lnpre_err, P.val, P.err, result);
      return GSL_ERROR_SELECT_2(stat_e, stat_P);
    }
  }
  else if(abs_lam > 1000.0*ell*ell) {
    /* Large degree.
     */
    gsl_sf_result P;
    double lm;
    int stat_P = gsl_sf_conicalP_xgt1_neg_mu_largetau_e(ell+0.5,
                                                           lambda,
                                                           cosh_eta, eta,
                                                           &P, &lm);
    if(P.val == 0.0) {
      result->val = 0.0;
      result->err = 0.0;
      return stat_P;
    }
    else {
      double lnN;
      gsl_sf_result lnsh;
      double ln_abslam;
      double lnpre_val, lnpre_err;
      int stat_e;
      gsl_sf_lnsinh_e(eta, &lnsh);
      legendre_H3d_lnnorm(ell, lambda, &lnN);
      ln_abslam = log(abs_lam);
      lnpre_val  = 0.5*(M_LNPI + lnN - M_LN2 - lnsh.val) - ln_abslam;
      lnpre_err  = lnsh.err;
      lnpre_err += GSL_DBL_EPSILON * (0.5*(M_LNPI + M_LN2 + fabs(lnN)) + fabs(ln_abslam));
      lnpre_err += 2.0 * GSL_DBL_EPSILON * fabs(lnpre_val);
      stat_e = gsl_sf_exp_mult_err_e(lnpre_val + lm, lnpre_err, P.val, P.err, result);
      return GSL_ERROR_SELECT_2(stat_e, stat_P);
    }
  }
  else {
    /* Backward recurrence.
     */
    const double coth_eta = 1.0/tanh(eta);
    const double coth_err_mult = fabs(eta) + 1.0;
    gsl_sf_result rH;
    int stat_CF1 = legendre_H3d_CF1_ser(ell, lambda, coth_eta, &rH);
    double Hlm1;
    double Hl    = GSL_SQRT_DBL_MIN;
    double Hlp1  = rH.val * Hl;
    int lp;
    for(lp=ell; lp>0; lp--) {
      double root_term_0 = hypot(lambda,lp);
      double root_term_1 = hypot(lambda,lp+1.0);
      Hlm1 = ((2.0*lp + 1.0)*coth_eta*Hl - root_term_1 * Hlp1)/root_term_0;
      Hlp1 = Hl;
      Hl   = Hlm1;
    }

    if(fabs(Hl) > fabs(Hlp1)) {
      gsl_sf_result H0;
      int stat_H0 = gsl_sf_legendre_H3d_0_e(lambda, eta, &H0);
      result->val  = GSL_SQRT_DBL_MIN/Hl * H0.val;
      result->err  = GSL_SQRT_DBL_MIN/fabs(Hl) * H0.err;
      result->err += fabs(rH.err/rH.val) * (ell+1.0) * coth_err_mult * fabs(result->val);
      result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
      return GSL_ERROR_SELECT_2(stat_H0, stat_CF1);
    }
    else {
      gsl_sf_result H1;
      int stat_H1 = gsl_sf_legendre_H3d_1_e(lambda, eta, &H1);
      result->val  = GSL_SQRT_DBL_MIN/Hlp1 * H1.val;
      result->err  = GSL_SQRT_DBL_MIN/fabs(Hlp1) * H1.err;
      result->err += fabs(rH.err/rH.val) * (ell+1.0) * coth_err_mult * fabs(result->val);
      result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
      return GSL_ERROR_SELECT_2(stat_H1, stat_CF1);
    }
  }
}

int
gsl_sf_legendre_H3d_array(const int lmax, const double lambda, const double eta, double * result_array)
{
  /* CHECK_POINTER(result_array) */

 if(eta < 0.0 || lmax < 0) {
    int ell;
    for(ell=0; ell<=lmax; ell++) result_array[ell] = 0.0;
    GSL_ERROR ("domain error", GSL_EDOM);
  }
  else if(eta > GSL_LOG_DBL_MAX) {
    /* cosh(eta) is too big. */
    int ell;
    for(ell=0; ell<=lmax; ell++) result_array[ell] = 0.0;
    GSL_ERROR ("overflow", GSL_EOVRFLW);
  }
  else if(lmax == 0) {
    gsl_sf_result H0;
    int stat = gsl_sf_legendre_H3d_e(0, lambda, eta, &H0);
    result_array[0] = H0.val;
    return stat;
  }
  else {
    /* Not the most efficient method. But what the hell... it's simple.
     */
    gsl_sf_result r_Hlp1;
    gsl_sf_result r_Hl;
    int stat_lmax   = gsl_sf_legendre_H3d_e(lmax,   lambda, eta, &r_Hlp1);
    int stat_lmaxm1 = gsl_sf_legendre_H3d_e(lmax-1, lambda, eta, &r_Hl);
    int stat_max = GSL_ERROR_SELECT_2(stat_lmax, stat_lmaxm1);

    const double coth_eta = 1.0/tanh(eta);
    int stat_recursion = GSL_SUCCESS;
    double Hlp1 = r_Hlp1.val;
    double Hl   = r_Hl.val;
    double Hlm1;
    int ell;

    result_array[lmax]   = Hlp1;
    result_array[lmax-1] = Hl;

    for(ell=lmax-1; ell>0; ell--) {
      double root_term_0 = hypot(lambda,ell);
      double root_term_1 = hypot(lambda,ell+1.0);
      Hlm1 = ((2.0*ell + 1.0)*coth_eta*Hl - root_term_1 * Hlp1)/root_term_0;
      result_array[ell-1] = Hlm1;
      if(!(Hlm1 < GSL_DBL_MAX)) stat_recursion = GSL_EOVRFLW;
      Hlp1 = Hl;
      Hl   = Hlm1;
    }

    return GSL_ERROR_SELECT_2(stat_recursion, stat_max);
  }
}

/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/

#include "eval.h"

double gsl_sf_legendre_H3d_0(const double lambda, const double eta)
{
  EVAL_RESULT(gsl_sf_legendre_H3d_0_e(lambda, eta, &result));
}

double gsl_sf_legendre_H3d_1(const double lambda, const double eta)
{
  EVAL_RESULT(gsl_sf_legendre_H3d_1_e(lambda, eta, &result));
}

double gsl_sf_legendre_H3d(const int l, const double lambda, const double eta)
{
  EVAL_RESULT(gsl_sf_legendre_H3d_e(l, lambda, eta, &result));
}


text:
/**
 */
#include "aXe_grism.h"
#include <gsl/gsl_spline.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_sort_vector.h>
#include "fringe_conf.h"
#include "fringe_model.h"
#include "fringe_utils.h"

#define MIN(x,y) (((x)<(y))?(x):(y))
#define MAX(x,y) (((x)>(y))?(x):(y))

/**
 * Function: fringe_correct_PET
 * Performs fringe corrections to PET pixel lists
 *
 * The main function to compute and perform the fringe
 * correction to a list of PET pixels. for every pixel in
 * the PET list the pixel throughput function is
 * computed, and the relevant information from
 * the fringing model is extracted. Then the fringe
 * amplitude is computed and the count value of the pixel
 * is corrected. Also the PET pixel list is simultaneously
 * parsed and corrected for the fringing.
 *
 * Parameters:
 * @param fconf    - the fringe configuration file
 * @param act_beam - the beam aperture
 * @param obj_pet  - the PET pixel
 * @param bck_pet  - the PET pixel
 */
void
fringe_correct_PET(const fringe_conf *fconf, const beam act_beam,
		   ap_pixel *obj_pet, ap_pixel *bck_pet)
{
  ap_pixel *obj_pix=NULL;
  ap_pixel *bck_pix=NULL;

  //d_point angles;
  pixel_tput *p_tput ;
  gsl_vector **tput_vectors=NULL;

  //int ii=0;
  //int jj=0;
  int index;
  int pix_index;

  double lambda_mean;
  double pixel_ampl;
  double fringe_factor;
  double fringe_tot;

  optical_property *optprops;

  // allocate the pixel throughput structure
  p_tput = alloc_pixel_tput();

  // allocate memory for the optical property structure
  optprops = alloc_optprops_list(fconf);

  // compute the mean wavelength
  lambda_mean = (gsl_vector_get(fconf->fringe_range, 0) +
    (gsl_vector_get(fconf->fringe_range, 1)
     - gsl_vector_get(fconf->fringe_range, 0))/2.0)/1.0e+04;

  // initialize some values in the optical property list
  init_optprops_list(fconf, lambda_mean, optprops);

  // initialize the PET pixel storages
  obj_pix = obj_pet;
  bck_pix = bck_pet;

  fringe_tot = 0.0;
  pix_index = 0;
  // go over the whole PET list
  while (obj_pix->p_x != -1)
    {
      // check whether the pixel has to be fringe corrected
      // at all, based on boundaries given in the configuration
      if ( obj_pix->lambda < gsl_vector_get(fconf->fringe_range, 0)
      	   || obj_pix->lambda > gsl_vector_get(fconf->fringe_range, 1)
      	   || obj_pix->dlambda > fconf->max_dispersion)
      	{

      	  // count up the object PET
      	  obj_pix++;

      	  // count up the background PET
      	  if (bck_pet)
      	    bck_pix++;

      	  // jump to the next PET pixel
      	  continue;
      	}

      // determine the filter throughput values
      //      tput_vectors = evaluate_pixel_throughput(fconf, act_beam, obj_pix);
      // for the ISR with the WFC value
      tput_vectors = get_gauss_throughput(fconf, obj_pix, 80.0);

      // fill the optical thickness of the layers
      // into the structure
      fill_optprops_thickness(fconf->opt_layers, obj_pix->p_x,
			      obj_pix->p_y, optprops);

      pixel_ampl = 0.0;
      for (index=0; index < (int)tput_vectors[0]->size; index++)
	{
	  // fill all information in the optical
	  // property list
	  fill_optprops_all(fconf->opt_layers,
			    gsl_vector_get(tput_vectors[0],index),
			    optprops);

	  // compute and add the contribution at a wavelength
	  pixel_ampl += gsl_vector_get(tput_vectors[1],index)*
	    fringe_contrib_single(optprops, fconf);
	}

      //put together the fringe correction factor
      fringe_factor = fconf->fringe_amp * pixel_ampl + 1.0;
#ifdef DEBUGFCONF
      fprintf(stderr, "Wavelength: %f, dispersion: %f, fringe factor: %f, (x,y): (%i, %i)\n",
	      obj_pix->lambda, obj_pix->dlambda, fringe_factor, obj_pix->p_x, obj_pix->p_y);
#endif

      // correct the object PET
      obj_pix->count /= fringe_factor;
      //obj_pix->count = fringe_factor;
      //fprintf(stdout, "%e ",obj_pix->count);

      fringe_tot += fringe_factor;
      pix_index++;

      // treat the background PET pixel
      if (bck_pet)
	{
	  if (bck_pix->p_x != obj_pix->p_x || bck_pix->p_y != obj_pix->p_y)
	    aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
			 "aXe_FRINGECORR: The sequence of pixels\n"
			 "in the object PET and background PET differs.\n"
			 "Something must be wrong!\n");

	  // correct also the background PET
	  bck_pix->count /= fringe_factor;

	  // count up the background PET
	  bck_pix++;
	}

      // count up the object PET
      obj_pix++;
      //      obj_pix->p_x = -1;
    }

  //  if (pix_index)
  //    fprintf(stdout,"Mean fringe factor: %f\n", fringe_tot/(double)pix_index);

  if (tput_vectors)
    {
      // release the memory in the vectors
      gsl_vector_free(tput_vectors[0]);
      gsl_vector_free(tput_vectors[1]);
      free(tput_vectors);
    }

  // free the optical property structure
  free_optprops_list(optprops);

  // free the memory
  free_pixel_tput(p_tput);
}

/**
 * Function: evaluate_pixel_throughput
 * Computes the wavelength and filter throughputs for a PET pixel.
 *
 * The function computes and returns two vector, one with wavelength
 * steps and the second with the pixel throughputs at those wavelengths.
 * This data is computed to a single PET pixel according to the settings
 * in the fringe configuration file.
 *
 * Parameters:
 * @param fconf    - the fringe configuration file
 * @param act_beam - the beam aperture
 * @param obj_pix  - the PET pixel
 *
 * Returns:
 * @return double_vector - the wavelength and filter throughput steps
 */
gsl_vector **
evaluate_pixel_throughput(const fringe_conf *fconf,
			  const beam act_beam, const ap_pixel *obj_pix)
{
  pixel_tput *p_tput ;

  gsl_vector **double_vector;
  gsl_vector *lambda_values;
  gsl_vector *through_values;

  interpolator *interp;

  double stepsize;
  double lambda_act;
  double through_act;
  double through_tot;

  double abs_min;
  double abs_max;

  int index;

  // allocate the pixel throughput structure
  p_tput = alloc_pixel_tput();

  // allocate the space for the vectors
  lambda_values  = gsl_vector_alloc(fconf->num_steps + 1);
  through_values = gsl_vector_alloc(fconf->num_steps + 1);

  // allocate space for the return vector
  double_vector  = (gsl_vector **)malloc(2*sizeof (gsl_vector *));

  compute_pixel_tput(act_beam, obj_pix, p_tput);
#ifdef DEBUGFCONF
  print_pixel_tput(p_tput);
#endif
  // convert the pixel throughput information
  // into an interplator
  interp = create_interp_tput(p_tput);

  abs_min = MAX(interp->xmin, gsl_vector_get(fconf->fringe_range, 0));
  abs_max = MIN(interp->xmax, gsl_vector_get(fconf->fringe_range, 1));

  // compute the steps in wavelength
  stepsize = (abs_max - abs_min)/(double)fconf->num_steps;

  // set the lower wavelength end
  lambda_act = abs_min;

  // loop over the wavelength
  through_tot = 0.0;
  for (index=0; index < fconf->num_steps; index++)
    {
      // set the wavelength value
      gsl_vector_set(lambda_values, index, lambda_act);

      // compute the throughput
      through_act = eval_interp(interp, lambda_act);

      // compute the total throughput
      through_tot += through_act;

      // set the thoughput value
      gsl_vector_set(through_values, index, through_act);

      lambda_act +=  stepsize;
    }

  // treat the last step separately, avoiding
  // to leave the allowed interpolation range
  // due to rounding errors
  gsl_vector_set(lambda_values, fconf->num_steps, abs_max);
  through_act = eval_interp(interp, abs_max);
  through_tot += through_act;
  gsl_vector_set(through_values, fconf->num_steps, through_act);

  for (index=0; index < (int)through_values->size; index++)
    {
      // normalize the filter throughput values
      gsl_vector_set(through_values, index,
		     gsl_vector_get(through_values, index)/through_tot);

      // convert the wavelength from AA to micron
      gsl_vector_set(lambda_values, index,
		     gsl_vector_get(lambda_values, index)*1.0e-04);
    }
#ifdef DEBUGFCONF
  through_tot = 0.0;
  for (index=0; index < through_values->size; index++)
    {
      fprintf(stderr, "Wavelength: %f, Throughput: %f\n",
	      gsl_vector_get(lambda_values, index),
	      gsl_vector_get(through_values, index));
      through_tot += gsl_vector_get(through_values, index);
    }
  fprintf(stderr, "Total throughput: %f\n", through_tot);
#endif

  // build up the output array
  double_vector[0] = lambda_values;
  double_vector[1] = through_values;

  // free the memory of the interpolator
  free_interp(interp);

  // free the memory of the pixel throughput
  free_pixel_tput(p_tput);

  // return the two gsl vectors
  return double_vector;
}

/**
 * Function: get_gauss_throughput
 *
 * Parameters:
 * @param fconf    - the fringe configuration file
 * @param act_beam - the beam aperture
 * @param obj_pix  - the PET pixel
 *
 * Returns:
 * @return double_vector - the wavelength and filter throughput steps
 */
gsl_vector **
get_gauss_throughput(const fringe_conf *fconf, const ap_pixel *obj_pix,
		     const double fwhm)
{
  pixel_tput *p_tput ;

  gsl_vector **double_vector;
  gsl_vector *lambda_values;
  gsl_vector *through_values;

  //interpolator *interp;

  double stepsize;
  double lambda_act;
  double through_act;
  double through_tot;

  double abs_min;
  double abs_max;

  double sigma;
  double arg;

  int index;

  // convert the FWHM into sigma
  sigma = fwhm / 2.35482;

  // allocate the pixel throughput structure
  p_tput = alloc_pixel_tput();

  // allocate the space for the vectors
  lambda_values  = gsl_vector_alloc(fconf->num_steps + 1);
  through_values = gsl_vector_alloc(fconf->num_steps + 1);

  // allocate space for the return vector
  double_vector  = (gsl_vector **)malloc(2*sizeof (gsl_vector *));

  // compute the minimum and maximum wavelength
  abs_min = MAX(obj_pix->lambda - 2.5 * fwhm,
		gsl_vector_get(fconf->fringe_range, 0));
  abs_max = MIN(obj_pix->lambda + 2.5 * fwhm,
		gsl_vector_get(fconf->fringe_range, 1));

  // compute the steps in wavelength
  stepsize = (abs_max - abs_min)/(double)fconf->num_steps;

  // set the lower wavelength end
  lambda_act = abs_min;

  // loop over the wavelength
  through_tot = 0.0;
  for (index=0; index < fconf->num_steps; index++)
    {
      // set the wavelength value
      gsl_vector_set(lambda_values, index, lambda_act);

      // z=(coord-g(2))/g(3)
      // sumgauss=sumgauss+DBLE(g(1)*EXP(-0.5*z**2))
      arg = (lambda_act - obj_pix->lambda) / sigma;

      // compute the throughput
      through_act = exp(-0.5*arg*arg);

      // compute the total throughput
      through_tot += through_act;

      // set the thoughput value
      gsl_vector_set(through_values, index, through_act);

      // increment the wavelength
      lambda_act += stepsize;
    }

  // treat the last step separately, avoiding
  // to leave the allowed interpolation range
  // due to rounding errors
  gsl_vector_set(lambda_values, fconf->num_steps, abs_max);
  arg = (abs_max - obj_pix->lambda) / sigma;
  through_act = exp(-0.5*arg*arg);
  through_tot += through_act;
  gsl_vector_set(through_values, fconf->num_steps, through_act);

  for (index=0; index < (int)through_values->size; index++)
    {
      // normalize the filter throughput values
      gsl_vector_set(through_values, index,
		     gsl_vector_get(through_values, index)/through_tot);

      // convert the wavelength from AA to micron
      gsl_vector_set(lambda_values, index,
		     gsl_vector_get(lambda_values, index)*1.0e-04);
    }

#ifdef DEBUGFCONF
  through_tot = 0.0;
  for (index=0; index < through_values->size; index++)
    {
      fprintf(stderr, "Wavelength: %f, Throughput: %f\n",
	      gsl_vector_get(lambda_values, index),
	      gsl_vector_get(through_values, index));
      through_tot += gsl_vector_get(through_values, index);
    }
  fprintf(stderr, "Total throughput: %f\n", through_tot);
#endif

  // build up the output array
  double_vector[0] = lambda_values;
  double_vector[1] = through_values;

  return double_vector;
}

/**
 * Function: create_interp_tput
 * The function creates and returns an interpolator
 * for a pixel throughput function. The interpolator
 * structure is allocated, data is filled in and
 * finally initialized.
 *
 * Parameters:
 * @param p_tput   - the pixel throughput structure
 *
 * Returns:
 * @return interp - the new interpolator
 */
interpolator *
create_interp_tput(const pixel_tput *p_tput)
{
  interpolator *interp;

  double *xvals;
  double *yvals;

  // allocate memory for the independent values
  xvals = (double *) malloc(4 * sizeof(double));
  if (!xvals) {
    aXe_message (aXe_M_ERROR, __FILE__, __LINE__,
		 "Memory allocation failed");
  }

  // allocate memory for the dependent values
  yvals = (double *) malloc(4 * sizeof(double));
  if (!yvals) {
    aXe_message (aXe_M_ERROR, __FILE__, __LINE__,
		 "Memory allocation failed");
  }

  // set the independent vector values
  xvals[0] = gsl_vector_get(p_tput->lambda_values, 0);
  xvals[1] = gsl_vector_get(p_tput->lambda_values, 1);
  xvals[2] = gsl_vector_get(p_tput->lambda_values, 2);
  xvals[3] = gsl_vector_get(p_tput->lambda_values, 3);

  // set the dependent vector values
  yvals[0] = 0.0;
  yvals[1] = p_tput->tp_max;
  yvals[2] = p_tput->tp_max;
  yvals[3] = 0.0;

  // create the interpolator
  interp = create_interp(4, FILTER_INTERP_TYPE, xvals, yvals);

  // return the interpolator
  return interp;
}

/**
 * Function: compute_pixel_tput
 * The function computes all relevant data to characterize
 * the pixel throughput function for an individual PET pixel.
 * The information is written into a pixel throughput structure.
 *
 * Parameters:
 * @param act_beam - the beam aperture
 * @param obj_pix  - the PET pixel
 * @param p_tput   - the pixel throughput structure
 */
void
compute_pixel_tput(const beam act_beam, const ap_pixel *obj_pix,
		   pixel_tput *p_tput)
{
  d_point angles;

  double theta_2a;
  double theta_3a;

  double theta_1b;
  double theta_2b;

  double sin_theta_2a;
  double sin_theta_3a;

  double sin_theta_1b;
  double cos_theta_1b;
  double sin_theta_2b;

  double dist1;
  double dist2;

  // compute the anngles alpha and beta
  angles = compute_tput_angles(act_beam, obj_pix);

  // the computation of p_min-p3 in the
  // design document
  theta_1b = M_PI_2 + angles.x - angles.y;
  theta_2b = angles.y;

  // compute the sin's and cos's neede
  sin_theta_1b = sin(theta_1b);
  cos_theta_1b = cos(theta_1b);
  sin_theta_2b = sin(theta_2b);

  // compute one delta
  dist2 = obj_pix->dlambda * sin_theta_1b / sin_theta_2b;

  // the computation of p_h-p0 in the
  // design document
  theta_2a = M_PI - angles.y;
  theta_3a = angles.y - angles.x - M_PI_4;

  // next line due to sin(angle) = sin(180-angle)
  sin_theta_2a = sin_theta_2b;
  sin_theta_3a = sin(theta_3a);

  // compute the second delta
  dist1 = obj_pix->dlambda * M_SQRT1_2 * sin_theta_3a / sin_theta_2a;

  // determine the maximum throughput
  p_tput->tp_max = 1.0/MAX(fabs(sin_theta_1b), fabs(cos_theta_1b));

  // put together the distances, using
  // the deltas and the dispersion
  gsl_vector_set(p_tput->lambda_values, 0, obj_pix->lambda - dist1 - dist2);
  gsl_vector_set(p_tput->lambda_values, 1, obj_pix->lambda - dist1);
  gsl_vector_set(p_tput->lambda_values, 2, obj_pix->lambda + dist1);
  gsl_vector_set(p_tput->lambda_values, 3, obj_pix->lambda + dist1 + dist2);

  // sort the distances
  gsl_sort_vector(p_tput->lambda_values);
}

/**
 * Function: compute_tput_angles
 * The function computes the angles which characterize
 * a pixel in the lambda-cross dispersion plane.
 * In detail this is the angle between the lambda-axis
 * and the x-axis, and the angle between the lambda abd
 * the cross dispersion axis.
 * The two angles are computed from the basic pixel
 * and beam data.
 *
 *
 * Parameters:
 * @param act_beam - the beam aperture
 * @param act_pet  - the PET pixel
 *
 * Returns:
 * @return angles - the relevant pixel angles
 */
d_point
compute_tput_angles(const beam act_beam, const ap_pixel *act_pet)
{
  d_point angles;

  // compute the angle 'alpha' between lambda-axis and
  // x-axis
  angles.x = -atan2(act_beam.spec_trace->deriv(act_pet->xi,
					       act_beam.spec_trace->data),1.0);
  // compute the angle beta between lambda and cross dispersion axis
  angles.y = act_beam.orient + angles.x;

  // return the two angles
  return angles;
}

/**
 * Function: alloc_pixel_tput
 * The function allocates and returns memory
 * for a pixel throughput function
 *
 * Returns:
 * @return p_tput -the new allocated pixel throughput structure
 */
pixel_tput *
alloc_pixel_tput()
{
  pixel_tput *p_tput;

  // allcoate the pointer
  p_tput = (pixel_tput *) malloc(sizeof(pixel_tput));

  // allocate the vector in the structure
  p_tput->lambda_values = gsl_vector_alloc(4);

  // return the structure
  return p_tput;
}

/**
 * Function: print_pixel_tput
 * The function prints the content of a pixel throughput
 * structure onto the screen.
 *
 * Parameters:
 * @param p_tput - the pixel throughput structure
 */
void
print_pixel_tput(pixel_tput *p_tput)
{
  fprintf(stdout, "l_min: %f, l_1: %f, l_2: %f, l_max: %f; tp_max: %f\n\n",
	  gsl_vector_get(p_tput->lambda_values,0),
	  gsl_vector_get(p_tput->lambda_values,1),
	  gsl_vector_get(p_tput->lambda_values,2),
	  gsl_vector_get(p_tput->lambda_values,3),
	  p_tput->tp_max);
}

/**
 * Function: free_pixel_tput
 * The function releases the memory allocated
 * in a pixel throughput function.
 *
 * Parameters:
 * @param p_tput - the pixel throughput structure
 */
void
free_pixel_tput(pixel_tput *p_tput)
{
  // free the gsl-vector
  gsl_vector_free(p_tput->lambda_values);

  // free everything
  free(p_tput);

  // set it to NULL
  p_tput = NULL;
}

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe


text:
/**
 *
 * @file qwrapper_dgemm_tile.c
 *
 *  PLASMA core_blas quark wrapper
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Mark Gates
 * @date 2010-11-15
 * @generated d Tue Jan  7 11:44:59 2014
 *
 **/
#include <cblas.h>
#include "common.h"

/***************************************************************************//**
 *
 * Version of dgemm for tile storage, to avoid dependency problem when
 * computations are done within the tile. alpha and beta are passed as
 * pointers so they can depend on runtime values.
 *
 * @param[in] Alock
 *          Pointer to tile owning submatrix A.
 *
 * @param[in] Block
 *          Pointer to tile owning submatrix B.
 *
 * @param[in] Clock
 *          Pointer to tile owning submatrix C.
 *
 **/
void QUARK_CORE_dgemm_tile(Quark *quark, Quark_Task_Flags *task_flags,
                           PLASMA_enum transA, PLASMA_enum transB,
                           int m, int n, int k, int nb,
                           const double *alpha, const double *A, int lda,
                                                            const double *B, int ldb,
                           const double *beta,        double *C, int ldc,
                           const double *Alock,
                           const double *Block,
                           const double *Clock)
{
    DAG_CORE_GEMM;
    QUARK_Insert_Task(quark, CORE_dgemm_tile_quark, task_flags,
        sizeof(PLASMA_enum),              &transA, VALUE,
        sizeof(PLASMA_enum),              &transB, VALUE,
        sizeof(int),                      &m,      VALUE,
        sizeof(int),                      &n,      VALUE,
        sizeof(int),                      &k,      VALUE,
        sizeof(double),       alpha,           INPUT,
        sizeof(double)*nb*nb, A,               NODEP,          /* input; see Alock */
        sizeof(int),                      &lda,    VALUE,
        sizeof(double)*nb*nb, B,               NODEP,          /* input; see Block */
        sizeof(int),                      &ldb,    VALUE,
        sizeof(double),       beta,            INPUT,
        sizeof(double)*nb*nb, C,                       NODEP,  /* inout; see Clock */
        sizeof(int),                      &ldc,    VALUE,
        sizeof(double)*nb*nb, Alock,           INPUT,
        sizeof(double)*nb,    Block,           INPUT,
        sizeof(double)*nb,    Clock,                   INOUT,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_dgemm_tile_quark = PCORE_dgemm_tile_quark
#define CORE_dgemm_tile_quark PCORE_dgemm_tile_quark
#endif
void CORE_dgemm_tile_quark(Quark *quark)
{
    PLASMA_enum transA, transB;
    int m, n, k, lda, ldb, ldc;
    const double *alpha, *beta;
    const double *A, *B;
    double *C;

    quark_unpack_args_13( quark, transA, transB, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc );
    cblas_dgemm(
        CblasColMajor,
        (CBLAS_TRANSPOSE)transA, (CBLAS_TRANSPOSE)transB,
        m, n, k,
        (*alpha), A, lda,
                             B, ldb,
        (*beta),  C, ldc );
}