First git commit of the project

gcc -O3 -Wall -o bl_pol main_program.c hcubature.c restframe_quantities.c polarization_angle.c numpar.c stringpar.c strrev.c trig_funct.c -lm

/* Body of convergence test, shared between hcubature.c and

   pcubature.c.  We use an #include file because the two routines use

   somewhat different data structures, and define macros ERR(j) and

   VAL(j) to get the error and value estimates, respectively, for

   integrand j. */

{

     unsigned j;

#    define SQR(x) ((x) * (x))

     switch (norm) {

	 case ERROR_INDIVIDUAL:

	      for (j = 0; j < fdim; ++j)

		   if (ERR(j) > reqAbsError && ERR(j) > fabs(VAL(j))*reqRelError)

			return 0;

	      return 1;

	 case ERROR_PAIRED:

	      for (j = 0; j+1 < fdim; j += 2) {

		   double maxerr, serr, err, maxval, sval, val;

		   /* scale to avoid overflow/underflow */

		   maxerr = ERR(j) > ERR(j+1) ? ERR(j) : ERR(j+1);

		   maxval = VAL(j) > VAL(j+1) ? VAL(j) : VAL(j+1);

		   serr = maxerr > 0 ? 1/maxerr : 1;

		   sval = maxval > 0 ? 1/maxval : 1;

		   err = sqrt(SQR(ERR(j)*serr) + SQR(ERR(j+1)*serr)) * maxerr;

		   val = sqrt(SQR(VAL(j)*sval) + SQR(VAL(j+1)*sval)) * maxval;

		   if (err > reqAbsError && err > val*reqRelError)

			return 0;

	      }

	      if (j < fdim) /* fdim is odd, do last dimension individually */

		   if (ERR(j) > reqAbsError && ERR(j) > fabs(VAL(j))*reqRelError)

			return 0;

	      return 1;

	 case ERROR_L1: {

	      double err = 0, val = 0;

	      for (j = 0; j < fdim; ++j) {

		   err += ERR(j);

		   val += fabs(VAL(j));

	      }

	      return err <= reqAbsError || err <= val*reqRelError;

	 }

	 case ERROR_LINF: {

	      double err = 0, val = 0;

	      for (j = 0; j < fdim; ++j) {

		   double absval = fabs(VAL(j));

		   if (ERR(j) > err) err = ERR(j);

		   if (absval > val) val = absval;

	      }

	      return err <= reqAbsError || err <= val*reqRelError;

	 }

	 case ERROR_L2: {

	      double maxerr = 0, maxval = 0, serr, sval, err = 0, val = 0;

	      /* scale values by 1/max to avoid overflow/underflow */

	      for (j = 0; j < fdim; ++j) {

		   double absval = fabs(VAL(j));

		   if (ERR(j) > maxerr) maxerr = ERR(j);

		   if (absval > maxval) maxval = absval;

	      }

	      serr = maxerr > 0 ? 1/maxerr : 1;

	      sval = maxval > 0 ? 1/maxval : 1;

	      for (j = 0; j < fdim; ++j) {

		   err += SQR(ERR(j) * serr);

		   val += SQR(fabs(VAL(j)) * sval);

	      }

	      err = sqrt(err) * maxerr;

	      val = sqrt(val) * maxval;

	      return err <= reqAbsError || err <= val*reqRelError;

	 }

     }

     return 1; /* unreachable */

}

/* Adaptive multidimensional integration of a vector of integrands.

 *

 * Copyright (c) 2005-2013 Steven G. Johnson

 *

 * Portions (see comments) based on HIntLib (also distributed under

 * the GNU GPL, v2 or later), copyright (c) 2002-2005 Rudolf Schuerer.

 *     (http://www.cosy.sbg.ac.at/~rschuer/hintlib/)

 *

 * Portions (see comments) based on GNU GSL (also distributed under

 * the GNU GPL, v2 or later), copyright (c) 1996-2000 Brian Gough.

 *     (http://www.gnu.org/software/gsl/)

 *

 * This program is free software; you can redistribute it and/or modify

 * it under the terms of the GNU General Public License as published by

 * the Free Software Foundation; either version 2 of the License, or

 * (at your option) any later version.

 *

 * This program is distributed in the hope that it will be useful,

 * but WITHOUT ANY WARRANTY; without even the implied warranty of

 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

 * GNU General Public License for more details.

 *

 * You should have received a copy of the GNU General Public License

 * along with this program; if not, write to the Free Software

 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

 *

 */

#ifndef CUBATURE_H

#define CUBATURE_H

#include <stdlib.h> /* for size_t */

#ifdef __cplusplus

extern "C"

{

#endif /* __cplusplus */

/* USAGE: Call hcubature or pcubature with your function as described

          in the README file. */

/* a vector integrand - evaluates the function at the given point x

   (an array of length ndim) and returns the result in fval (an array

   of length fdim).   The void* parameter is there in case you have

   to pass any additional data through to your function (it corresponds

   to the fdata parameter you pass to cubature).  Return 0 on

   success or nonzero to terminate the integration. */

typedef int (*integrand) (unsigned ndim, const double *x, void *,

                          unsigned fdim, double *fval);

/* a vector integrand of a vector of npt points: x[i*ndim + j] is the

   j-th coordinate of the i-th point, and the k-th function evaluation

   for the i-th point is returned in fval[i*fdim + k].  Return 0 on success

   or nonzero to terminate the integration. */

typedef int (*integrand_v) (unsigned ndim, size_t npt,

			    const double *x, void *,

			    unsigned fdim, double *fval);

/* Different ways of measuring the absolute and relative error when

   we have multiple integrands, given a vector e of error estimates

   in the individual components of a vector v of integrands.  These

   are all equivalent when there is only a single integrand. */

typedef enum {

     ERROR_INDIVIDUAL = 0, /* individual relerr criteria in each component */

     ERROR_PAIRED, /* paired L2 norms of errors in each component,

		      mainly for integrating vectors of complex numbers */

     ERROR_L2, /* abserr is L_2 norm |e|, and relerr is |e|/|v| */

     ERROR_L1, /* abserr is L_1 norm |e|, and relerr is |e|/|v| */

     ERROR_LINF /* abserr is L_\infty norm |e|, and relerr is |e|/|v| */

} error_norm;

/* Integrate the function f from xmin[dim] to xmax[dim], with at most

   maxEval function evaluations (0 for no limit), until the given

   absolute or relative error is achieved.  val returns the integral,

   and err returns the estimate for the absolute error in val; both

   of these are arrays of length fdim, the dimension of the vector

   integrand f(x). The return value of the function is 0 on success

   and non-zero if there  was an error. */

/* adapative integration by partitioning the integration domain ("h-adaptive")

   and using the same fixed-degree quadrature in each subdomain, recursively,

   until convergence is achieved. */

int hcubature(unsigned fdim, integrand f, void *fdata,

	      unsigned dim, const double *xmin, const double *xmax, 

	      size_t maxEval, double reqAbsError, double reqRelError, 

	      error_norm norm,

	      double *val, double *err);

/* as hcubature, but vectorized integrand */

int hcubature_v(unsigned fdim, integrand_v f, void *fdata,

		unsigned dim, const double *xmin, const double *xmax, 

		size_t maxEval, double reqAbsError, double reqRelError, 

		error_norm norm,

		double *val, double *err);

/* adaptive integration by increasing the degree of (tensor-product

   Clenshaw-Curtis) quadrature rules ("p-adaptive"), rather than

   subdividing the domain ("h-adaptive").  Possibly better for

   smooth integrands in low dimensions. */

int pcubature_v_buf(unsigned fdim, integrand_v f, void *fdata,

		    unsigned dim, const double *xmin, const double *xmax,

		    size_t maxEval, 

		    double reqAbsError, double reqRelError,

		    error_norm norm,

		    unsigned *m,

		    double **buf, size_t *nbuf, size_t max_nbuf,

		    double *val, double *err);

int pcubature_v(unsigned fdim, integrand_v f, void *fdata,

		unsigned dim, const double *xmin, const double *xmax, 

		size_t maxEval, double reqAbsError, double reqRelError, 

		error_norm norm,

		double *val, double *err);

int pcubature(unsigned fdim, integrand f, void *fdata,

	      unsigned dim, const double *xmin, const double *xmax, 

	      size_t maxEval, double reqAbsError, double reqRelError, 

	      error_norm norm,

	      double *val, double *err);

#ifdef __cplusplus

}  /* extern "C" */

#endif /* __cplusplus */

#endif /* CUBATURE_H */

#include <ctype.h>

#include <math.h>

#include <stdio.h>

#include <stdlib.h>

#include <stdint.h>

#include <string.h>

#include "cubature.h"

#define TRUE 1

#define FALSE 0

#define VERBOSE 0

#define NR_END 1

#define FREE_ARG char*

#define C_LIGHT 3e10

#define NS_RADIUS 1E6

#define PI 3.14159

#define enep 2.73

extern char* poltype;

extern char* simulfile;

double blackbody(double energy, double ktbb);

double IntensityRestFrame(double E_rest, double CosAlphaPrime, double ktbb);

double FunctERest(double E_obs, double delta, double u);

double FunctCosPsi(double i_obs, double theta, double phi);

double lensing_factor(double compact, double xval);

double dCosapha_dCosPsi(double compact, double xval);

double FunctCosAlpha(double u, double CosPsi);

double FunctCosPsi(double i_obs, double theta, double phi);

int observed_flux(unsigned dim, const double* x, void* params, unsigned fdim, double* retval);

double doppler_factor(double spin, double compact, double cos_alpha, double Psi, double i_obs, double theta, double Phi);

double polarization_angle(double alpha, double Psi, double theta, double phi, double i_obs, double beta);

double PolarizationRestFrame(double E_rest, double CosAlphaPrime);

double FunctBetaValue (double spin, double u, double theta);

double numpar(char* parameter);

char* stringpar(char* input_line);

char* strrev(char* string);

double cot(double a);

double csc(double a);

double arctan(double yq, double xq);

double reduce_azim(double phi);