text:
/*GENLIB library routines */
/* Compile progs: gcc -O3 -o prog prog.c ~/genlib.c -lm -lgsl -lgslcblas */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

#define pi 3.14159265358979
#define infinity 999999
#define true 1
#define false 0
#define maxranges 5                    /* for numerical integration */
#define maxsimpsonpoints 1025          /* change libhdr if these changed */

#define MBIG 1000000000
#define MSEED 161803398
#define FAC (1.0/MBIG)

#define ib1 1
#define ib2 2
#define ib5 16
#define ib18 131072

#define undefined -99999999999999.9999999

typedef double (*fun1)();

FILE *seedfileptr, *tmoutfile, *fopen();

int inputcounter=0;

/* GSL random number definitions */
gsl_rng * gsl_rng_r;
unsigned long int seed;
int random_number_generator_initialised_flag = 0;

long sseed;  /* seed for the slow generator */
int tracelevel;

struct acc
{
   int n;
   double sum;
   double sumsq;
};

struct covacc
{
   int n;
   double sumxy, sumx, sumy, sumx2, sumy2;
};

gabort(char *s, int r)
{
   printf("ERROR %s %d\n", s, r);
   exit(1);
}

/* This is the random number generator ran3 from Numerical Recipes from
Knuth.  It's the fastest of the ``non-quick-and-dirty'' generators, and
is supposed to have decent properties */
double uniform_old()
{
	static int inext,inextp;
	static long ma[56];
	static int iff=0;
	long mj,mk;
        double res;
	int i,ii,k;

	if(sseed<0||iff==0)
	{
		iff=1;
		mj=MSEED-(sseed<0?-sseed:sseed);
		mj%=MBIG;
		ma[55]=mj;
		mk=1;
		for(i=1;i<=54;i++)
		{
			ii=(21*i)%55;
			ma[ii]=mk;
			mk=mj-mk;
			if(mk<0) mk+=MBIG;
			mj=ma[ii];
		}
		for(k=1;k<=4;k++)
			for(i=1;i<=55;i++)
			{
				ma[i] -=ma[1+(i+30)%55];
				if(ma[i]<0) ma[i]+=MBIG;
			}
		inext=0;
		inextp=31;
		sseed=1;
	}
	if(++inext==56) inext=1;
	if(++inextp==56) inextp=1;
	mj=ma[inext]-ma[inextp];
	if(mj<0) mj+=MBIG;
	ma[inext]=mj;
	sseed=mj;
        res = mj*FAC;
        if (res > 1.0) res = 1.0;
        else if (res < 0.0) res = 0.0;
	return res;
}

double uniform()
{
   double res;
   if (!random_number_generator_initialised_flag)
      gabort("uniform() called, but !random_number_generator_initialised_flag", 0);
   res = gsl_ran_flat(gsl_rng_r,0.0,1.0);
   return(res);
}

int irbit1(unsigned long *iseed)
/* routine to return random bits */
{
  unsigned long newbit;
  newbit = (*iseed &ib18) >> 17
  ^ (*iseed &ib5) >> 4
  ^ (*iseed &ib2) >> 1
  ^ (*iseed &ib1);
  *iseed = (*iseed<<1) | newbit;
  return (int)newbit;
}

initialise_uniform_generator()
{
/* create a generator chosen by the environment variable GSL_RNG_TYPE */
   const gsl_rng_type * T;
//   printf("initialise_uniform_generator\n");
   gsl_rng_env_setup();
   T = gsl_rng_default;
   gsl_rng_r = gsl_rng_alloc (T);
   random_number_generator_initialised_flag = 1;
}

getseed()
{
   FILE *seedfileptr;
   static int c = 1;
   if (c==1)         // first time routine called
   {
      initialise_uniform_generator();
      c = 0;
   }
   printf("Enter seed (-99 to read from file) ");
   manual: scanf("%lu", &seed);
   if (seed == -99)
   {
      seedfileptr = fopen("seedfile", "r");
      if (seedfileptr==0)
      {
         printf("No seedfile, enter seed please ");
         goto manual;
      }
      fscanf(seedfileptr, "%lu", &seed);
//      printf("Seed read %lu\n", seed);
      fclose(seedfileptr);
   }
   gsl_rng_set(gsl_rng_r, seed);
}

getseedquick()
{
   FILE *seedfileptr;
   static int c = 1;
   if (c==1)         // first time routine called
   {
      initialise_uniform_generator();
      c = 0;
   }
   seedfileptr = fopen("seedfile", "r");
   if (seedfileptr==0)
   {
      printf("No seedfile, enter seed please ");
      scanf("%lu", &seed);
   }
   else
   {
      fscanf(seedfileptr, "%lu", &seed);
      fclose(seedfileptr);
   }
//   printf("getseedquick: Seed read %lu\n", seed);
   gsl_rng_set(gsl_rng_r, seed);
}

terminate_uniform_generator()
{
   gsl_rng_free(gsl_rng_r);
}

writeseed()
{
   FILE *seedfileptr;
   double temp, x;
   x = uniform();
   temp = floor(x*100000000);
   seed = (unsigned long int)temp;
//   printf("Write seed %lu\n", seed);
   seedfileptr = fopen("seedfile", "w");
   fprintf(seedfileptr, "%lu\n", seed);
   fclose(seedfileptr);
   terminate_uniform_generator();
}

initacc(a)
struct acc *a;
{
   a -> n = 0;
   a -> sum = 0;
   a -> sumsq = 0;
}

accum(struct acc *a, double x)
{
   a->n = a->n + 1;
   a->sum = a->sum + x;
   a->sumsq = a->sumsq + x*x;
}

double accmean(struct acc *a)
{
   return(a->sum/(double)(a->n));
}

double variance(struct acc *a)
{
   double num, denom;
   if (a->n == 0) return((double)-infinity);
   num = a->sumsq - (a->sum)*(a->sum)/((double)(a->n));
   denom = (double)(a->n) - (double)1;
   return(num/denom);
}

/* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! normal
! ------
! Returns 2 random long reals from the standard normal distribution.
!
! Parameters: mu = mean
!           sdev = standard deviation
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*/

double normal(double mu, double sdev)
{
  double u1, u2, r;
  u1= uniform();
  if (u1==0.0) u1 = 0.00001;            /*prevent fatal error*/
  if (u1==1.0) u1 = .999999;
  u2= uniform();
  r = sqrt (-2.0*log(u1)) * cos(2.0*pi*u2);
  return(r*sdev + mu);
}

/*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! quick normal
! ------------
! Returns two uniform long reals from the normal distribution.
!
! Parameters: x1, x2 -> return values
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*/
quicknormal (double *x1, double *x2)
{
  double u1, u2, r;
  u1= uniform();
  if (u1==0.0) u1 = 0.00001;            /*prevent fatal error*/
  if (u1==1.0) u1 = .999999;
  u2 = uniform();
  *x1 = sqrt (-2.0*log(u1)) * cos(2.0*pi*u2);
  *x2 = sqrt (-2.0*log(u1)) * sin(2.0*pi*u2);
}

/* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! wishart
! -------
!
! Returns two uniform real numbers from the wishart (bivariate) distribution.
! These are obtained by squaring correlated uniform numbers from the bivariate
! normal distribution.
!
! Parameters : z1 = return value for symmetrical distribution (i.e. metric)
!              z2 =  ,,     ,,    ,, negative sided distribution (i.e. fitness)
!              epsilon1 = sqrt(E(z1*z1))
!              epsilon2 = sqrt(E(z2*z2))
!              rho = correlation of absolute values.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
wishart(double *z1, double *z2, double epsilon1, double epsilon2, double rho)
{
   double sq;
   rho = sqrt(rho);
   sq = 1/sqrt(3.0);
   quicknormal(z1, z2);
   *z2 = *z1*rho + *z2*sqrt(1.0 - rho*rho);
   *z1 = sq*epsilon1*(*z1)*(*z1);
   *z2 = -sq*epsilon2*(*z2)*(*z2);
   if (uniform() > 0.5) *z1 = -(*z1);
}

/* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! bivnormal
! -------
!
! Returns two real numbers from the bivariate normal distribution.
!
! Parameters : z1 = return value for symmetrical distribution (i.e. metric)
!              z2 =  ,,     ,,    ,, negative sided distribution (i.e. fitness)
!              epsilon1 = sqrt(E(z1*z1))
!              epsilon2 = sqrt(E(z2*z2))
!              rho = correlation of absolute values.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
bivnormal(double *z1, double *z2, double epsilon1, double epsilon2, double rho)
{
   rho = sqrt(rho);
   quicknormal(z1, z2);
   *z2 = *z1*rho + *z2*sqrt(1.0 - rho*rho);
   *z1 = epsilon1*(*z1);
   *z2 = -epsilon2*(*z2);
   if (*z2 > 0.0) *z2 = -*z2;
}

cubenormal(z1, z2, epsilon1, epsilon2, rho)
double *z1, *z2, epsilon1, epsilon2, rho;
{
   double sq;
   sq = 1.0/sqrt(16.0);
   rho = sqrt(rho);
   quicknormal(z1, z2);
   *z2 = *z1*rho + *z2*sqrt(1.0 - rho*rho);
   *z1 = sq*epsilon1*(*z1)*(*z1)*(*z1);
   *z2 = -sq*epsilon2*(*z2)*(*z2)*(*z2);
   if (*z2 > 0.0) *z2 = -*z2;
}

/* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! doublegamma
! -----------
!
! Returns uniform long real from the double gamma distribution.
!
! Parameter: epsilon = sqrt(E(a*a))
!            proportion positive = proportion of distribution +ve
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*/
double doublegamma (epsilon, proppositive)
double epsilon, proppositive;
{
  double r, sigma;
  sigma = sqrt(epsilon*1.0/sqrt(3.0));
  r= normal(0.0, sigma);
  r= r*r;
  if (uniform() > proppositive) r= -r;       /*uniformly allocate sign*/
  return(r);
}

double gamdev(int ia, double epsilon)
/* returns random number from gamma distribution with integer parameter
taken from Numerical Recipes in C */
{
   int j;
   double x;
   if ((ia<1)||(ia>=6)) gabort("Parameter to gamdev() invalid ", (double)ia);
   x = 1.0;
   for (j=1; j<=ia; j++)
   {
      x *= uniform();
   }
   x = -log(x);
   x = (epsilon*x)/sqrt((double)(ia*(ia+1)));
   return x;
}

initcovacc(struct covacc *a)
{
   a->n = 0;
   a->sumx = 0;
   a->sumy = 0;
   a->sumxy = 0;
   a->sumx2 = 0;
   a->sumy2 = 0;
}

covaccum(struct covacc *a, double x, double y)

{
   a->n = a->n + 1;
   a->sumx = a->sumx + x;
   a->sumy = a->sumy + y;
   a->sumxy = a->sumxy + x*y;
   a->sumx2 = a->sumx2 + x*x;
   a->sumy2 = a->sumy2 + y*y;
}

double covariance(struct covacc *a)
{
   if (a->n < 2)
   {
      return(-infinity);
   }
   return((a->sumxy - (a->sumx*a->sumy/(double)a->n))/((double)a->n - 1.0));
}

double correlation(struct covacc *a)
{
   double num, denom, cov, v1, v2;
   if (a->n < 2)
   {
      return(-infinity);
   }
   cov = covariance(a);
   num = a->sumx2 - (a->sumx)*(a->sumx)/((double)(a->n));
   denom = (double)(a->n) - 1.0;
   v1 = num/denom;
   num = a->sumy2 - (a->sumy)*(a->sumy)/((double)(a->n));
   denom = (double)(a->n) - 1.0;
   v2 = num/denom;
   return(cov/sqrt(v1*v2));
}

double se(a)
struct acc *a;
{
   if (a->n < 1) return(-infinity);
   if (variance(a)/a->n < 0) return(-infinity);
   return(sqrt(variance(a)/a->n));
}

printmse(s, r)
char *s;
struct acc *r;
{

   printf("%s %f +/- %f\n", s, accmean(r), se(r));
}

/*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! generate poisson table
! ----------------------
!
! generates a table of probabilities for each whole number given a poisson
! distribution.
!
! Parametes: mean = mean of poisson.
!            last number = max. elements in table
!            table -> table of probabilities.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*/

generatepoissontable (double mean, int *lastnumber, double *table, int max)
{
   int j;
   double prev;
   prev = exp(-mean);
   if (prev == 0) gabort("Insufficient resolution in generate poisson table", 0);
   table[0] = prev;
   for (j = 1; j <= max; j++)
   {
      *lastnumber = j;
      prev = prev*mean/(double)j;
      table[j] = table[j-1] + prev;
      if (1 - table[j] < 0.000001) break;
   }
}

/*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! poisson
! -------
!
! Looks up a table of probabilities of the poisson distribution already
! set up.  Returns the index of this probability.  This discrete vaue
! will come from the posson distribution.
!
! Parameters: last number = last entry in the table
!             table      -> table of probabilities
!
! Returns   : integer from poisson distribution.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*/

int poisson (int lastnumber, double *table)
{
   int j;
   float u;
   u = uniform();
   for (j = 0; j<=lastnumber; j++)
   {
      if ((double)u < table[j]) return(j);
   }
   return(lastnumber);
}

/*-------------------------------------------------------------------*/
/* binarysearchcumulative                                            */
/* ----------------------                                            */

/* Assumes array contains a cumulative density function of a discrete*/
/* distribution.  Returns an                                         */
/* integer randomly drawn from the density function.                 */

/* Parameters array -> cdf.                                          */
/*            size = no. of elements in array [1..size]              */
/*-------------------------------------------------------------------*/

int binarysearchcumulative(array, size)
double *array;
int size;
{
   int cur, pre, post;
   double r;
   r = uniform();
/*   printf("Uniform %f\n", r);*/
   pre = 1;
   post = size;
   cur = (size + 1)/2;
   do
   {
      if (array[cur] < r) pre = cur; else post = cur;
      cur = (pre + post)/2;
      if ((cur==post) || (cur==pre))
      {
         if (r < array[pre])
         {
/*            printf("Returning pre %d\n", pre);*/
            return(pre);
         }
         else
         {
/*            printf("Returning post %d\n", post);*/
            return(post);
         }
      }
   } while (size > 0);
}

/*----------------------------------------------------------------------*/
/*									*/
/* discrete 								*/
/* --------								*/
/*									*/
/* Returns random integer in range 1..n with equal probability.          */
/* Parameter: n = max. integer to return.				*/
/*									*/
/*----------------------------------------------------------------------*/
int discrete(int n)
{
   int res;
   res = (int)(uniform()*(double)n) + 1;
   if (res<1) res = 1;
   else if (res>n) res = n;
   return(res);
}

/*----------------------------------------------------------------------*/
/*                            						*/
/* samplewithoutreplacement						*/
/* ------------------------						*/
/*									*/
/* Returns a random integer from those present in array, which is valid */
/* from 1 to limit.  Overwrites sampled integer with array[limit],      */
/* then decrements limit.						*/
/* Parameters: limit -> no. valid integers in array                     */
/*             array -> array[1..limit] of integers to sample           */
/* 									*/
/*----------------------------------------------------------------------*/
int samplewithoutreplacement(limit, array)
int *limit, *array;
{
   int index, res;
   index = discrete(*limit);
   res = array[index];
   array[index] = array[*limit];
   *limit = *limit - 1;
   return(res);
}

/* trap - asks for input to stop program */
trap()
{
   int i;
   printf("Enter any int to continue\n");
   scanf("%d", &i);
}

/*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! getint;
! -------
!
! Displays a prompt and reads in an integer.  Checks the range of the integer
! and stops the program if it is out of range.
!
! Parameters: s -> string prompt
!             i -> integer variable to input.
!           min =  max value
!           max = min value.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
*/
getint(char *s, int *i, int min, int max)
{
   printf("%s", s);
   scanf("%d", i);
   if ((min != -infinity) && (*i < min))
   {
      printf("Value too small. Reading %s. Terminating.", s);
      printf("\n");
      printf("Min allowable: ");
      printf("%d", min);
      printf("\n");
      gabort("", 0);
   }
   if ((max != infinity) && (*i > max))
   {
      printf("Value too large. Reading %s. Terminating.", s);
      printf("\n");
      printf("Max allowable: ");
      printf("%d", max);
      printf("\n");
      gabort("", 0);
   }
}

/*
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! getint and skip;
! ----------------
!
! Displays a prompt and reads in an integer, skipping
! to the end of the line.  Checks the range of the integer
! and stops the program if it is out of range.
!
! Parameters: s -> string prompt
!             i -> integer variable to input.
!           min =  max value
!           max = min value.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
*/
getintandskip(char *s, int *i, int min, int max)
{
   char ch;
   printf(s);
   printf(" ");
   scanf("%d", i);
   do
   {
     ch = getchar();
   }
   while (ch != '\n');
   if ((min != -infinity) && (*i < min))
   {
      printf("Value too small. Reading %s. Terminating.", s);
      printf("\n");
      printf("Min allowable: ");
      printf("%d", min);
      printf("\n");
      gabort("", 0);
   }
   if ((max != infinity) && (*i > max))
   {
      printf("Value too large. Reading %s. Terminating.", s);
      printf("\n");
      printf("Max allowable: ");
      printf("%d", max);
      printf("\n");
      gabort("", 0);
   }
}

/*
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
*/
printreal(s,r,  x, y)
char *s;
double r;
int x, y;
{
   printf("%s %f\n", s, r);

}

/*
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! printint;
! ---------
!
! Prints out the integer with leading spaces x after the string s.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
*/
printint (s, i, x)
char *s;
int i, x;
{
   printf("%s %d", s, i);
   printf("\n");
}

/*
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! print tail
! ----------
!
! Prints a line of dashes.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
*/

printtail()
{
   int j;
   printf("\n");
   for (j=1; j<=80; j++)
   {
      printf("-");
   }
   printf("\n\n");
}

/*
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! cointoss
! --------
!
! Uses random number to return false with probability 0.5 otherwise
! returns a non zero integer.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
*/
int cointoss()
{
   if (uniform() > 0.5) return(0);
   return(1);
}

/*
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! printline
! ----------
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
*/
printline (s)
char *s;
{
   printf(s);
   printf("\n");
}

spaces(n)
int n;
{
   int i;
   for (i=1; i<=n; i++) printf(" ");
}

/*
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! getreal and skip;
! ----------------
!
! Displays a prompt and reads in an real, skipping
! to the end of the line.  Checks the range of the integer
! and stops the program if it is out of range.
!
! Parameters: s -> string prompt
!             i -> integer variable to input.
!           min =  max value
!           max = min value.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
*/
getrealandskip(char *s, double *i, double min, double max)
{
   char ch;
   printf(s);
   printf(" ");
   scanf("%lf", i);
   do
   {
     ch = getchar();
   }
   while (ch != '\n');
   if ((min != -infinity) && (*i < min))
   {
      printf("Value too small. Reading %s. Terminating.", s);
      printf("\n");
      printf("Min allowable: ");
      printf("%f", min);
      printf("\n");
      gabort("", 0);
   }
   if ((max != infinity) && (*i > max))
   {
      printf("Value too large. Reading %s. Terminating.", s);
      printf("\n");
      printf("Max allowable: ");
      printf("%f", max);
      printf("\n");
      gabort("", 0);
   }
}

double calculaterealmean(double *a, int n)
{
   double sum;
   int i;
   sum = 0.0;
   for (i=1; i<=n; i++)
   {
      sum = sum + a[i];
   }
   return(sum/(double)n);
}

tracestart()
{
printf("Enter trace level ");
scanf("%d", &tracelevel);
}

trace(s, i)
{
   if (tracelevel==0) return(0);
   printf("%s %d\n", s, i);
}

outputrep(int j, int howoften)
{
  if ((double)j/(double)howoften  - floor((double)j/(double)howoften) == 0.0)
     printf("Completed iteration %d\n", j);
}

outputrepdot(int j, int reps)
{
  int howoften;
  howoften = reps/10;
  if ((double)j/(double)howoften  - floor((double)j/(double)howoften) == 0.0)
  {
     printf("%d.", j/howoften);
     if (j==reps) printf("\n");
     fflush(stdout);
  }
}

monitorinput()
{
   inputcounter--;
   if (inputcounter <= 0)
   {
      printf("Enter int to continue ");
      scanf("%d", &inputcounter);
   }
}

double normalheight(double x, double mu, double var)
{
   double res;
   res  = (1.0/sqrt(2.0*pi*var))*exp(-(x-mu)*(x-mu)/(2.0*var));
   return(res);
}

/*
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! simpson
! -------
!
! Computes the area under the function using simpsons rule for the set of
! points given.
!
! Parameters: x -> array of pairs of points.
!              points = no. of points.
!                   a = lower value
!                   v = upper value
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
*/

double simpson (double *x, int points, double a, double b)
{
   int j, coeff;
   double  res, h;
   if (2*(points/2) == points)
   {
      printf("FATAL ERROR: Even no. of points for SIMPSON.\n");
      monitorinput();
   }
   h = (b-a) / ((double)points-1.0);
   res = 0.0;
   for (j = 1; j<=points; j++)
   {
      if ((j == points) || (j==1 ))
      {
         coeff = 1;
      }
      else if (j == 2)
      {
         coeff = 4;
      }
      else if (coeff == 2) coeff = 4; else coeff = 2;
/*      %if trace level > 0 %then printstring("Coeff, x(j) ") %and write(coeff, 3) %and print(x(j), 2, 4) %and newline*/
      res = res + (double)coeff*x[j];
   }
   return((h/3.0) * res);
}

int odd(int i)
{
   int x;
   x = i/2;
   if ((double)(i)/2.0 - (double)x == 0.0) return(false); else return(true);
}

/*
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! solve quadratic
! ---------------
!
! Returns the roots of quadratic of parameters given.
!
! Parameters: a, b, c = parameters of quadratic.
!             root1, root2 = roots of quadratic.
!
! Returns : true => real roots.
!          false => no real root.
!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
*/
int solvequadratic(double a, double b, double c, double *r1, double *r2)
{
   double x;
   x = b*b - 4.0*a*c;
   if (x < 0.0) return(false);
   *r1 = (-b + sqrt(x))/(2.0*a);
   *r2 = (-b - sqrt(x))/(2.0*a);
   return(true);
}

int skiptoendofline(FILE *fptr)
{
   int status;
   char ch;
   do
   {
      status = fscanf(fptr, "%c", &ch);
      if (status==EOF) break;
/*      printf("skiptoendofline ch %c\n", ch);*/
   }
   while (ch!='\n');
   return(status);
}

/*
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
!
! aitken
! ------
!
! Uses Aitken's formula to predict the asymptote from the three last points
! in the array given.
!
! Parameters: x-> array of points.
!               t = max. point in array.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
*/
double aitken (double *a, int t)
{
   double x;
   x = a[t-2] - 2.0*a[t-1] + a[t];
   if (x == 0) return(a[t]);
   return(-((a[t-1] - a[t])*(a[t-1] - a[t]))/x + a[t]);
}

int factorial(int i)
{
   int res;
   res = 1;
   for (;;)
   {
      if (i==0) break;
      res = res*i;
      i--;
   }
   return(res);
}

double doublefactorial(double x)
{
   double res;
   x = floor(x);
   res = 1.0;
   for (;;)
   {
      if (x==0.0) break;
      res = res*x;
      x = x - 1.0;
   }
   return(res);
}

double logdoublefactorial(double x)
{
   double res;
   x = floor(x);
   res = 0.0;
   for (;;)
   {
      if (x==0.0) break;
      res = res + log(x);
      x = x - 1.0;
   }
   return(res);
}

/*
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! gammaalpha
! ----------
!
! Returns the value of the parameter alpha of the gamma function.
!
!
! Parameters: beta = parameter of gamma func.

!            epsilon = sqrt(E(a*a))
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
*/
double gammaalpha( beta, epsilon)
double beta, epsilon;
{
   double res;
   res = sqrt(beta*(beta + 1.0))/epsilon;
   return(res);
}

double asymptote_gamma(double z)
{
   double res;
   res = (z - 0.5)*log(z) -z +0.5*log(2.0*pi) + 1.0/(12.0*z)
        - 1.0/(360.0*z*z*z) + 1.0/(1260.0*z*z*z*z*z) - 1.0/(1680.0*z*z*z*z*z*z*z);
   return(exp(res));
}

double series_gamma(double z)
/* return gamma(z) using series expansion 6.1.34 of Abramowitz and Stegun */
{
   #define maxc 26
   double c[maxc+1], res;
   int i;
   if (z > 2.0) return(asymptote_gamma(z));
   c[1] = 1.0;
   c[2] = 0.5772156649015329;
   c[3] =-0.6558780715202538;
   c[4] =-0.0420026350340952;
   c[5] = 0.1665386113822915;
   c[6] =-0.0421977345555443;
   c[7] =-0.0096219715278770;
   c[8] = 0.0072189432466630;
   c[9] =-0.0011651675918591;
   c[10]=-0.0002152416741149;
   c[11]= 0.0001280502823882;
   c[12]=-0.0000201348547807;
   c[13]=-0.0000012504934821;
   c[14]= 0.0000011330272320;
   c[15]=-0.0000002056338417;
   c[16]= 0.0000000061160950;
   c[17]= 0.0000000050020075;
   c[18]=-0.0000000011812746;
   c[19]= 0.0000000001043427;
   c[20]= 0.0000000000077823;
   c[21]=-0.0000000000036968;
   c[22]= 0.0000000000005100;
   c[23]=-0.0000000000000206;
   c[24]=-0.0000000000000054;
   c[25]= 0.0000000000000014;
   c[26]= 0.0000000000000001;
   res = 0.0;
   for (i=1; i<=26; i++)
   {
      res += c[i]*pow(z, (double)i);
/*      printf("gamma(z) %15.12f\n", 1.0/res); */
   }
   return(1.0/res);
}

/*
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! gammaht
! -------
!
! Returns the height at point a of a gamma distribution.
!
! Parameters: epsilon, beta = gamma parameters.
!             a = value to evaluate.
!
! Returns gamma fn.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
*/
double gammaht(double epsilon, double beta, double a)
{
   double gammabeta, res, a1, a2, a3, alpha;
   if (a==0.0)
   {
      printf("Illegal value (0.0) for gammaht()\n");
      monitorinput();
   }
   if (beta < 0.0) gabort("gammaht: invalid beta\n", beta);
   alpha = gammaalpha(beta, epsilon);
   gammabeta = series_gamma(beta);
   a1 = pow(alpha, beta);
   a2 = exp(-alpha*a);
   a3 = pow(a, beta - 1.0);
   res = a1*a2*a3/gammabeta;
   return(res);
}

/*
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! half points
! -----------
!
! This routine is called to check that the simpson numerical integration
! is converging satisfactorily.  It halves the number of points in
! the array specified.
!

! Parameters: a -> array of points.
!              points -> no. of points.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
*/
halfpoints(double *a, int *points)
{
   int ptr1, ptr2, newpoints;
   if (odd(*points)==false)
   {
      printf("Half points: ERROR: no. of points for simpson must be odd.35\n");
      monitorinput();
   }
   ptr1 = 1;
   ptr2 = 1;
   do
   {
      a[ptr1] = a[ptr2];
      newpoints = ptr1;
      ptr1++;
      ptr2 = ptr2 + 2;
   }
   while (ptr2 <= *points);
   *points = newpoints;
}

setuppoints(int pts, double *points, double lower, double upper, fun1 fn)
{
      int i;
      double interval, x, a;
      interval = (upper - lower)/((double)pts - 1.0);
      x = lower;
      for (i = 1; i<=pts; i++)
      {
         a = fn(x);
         points[i] = a;
         x = x + interval;
      }

}

/*
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
!
! numerical integ
! ---------------
!
! Numerically integrates the function fn and returns the result in res.
! Integrates over a number of ranges with different numbers of points per
! range.  Points must be a number like 65, 129, 257, 513... in each
! range.  Halves number of points to do integration 3 times as a check on
! convergence
!
! Parameters:
!
!
! resvec: contains vector of result for different numbers of points
! start, finish: values limiting the function
! fn(x): is the function to integrate
! ranges: no. of ranges of points to evaluate.
! points per range: no of points in each range.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
*/

double numericalinteg(double *resvec, double start, double finish, fun1 fn, int ranges,
   int *pointsperrange)
{
   double lower, upper, r;
   double res[maxranges+2][4];
/*        The first dimension is contains the integral for each range, the
       last element containing the total.  The second dimension contains
       results for each halving of the points.
*/
   double points[maxsimpsonpoints+1];
   int i, j, k, fac, p;
   static int tmoutfileflag = 0;
   if (tmoutfileflag==0)
   {
/*      tmoutfile = fopen("numericalinteg.out", "w");*/
      tmoutfileflag = 1;
   }
/*   printf("Start %f, finish %f, ranges %d\n", start, finish, ranges);*/
   upper = start;
   for (i = 1;  i<= ranges; i++)
   {
      lower = upper;
      upper = lower + (finish - start)/(double)ranges;
      p = pointsperrange[i];
      if (p > maxsimpsonpoints) gabort("too many points for simpson",p);
      setuppoints(p, points, lower, upper, fn);
      for (j = 1; j<=3; j++)
      {
         r = simpson(points, p, lower, upper);
         res[i][j] = r;
         if (j != 3) halfpoints(points, &p);
      }
   }
   for (i = 1; i<=3; i++)
   {
      res[ranges+1][i] = 0.0;
      for (j = 1; j<=ranges; j++)
         res[ranges+1][i] = res[ranges+1][i] + res[j][i];
      resvec[i] = res[ranges+1][i];
   }
/*   fprintf(tmoutfile, "Results of numerical integration\n");
   for (i=1; i<=ranges; i++) fprintf(tmoutfile, " Pts   Range %d", i);
   fprintf(tmoutfile, "     All");
   fprintf(tmoutfile, "\n");
*/
   fac = 1;
/*   for (i = 1; i<=3; i++)
   {
      for (j = 1; j<=ranges+1; j++)
      {
         if (j!=ranges+1) fprintf(tmoutfile, "%4d", (pointsperrange[j]-1)/fac + 1);
         fprintf(tmoutfile, "%10.4f", res[j][i]);
      }
      fac = fac*2;
      fprintf(tmoutfile, "\n");
   }
*/
   return(res[ranges+1][1]);
}

FILE *openforread(char *str)
{
   FILE *f;
   f = fopen(str, "r");
   if (f==0)
   {
      printf("ERROR: File %s not found.\n", str);
      exit(0);
   }
   else //printf("Opened file %s for read.\n", str);
   return(f);
}

FILE *openforwrite(char *str, char *mode)
{
   FILE *f;
   f = fopen(str, mode);
   if (f==0)
   {
      printf("ERROR: Unable to open file %s for write.\n", str);
      exit(0);
   }
   else //printf("Opened file %s for write, ", str);
   if (mode[0]=='a') printf("append mode.\n");
   else if (mode[0]=='w') printf("overwrite mode.\n");
   else gabort("Invalid mode given in openforwrite\n", 0);
   return(f);
}

FILE *openforreadsilent(char *str)
{
   FILE *f;
   f = fopen(str, "r");
   if (f==0)
   {
      printf("ERROR: File %s not found.\n", str);
      exit(0);
   }
   return(f);
}

FILE *openforwritesilent(char *str, char *mode)
{
   FILE *f;
   f = fopen(str, mode);
   if (f==0)
   {
      printf("ERROR: Unable to open file %s for write.\n", str);
      exit(0);
   }
   if (mode[0]=='a')
   {
   }
   else if (mode[0]=='w')
   {
   }
   else gabort("Invalid mode given in openforwrite\n", 0);
   return(f);
}

printmsefile(f, s, r)
FILE *f;
char *s;
struct acc *r;
{

   fprintf(f, "%s %f +/- %f\n", s, accmean(r), se(r));
}

quadratic_regression(double n, double sumx, double sumx2, double sumy, double sumx3,
  double sumxy, double sumx4, double sumx2y, double *b1, double *b2, double *b3)
{
/* Solve system of 3 simultaneous equations:

a*b1 + b*b2 + c*b3 = d         (1)
e*b1 + f*b2 + g*b3 = h         (2)
i*b1 + j*b2 + k*b3 = l         (3)

The actual parameters are appropriate for solving a quadratic regression -

y = b3*x^2 + b2*x + b3

See assignments below.

*/
   double a, b, c, d, e, f, g, h, i, j, k, l;
   double z1, z2, z3, z4, z5, z6;
   a = n;
   b = sumx;
   c = sumx2;
   d = sumy;
   e = b;
   f = c;
   g = sumx3;
   h = sumxy;
   i = c;
   j = g;
   k = sumx4;
   l = sumx2y;
/*   printf("\n%lf %lf %lf %lf\n", a, b, c, d);
   printf("%lf %lf %lf %lf\n", e, f, g, h);
   printf("%lf %lf %lf %lf\n\n", i, j, k, l);
*/

/* (1) - (2) -> b1*z1 + b2*z2 = z3  (4)     */
   z1 = a/c - e/g;
   z2 = b/c - f/g;
   z3 = d/c - h/g;

/* (1) - (3) -> b1*z4 + b2*z5 = z6  (5)     */
   z4 = a/c - i/k;
   z5 = b/c - j/k;
   z6 = d/c - l/k;

/* (4) - (5) */

   *b1 = (z3/z2 - z6/z5)/(z1/z2 - z4/z5);

   *b3 = ((d - a*(*b1))/b - (h - e**b1)/f) / (c/b - g/f);

   *b2 = (d - c*(*b3) - a*(*b1))/b;
}

int findtext(char *s, FILE *inptr)
{
   int len, i, curchar, dum;
   char c;
   curchar = 0;
   len = 0;
   for (i=0; i<=100; i++)
   {
      if (s[i] == '\0') break;
      len++;
   }
   for (;;)
   {
      c = getc(inptr);
      if (c == EOF) return(false);
      if (c==s[curchar])
      {
         curchar++;
         if (curchar==len) return(true);
      }
      else curchar = 0;
   }
}

/* FIND1DMAX */
/* Finds the maximum of 3 points using a quadratic approximation */

find1dmax(double *xvec, double *yvec, double *xmax, double *ymax)
{
   int i;
   double a, b, c, x1, x2, x3, y1, y2, y3, temp;
   double num, denom;
/*   for (i=1; i<=3; i++)
   {
      printf("xvec[%1d] %lf yvec[%1d] %lf\n", i, xvec[i], i, yvec[i]);
   }
*/
   x1 = xvec[1];
   x2 = xvec[2];
   x3 = xvec[3];
   y1 = yvec[1];
   y2 = yvec[2];
   y3 = yvec[3];
   if (y1==y3) y1+= 0.000000000000001;
   if (y3==y2) y2+= 0.0000000000000001;
/*   printf("x1 %f x2 %f x3 %f y1 %f y2 %f y3 %f\n", x1, x2, x3, y1, y2, y3);*/
   num = (y1 - y3)*(x1 - x2)/(x1 - x3) - y1 + y2;
   denom = (x1*x1 - x3*x3)*(x1 - x2)/(x1 - x3) - (x1*x1 - x2*x2);
   a = num/denom;
   num = y1 - y2 - a*(x1*x1 - x2*x2);
   b = num/(x1 - x2);
   c = y1 - a*x1*x1 - b*x1;
   temp = -b/(2.0*a);
   *xmax = temp;
   *ymax = a*temp*temp + b*temp + c;
}

/* FIND2DMAX */
/* Finds the maximum of a 3x3 grid using a sequential quadratic approximation */

double find2dmax(double *xvec, double *yvec, double zmatrix[4][4], double *xmax, double *ymax)
{
   int i, j;
   double vec[4], vec1[4], temp, max, max1, max2;
/*   for (i=1; i<=3; i++)
   {
      for (j=1; j<=3; j++)
      {
         if (zmatrix[i][j]==undefined) return(undefined);
      }
   }
*/
/* find the maximum for columns */
   for (i=1; i<=3; i++)
   {
      for (j=1; j<=3; j++)
      {
         vec[j] = zmatrix[i][j];
      }
      find1dmax(yvec, vec, &temp, &max);
/*      printf("Find1smax: zmax %lf ymax %lf\n", max, temp);*/
      vec1[i] = max;
   }
   find1dmax(xvec, vec1, xmax, &max1);
/*   printf("Columns: max %lf xmax %lf\n", max1, *xmax);*/
/* find the maximum for rows */
   for (i=1; i<=3; i++)
   {
      for (j=1; j<=3; j++)
      {
         vec[j] = zmatrix[j][i];
      }
      find1dmax(xvec, vec, &temp, &max);
/*      printf("Find1dmax: zmax %lf ymax %lf\n", max, temp);*/
      vec1[i] = max;
   }
   find1dmax(yvec, vec1, ymax, &max2);
/*   printf("Rows: max %lf ymax %lf\n", max2, *ymax);*/
   max = (max1 + max2)/2.0;
   return(max);
}

int bnldev(double pp, int n)
{
   int j, bnl;
   double p;
   p = (pp < 0.5 ? pp : 1.0 - pp);
   bnl = 0;
   for (j=1; j<=n; j++)
      if (uniform() < p) ++bnl;
   if (p!=pp) bnl = n - bnl;
   return bnl;
}

void metropolis(double y[], int ndim, double(*func)(double []),
   int limit, double delta[], double prevres, double *temperature)
{
   #define step_change 0.2
   int i, j, param, accept;
   double step, res, diff, u;
   param = 0;
   *temperature = 1.0;
   for (i=1; i<=limit; i++)
   {
/*      printf("%d ", i);*/
      param++;
      if (param > ndim) param = 1;
      step = normal(0.0, delta[param]);
/*      printf("param %d delta[param] %lf, step %lf\n", param, delta[param], step);*/
      y[param] += step;
      res = (*func)(y);
      if (res > prevres)
      {
         diff = prevres - res;
         u = uniform();
/*         printf("diff %lf, exp(diff) %lf, u %lf\n", diff, exp(diff), u);*/
         if (u < pow(exp(diff), 1.0/(*temperature)))
         {
            accept = true;
/*            printf("Accepted\n");*/
         }
         else
         {
            accept = false;
/*            printf("Rejected\n");*/
         }
      }
      else
      {
         accept = true;
/*         printf("Accepted unconditionally\n");*/
      }
      if (accept == true)      /* Note assume to be minimizing not maximizing */
      {
         prevres = res;      /* accept the change */
         delta[param] *= (1.0 + step_change); /* Increase the step size */
      }
      else
      {
         y[param] -=step;     /* reject the change */
         delta[param] *= (1.0 - step_change); /* Decrease the step size */
      }
/*      printf("Dump of delta vector: ");
      for (j=1; j<=ndim; j++)
      {
         printf("%lf ", delta[j]);
      }
      printf("\n"); monitorinput();
*/
      *temperature -= 2.0/(double)limit;
      if (*temperature < 0.001) *temperature = 0.001;
   }
}

/* A slower way of generating poissons that the tabular method in poisson */
int genpoisson(double xm)
{
   static double sq, alxm, g, oldm = (-1.0);
   double em, t, y;
//   if (xm>300) gabort("genpoisson: xm too large", 0);
   if (xm>=200.0)   /* Use normal generator for very high xm */
   {
      return normal(0.0, sqrt(xm)) + xm;
   }
   if (xm!=oldm)
   {
      oldm = xm;
      g = exp(-xm);
   }
   em = -1;
   t = 1.0;
   do
   {
      ++em;
      t *= uniform();
   }
   while (t > g);
   return em;
}

void qsortreals (double *a, int from, int to)
   {
      int l, v;
      double d;
      if (from >= to) return;
      l = from;
      v = to;
      d = a[v];
      for (;;)
      {
         while ((l < v) && (a[l] <= d))
         {
            l = l + 1 ;
         }
         if (l==v) break;
         a[v] = a[l];
         while ((v > l) && (a[v] >= d))
         {
            v = v - 1;
         }
         if (v==l) break;
         a[l] = a[v];
      }

/*      ! now l=v*/
      a[v] = d;
      l = l-1;
      v = v+1;
      qsortreals(a, from, l);
      qsortreals(a, v, to);
}

#define maxofrs 100

FILE *openfilecheck2dir(char *str)
/*
 Attempts to open file str, if fails first time attempts to open
../str or ../../str
*/
{
   FILE *f;
   static char s[maxofrs];
   int i;
   f = fopen(str, "r");
   if (f==0)
   {
      s[0] = '.';
      s[1] = '.';
      s[2] = '/';
      i=0;
      for (;;)
      {
         if (i+3 > maxofrs) gabort("String too long", i);
         s[i+3] = str[i];
         if (str[i] == '\0') break;
         i++;
      }
      f = fopen(s, "r");
      if (f==0)
      {
         s[0] = '.';
         s[1] = '.';
         s[2] = '/';
         s[3] = '.';
         s[4] = '.';
         s[5] = '/';
         i=0;
         for (;;)
         {
            if (i+6 > maxofrs) gabort("String too long", i);
            s[i+6] = str[i];
            if (str[i] == '\0') break;
            i++;
         }
         f = fopen(s, "r");
         if (f==0)
         {
            printf("ERROR: File %s not found.\n", str);
            exit(0);
         }
         else  printf("Opened file %s for read.\n", s);
      }
      else  printf("Opened file %s for read.\n", s);
   }
   else printf("Opened file %s for read.\n", str);
   return(f);
}

float rtsafe(void (*funcd)(float, float *, float *), float x1, float x2, float xacc)
{
   #define MAXIT 100
   int j;
   float df, dx, dxold, f, fh, fl;
   float temp, xh, xl, rts;
   (*funcd)(x1, &fl, &df);
   (*funcd)(x2, &fh, &df);
   if ((fl > 0.0 && fh > 0.0) || (fl < 0.0 && fh < 0.0))
      gabort("Rtsafe: roots must be bracketed", 0);
   if (fl == 0.0) return(x1);
   if (fh == 0.0) return(x2);
   if (fl < 0.0)
   {
      xl = x1;
      xh = x2;
   }
   else
   {
      xh = x1;
      xl = x2;
   }
   rts = 0.5*(x1+x2);
   dxold = fabs(x2-x1);
   dx = dxold;
   (*funcd)(rts, &f, &df);
   for (j=1; j<=MAXIT; j++)
   {
      if ((((rts-xh)*df-f)*((rts-xl)*df-f) >=0.0)
        || (fabs(2.0*f) > fabs(dxold*df)))
      {
         dxold = dx;
         dx = 0.5*(xh-xl);
         rts = xl +dx;
         if (xl==rts) return rts;
      }
      else
      {
         dxold = dx;
         dx = f/df;
         temp = rts;
         rts -= dx;
         if (temp==rts) return (rts);
      }
      if (fabs(dx) < xacc) return rts;
      (*funcd)(rts, &f, &df);
      if (f < 0.0)
         xl = rts;
      else
         xh = rts;
   }
   printf("warning: rtsafe: Maximum interations exceeded\n"); monitorinput();
   return (rts);
}

/* Algorith to generate bivariate normal deviates provided by Ian
   White
*/

bivariatenormal(double *z1, double *z2, double var1, double var2, double cov)
{
   double s;
   quicknormal(z1, z2);     /* uncorrelated normal deviates mean 0 SD 1 */
   *z1 *= sqrt(var1);       /* Scale Z1 according to the standard deviation */
   s = var2 - cov*cov/var1;
   if (s<0) gabort("Invalid parameters for bivnormal", s);
   *z2 = (cov/var1)*(*z1) + sqrt(s)*(*z2);
}

get_silent_mode(int argc, char *argv[], int *silent_mode)
{
   int n_command_line_param, i;
   *silent_mode = 0;
   n_command_line_param = argc;
//   printf("n_command_line_param %d\n", n_command_line_param);
   for (i=1; i<n_command_line_param; i++)
   {
//      printf("argv[%d] %s\n", i, argv[i]);
   }
   if (n_command_line_param==2)
   {
      if (strcmp(argv[1], "-s")==0)
      {
//         printf("Silent mode switched on\n");
         *silent_mode = 1;
      }
   }
}

make_numbered_file_name(char *outfilename, char *filestr, int ind)
{
   int len, i;
   int i1;
   char ch;
   len = strlen(filestr);
   for (i=0; i<len; i++) outfilename[i] = filestr[i];
   if (ind > 99)
      gabort("makefilename: file index too large", ind);
   if (ind > 9)
   {
      i1 = ind/10;
      ch = i1 + '0';
      outfilename[len] = ch;
//      printf("Char 1 ind %d i1 %d ch %c\n", ind, i1, ch); monitorinput();
      i1 = ind - (ind/10)*10;
      ch = i1 + '0';
//      printf("Char 2 ind %d i1 %d ch %c\n", ind, i1, ch); monitorinput();
      outfilename[len+1] = ch;
      outfilename[len+2] = '\0';
   }
   else
   {
      ch = ind + '0';
//      printf("Char 1 ind %d i1 %d ch %c\n", ind, i1, ch); monitorinput();
      outfilename[len] = ch;
      outfilename[len+1] = '\0';
   }
//   printf("makefilename: results %s\n", outfilename);
}

// --------------------------------------------------------------------
// Routine for grid searching a function f that takes integer values as
// it's paramters
// --------------------------------------------------------------------

double fill_out_grid(double **grid, int x0, int x1, int y0, int y1,
   int step_size_x, int step_size_y,
   double (*f)(int, int), int *max_x, int *max_y, int *nfunc)
{
   double max_fun, f_result;
   int i, j, evaluated;
   *max_x = -1;
   *max_y = -1;
   max_fun = undefined;
   i = x0;
   for (;;)
   {
      j = y0;
      for (;;)
      {
         evaluated = 0;
         if (grid[i][j]==undefined)
         {
            f_result = (*f)(i, j);
            (*nfunc)++;
            grid[i][j] = f_result;
            evaluated = 1;
         }
         else f_result = grid[i][j];
//         printf("i %d j %d f_result %lf ", i, j, f_result);
//         if (evaluated == 1) printf("evaluated\n");
//         else printf("not evaluated\n");
         if ((f_result > max_fun)||(max_fun==undefined))
         {
            max_fun = f_result;
            *max_x = i;
            *max_y = j;
         }
         if (j==y1) break;
         j += step_size_y;
         if (j > y1) j = y1;
      }
      if (i==x1) break;
      i += step_size_x;
      if (i > x1) i = x1;
   }
   return max_fun;
}

int_to_string(int v, char *str, int max_len)
{
   int divisor = 10000000, i1, ind = 0, digit_encountered = 0;
   char ch;
   if (v > divisor) gabort("int_to_string: integer parameter to large", v);
   for (;;)
   {
      i1 = v/divisor;
//      printf("i1 %d divisor %d\n", i1, divisor);
      if (i1 >= 0)
      {
         ch = i1 + '0';
         if (ch!='0') digit_encountered = 1;
//         printf("ch %c\n", ch);
         if (digit_encountered)
         {
            str[ind] = ch;
            ind++;
            if (ind == max_len) gabort("int_to_string: string too long", ind);
         }
         v = v - i1*divisor;
      }
      divisor /= 10;
      if (divisor == 1) break;
//      monitorinput();
   }
   ch = v + '0';
//   printf("ch %c\n", ch);
   str[ind] = ch;
   ind++;
   if (ind == max_len) gabort("int_to_string: string too long", ind);
   str[ind] = '\0';
}

// bilinear_interpolation:
// First 4 parameters are positions on x and y axes
// Second 4 are function values
// x and y are positions to be evaluated.

double bilinear_interpolation(double x1, double x2, double y1, double y2,
   double x1y1, double x1y2, double x2y1, double x2y2,
   double x, double y)
{
   double res;
   res = (x2 - x)*(y2 - y)*x1y1 +
         (x - x1)*(y2 - y)*x2y1 +
         (x2 - x)*(y - y1)*x1y2 +
         (x - x1)*(y - y1)*x2y2;

   res /= ((x2 - x1)*(y2 - y1));
   return res;
}

// --------------------
// kimura_fixation_prob
// --------------------
// Returns fixation probability of an additive mutation, assuming that s
// is small in a diploid population
//
// Parameters
// ---------
// s - selection coefficient, difference between the homozygotes.
// N - Population effective size.
//

double kimura_fixation_prob(double s, int N)
{
   double num, denom;
   if ((s < 0.0)&&((double)N*s > -0.0000001))  // Trap numerical problems
   {
      num = 1.0;
      denom = 2.0*(double)N;
   }
   else if ((s > 0.0)&&((double)N*s < 0.0000001))
   {
      num = 1.0;
      denom = 2.0*(double)N;
   }
   else
   {
      num = 1 - exp(-s);
      denom = 1 - exp(-2*(double)N*s);
   }
//   printf("s %14.14lf N %d num %lf denom %lf num/denom %lf\n",
///   s, N, num, denom, num/denom);
//   monitorinput();
   return num/denom;
}

// --------------------
// kimura_fixation_prob
// --------------------
// Returns fixation probability of an additive mutation, assuming that s
// is small in a diploid population
//
// Parameters
// ---------
// s - selection coefficient, difference between the homozygotes.
// N - Population effective size.
//

double kimura_fixation_prob_X_loci(double s, int N)
{
   double num, denom;
   if ((s < 0.0)&&((double)N*s > -0.0000001))  // Trap numerical problems
   {
      num = 1.0;
      denom = 2.0*(double)N;
   }
   else if ((s > 0.0)&&((double)N*s < 0.0000001))
   {
      num = 1.0;
      denom = 2.0*(double)N;
   }
   else
   {
      num = 1 - exp(-s);
      denom = 1 - exp(-2*(double)N*s);
   }
//   printf("s %14.14lf N %d num %lf denom %lf num/denom %lf\n",
///   s, N, num, denom, num/denom);
//   monitorinput();
   return num/denom;
}

//----------------------------------------
// normal_cumulative_probability_function
// ---------------------------------------
//
// Returns the area under the normal pdf from -infinity to x.
// From algorithm described in Abramowitz and Stegun.

double normal_cumulative_probability_function(const double x)
{
  const double b1 =  0.319381530;
  const double b2 = -0.356563782;
  const double b3 =  1.781477937;
  const double b4 = -1.821255978;
  const double b5 =  1.330274429;
  const double p  =  0.2316419;
  const double c  =  0.39894228;

  if(x >= 0.0) {
      double t = 1.0 / ( 1.0 + p * x );
      return (1.0 - c * exp( -x * x / 2.0 ) * t *
      ( t *( t * ( t * ( t * b5 + b4 ) + b3 ) + b2 ) + b1 ));
  }
  else {
      double t = 1.0 / ( 1.0 - p * x );
      return ( c * exp( -x * x / 2.0 ) * t *
      ( t *( t * ( t * ( t * b5 + b4 ) + b3 ) + b2 ) + b1 ));
    }
}

//----------------------------------------
// dump_matrix
//----------------------------------------
//
// Prints out the matrix X with 1..rows, 1..cols rows and columns, respectively.
//
// s -> string to print
// X -> matrix, usually allocated using dmatrix(), dimension rows, cols
// rows = rows, 1 relative
// cols = cols, 1 relative

dump_matrix(char *s, double **X, int rows, int cols)
{
   int i, j;
   printf ("%s\n", s);
   for (i = 1; i <=rows; i++)
   {
      for (j = 1; j <= cols; j++)
      {
         printf("%lf ", X[i][j]);
      }
      printf("\n");
   }
}

//----------------------------------------
// compute_matrix_transpose
//----------------------------------------
//
// Generates the transpose of X.
//
// X -> matrix, usually allocated using dmatrix(), dimension rows, cols
// rows = rows of X, 1 relative
// cols = cols of X, 1 relative
// XT -> result matrix, allocated using dmatrix(), dimension cols, rows

compute_matrix_transpose(double **X, double **XT, int rows, int cols)
{
   int i, j;
   for (i=1; i<=rows; i++)
   {
      for (j=1; j<=cols; j++)
      {
         XT[j][i] = X[i][j];
      }
   }
}

//----------------------------------------
// multiply_matrices
//----------------------------------------
//
// Carry out matrix operation XY, putting result in RES.
//
// RES -> result matrix dimension rows_X, cols_X, usually allocated using
//          dmatrix()
// X -> matrix dimension rows_X, cols_X, usually allocated using dmatrix()
// Y -> matrix dimension rows_Y, cols_Y, usually allocated using dmatrix()
// rows_X = rows of X, 1 relative
// cols_X = cols of X, 1 relative
// rows_Y = rows of Y, 1 relative
// cols_Y = cols of Y, 1 relative

multiply_matrices(double **RES, double **X, double **Y, int rows_X, int cols_X,
   int rows_Y, int cols_Y)
// Matrix X has dimension rows_X, cols_X and Y has dimension rows_Y, cols_Y. The
// result has dimension rows_X, cols_Y.
{
   int i, j, k;
   double element;
   if (cols_X != rows_Y)
   {
      gabort("multiply_matrices: cols_X != rows_Y", 0);
   }
//   printf("multiply_matrices: rows_X %d cols_X %d\n", rows_X, cols_X);
//   dump_matrix("multiply_matrices: Matrix X", X, rows_X, cols_X);
//   monitorinput();
//   printf("multiply_matrices: rows_Y %d cols_Y %d\n", rows_Y, cols_Y);
//   dump_matrix("multiply_matrices: Matrix Y", Y, rows_Y, cols_Y);
//   monitorinput();
   for (i=1; i<=rows_X; i++)
   {
      for (j=1; j<=cols_Y; j++)
      {
         element = 0.0;
         for (k=1; k<=cols_X; k++)
         {
//            printf("X[%d][%d] %lf X[%d][%d] %lf\n",
//                 i, k, X[i][k], k, j, Y[k][j]);
            element += X[i][k]*Y[k][j];
//            printf("i %d j %d k %d element %lf\n", i, j, k, element);
//            monitorinput();
         }
//         printf("RES[%d][%d] %lf\n", i, j, element); monitorinput();
         RES[i][j] = element;
      }
   }
}

// ----------------------------------------------
// normal_approx_to_binomial
// ----------------------------------------------
//
// Uses the normal approximation to calculate the probability of i successes
// from n trials if the probabilituy of one sucess is p

double normal_approx_to_binomial(int i, int n, double p)
{
   double m, s, z_upper, z_lower, area_lower, area_upper;
   m = p*(double)n;
   s = sqrt((1.0-p)*p*(double)n);
//   printf("m %lf s %lf\n", m, s);
   z_lower = ((double)i - 0.5 - m)/s;
   area_lower = normal_cumulative_probability_function(z_lower);
   z_upper = ((double)i + 0.5 - m)/s;
   area_upper = normal_cumulative_probability_function(z_upper);
//   printf("z_lower %lf area_lower %lf z_upper %lf area_upper %lf\n",
//      z_lower, area_lower, z_upper, area_upper);
   return(area_upper - area_lower);
}

// ----------------------------------------------

// loadparam
// paramvec = vector with 2 elements: element 0 = 0 => parameter not variable
//                                    element 0 = 1 => parameter variable
//                                    element 1 location to get parameter value
// curpos = index - 1 from which to load parameter
// p: matrix into which parameter loaded
//

loadparam(double *paramvec, int *curpos, double **p)
{
   if (paramvec[0]!=0.0)
   {
      *curpos = *curpos + 1;
      p[1][*curpos] = paramvec[1];
   }
}

// ----------------------------------------------
// unloadparam
// paramvec = vector with 2 elements: element 0 = 0 => parameter not variable
//                                    element 0 = 1 => parameter variable
// curpos = index -1 in which to unload parameter
// x = vector containing parameter value

unloadparam(double *paramvec, int *curpos, double x[])
{
   if (paramvec[0]!=0.0)
   {
      *curpos = *curpos + 1;
      paramvec[1] = x[*curpos];
   }
}


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

#ifndef __GSL_BLOCK_UINT_H__
#define __GSL_BLOCK_UINT_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

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

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

struct gsl_block_uint_struct
{
  size_t size;
  unsigned int *data;
};

typedef struct gsl_block_uint_struct gsl_block_uint;

GSL_FUN gsl_block_uint *gsl_block_uint_alloc (const size_t n);
GSL_FUN gsl_block_uint *gsl_block_uint_calloc (const size_t n);
GSL_FUN void gsl_block_uint_free (gsl_block_uint * b);

GSL_FUN int gsl_block_uint_fread (FILE * stream, gsl_block_uint * b);
GSL_FUN int gsl_block_uint_fwrite (FILE * stream, const gsl_block_uint * b);
GSL_FUN int gsl_block_uint_fscanf (FILE * stream, gsl_block_uint * b);
GSL_FUN int gsl_block_uint_fprintf (FILE * stream, const gsl_block_uint * b, const char *format);

GSL_FUN int gsl_block_uint_raw_fread (FILE * stream, unsigned int * b, const size_t n, const size_t stride);
GSL_FUN int gsl_block_uint_raw_fwrite (FILE * stream, const unsigned int * b, const size_t n, const size_t stride);
GSL_FUN int gsl_block_uint_raw_fscanf (FILE * stream, unsigned int * b, const size_t n, const size_t stride);
GSL_FUN int gsl_block_uint_raw_fprintf (FILE * stream, const unsigned int * b, const size_t n, const size_t stride, const char *format);

GSL_FUN size_t gsl_block_uint_size (const gsl_block_uint * b);
GSL_FUN unsigned int * gsl_block_uint_data (const gsl_block_uint * b);

__END_DECLS

#endif /* __GSL_BLOCK_UINT_H__ */

ext: h

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_plain.h>
#include <gsl/gsl_monte_miser.h>
#include <gsl/gsl_monte_vegas.h>

double calculateIntegralValue(int sampleCount)
{
    int correctSamples = 0;

    srand(time(NULL));
    FILE* file = fopen("sqrt.txt", "w");

    for (int i = 0; i < sampleCount; ++i){
        double x = (double)rand() / RAND_MAX; // double [0, 1]
        double y = (double)rand() / RAND_MAX; // double [0, 1000]

        if (y <= 1.0/sqrt(x)){
            fprintf (file,"%g %g\n", x, y); //Save correct samples to file
            correctSamples++;
        }
    }
    return (double)correctSamples / (double)sampleCount; // (b-a) is (1.0-0.0) so we can ignore this multiply
}

void calculateIntegralError(int maxSamples)
{
    srand(time(NULL));
    FILE* file = fopen("sqrt_err.txt", "w");

    for (int sampleCount = 1; sampleCount < maxSamples; ++sampleCount){
        int correctSamples = 0;
        for (int i = 0; i < sampleCount; ++i){
            double x = (double)rand() / RAND_MAX; // double [0, 1]
            double y = (double)rand() / RAND_MAX; // double [0, 1]

            if (y <= x * x){
                correctSamples++;
            }
        }

        double sum = (double)correctSamples / (double)sampleCount; // (b-a) is (1.0-0.0) so we can ignore this multiply
        double sumError = sum - 1.0/3.0;
        fprintf(file, "%d %g\n", sampleCount, sqrt(sumError * sumError));
    }
}

void monteCarloPLAIN()
{
    double res, err;
    double xl[3] = {0,0,0};
    double xu[3] = {M_PI, M_PI, M_PI};

    const gsl_rng_type *T;
    gsl_rng *r;

    gsl_monte_function G;

}

int main(void)
{
    // printf("Sum: %g\n", sum);

    return 0;
}

ext: c


text:
#pragma once

#include <gsl/gsl>

#include <mutex>
#include <condition_variable>

#include <arcana/threading/affinity.h>

namespace Babylon
{
    class SafeTimespanGuarantor
    {
    public:
        SafeTimespanGuarantor();

        void BeginSafeTimespan();
        void EndSafeTimespan();

        using SafetyGuarantee = gsl::final_action<std::function<void()>>;
        SafetyGuarantee GetSafetyGuarantee();

    private:
        arcana::affinity m_affinity{};
        uint32_t m_count{};
        std::mutex m_mutex{};
        std::unique_lock<std::mutex> m_lock{};
        std::condition_variable m_condition{};
    };
}


text:
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <stdbool.h>
#include <unistd.h>
#include <dirent.h>
#include "hdf5.h"
#include <math.h>
#include <time.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_sf_bessel.h>
#include "mclib_3d.h"
#include <omp.h>

#define R_DIM 1260
#define THETA_DIM 280
#define PHI_DIM 280

void read_hydro(char hydro_prefix[200], int frame, double r_inj, double **x, double **y,  double **z, double **szx, double **szy, double **r, double **theta, double **phi,\
  double **velx, double **vely, double **velz, double **dens, double **pres, double **gamma, double **dens_lab, double **temp, int *number,  int ph_inj_switch, double min_r, double max_r, double fps, FILE *fPtr)
{
    FILE *hydroPtr=NULL;
    char hydrofile[200]="", file_num[200]="", full_file[200]="",file_end[200]=""  ;
    char buf[10]="";
    int i=0, j=0, k=0, elem=0, elem_factor=0;
    int phi_min_index=0, phi_max_index=0, r_min_index=0, r_max_index=0, theta_min_index=0, theta_max_index=0; //all_index_buffer contains phi_min, phi_max, theta_min, theta_max, r_min, r_max indexes to get from grid files
    int r_index=0, theta_index=0, phi_index=0, hydro_index=0, all_index_buffer=0, adjusted_remapping_index=0, dr_index=0;
    int *remapping_indexes=NULL;
    float buffer=0;
    float *dens_unprc=NULL;
    float *vel_r_unprc=NULL;
    float *vel_theta_unprc=NULL;
    float *vel_phi_unprc=NULL;
    float *pres_unprc=NULL;
    double ph_rmin=0, ph_rmax=0;
    double r_in=1e10, r_ref=2e13;
    double *r_edge=NULL;
    double *dr=NULL;
    double *r_unprc=malloc(sizeof(double)*R_DIM);
    double *theta_unprc=malloc(sizeof(double)*THETA_DIM);
    double *phi_unprc=malloc(sizeof(double)*PHI_DIM);

    if (ph_inj_switch==0)
    {
        ph_rmin=min_r;
        ph_rmax=max_r;
    }

    snprintf(file_end,sizeof(file_end),"%s","small.data" );

    //density
    snprintf(hydrofile,sizeof(hydrofile),"%s%s%d%s",hydro_prefix,"u0", 1,"-" );
    modifyFlashName(file_num, hydrofile, frame,1);

    fprintf(fPtr,">> Opening file %s\n", file_num);
    fflush(fPtr);

    snprintf(full_file, sizeof(full_file), "%s%s", file_num, file_end);
    /*
    fprintf(fPtr,"Reading Density: %s\n", full_file);
    fflush(fPtr);
    */
    hydroPtr=fopen(full_file, "rb");
    fread(&buffer, sizeof(float), 1,hydroPtr); //random stuff about the file from fortran
    fread(&phi_min_index, sizeof(int)*1, 1,hydroPtr); //min and max indexes for the grid
    fread(&phi_max_index, sizeof(int)*1, 1,hydroPtr);
    fread(&theta_min_index, sizeof(int)*1, 1,hydroPtr);
    fread(&theta_max_index, sizeof(int)*1, 1,hydroPtr);
    fread(&r_min_index, sizeof(int)*1, 1,hydroPtr);
    fread(&r_max_index, sizeof(int)*1, 1,hydroPtr);
    fclose(hydroPtr);

    //fortran indexing starts @ 1, but C starts @ 0
    r_min_index--;
    r_max_index--;
    theta_min_index--;
    theta_max_index--;
    phi_min_index--;
    phi_max_index--;

    //number of elements defined by this now
    elem=(r_max_index+1-r_min_index)*(theta_max_index+1-theta_min_index)*(phi_max_index+1-phi_min_index); //add 1 b/c max_index is 1 less than max number of elements in file
    /*
    fprintf(fPtr,"Elem %d\n", elem);
    fprintf(fPtr,"Limits %d, %d, %d, %d, %d, %d\n", phi_min_index, phi_max_index, theta_min_index, theta_max_index, r_min_index, r_max_index);
    fflush(fPtr);
    */
    //now with number of elements allocate data, remember last element is some garbage that only fortran uses
    dens_unprc=malloc(elem*sizeof(float));
    vel_r_unprc=malloc(elem*sizeof(float));
    vel_theta_unprc=malloc(elem*sizeof(float));
    pres_unprc=malloc(elem*sizeof(float));
    vel_phi_unprc=malloc(elem*sizeof(float));

    hydroPtr=fopen(full_file, "rb");
    fread(&buffer, sizeof(float), 1,hydroPtr); //random stuff about the file from fortran
    fread(&all_index_buffer, sizeof(int)*1, 1,hydroPtr); //min and max indexes for the grid, dont need anymore so just save to dummy variable
    fread(&all_index_buffer, sizeof(int)*1, 1,hydroPtr);
    fread(&all_index_buffer, sizeof(int)*1, 1,hydroPtr);
    fread(&all_index_buffer, sizeof(int)*1, 1,hydroPtr);
    fread(&all_index_buffer, sizeof(int)*1, 1,hydroPtr);
    fread(&all_index_buffer, sizeof(int)*1, 1,hydroPtr);
    fread(&buffer, sizeof(float), 1,hydroPtr); //random stuff about the file from fortran
    fread(&buffer, sizeof(float), 1,hydroPtr); //random stuff about the file from fortran

    fread(dens_unprc, sizeof(float),elem, hydroPtr); //data
    fclose(hydroPtr);

    /*
    for (i=0;i<R_DIM*THETA_DIM*PHI_DIM;i++)
    {
        if ((i>98784000-5) || (i<5))
        {
            fprintf(fPtr,"Density %d: %0.7e\n", i, *(dens_unprc+i));
            fflush(fPtr);
        }
    }
    */

    //velocities divided by c
    //v_r
    snprintf(hydrofile,sizeof(hydrofile),"%s%s%d%s",hydro_prefix,"u0", 2,"-" );
    modifyFlashName(file_num, hydrofile, frame,1);
    snprintf(full_file, sizeof(full_file), "%s%s", file_num, file_end);
    /*
    fprintf(fPtr,"Reading v_r: %s\n", full_file);
    fflush(fPtr);
    */
    hydroPtr=fopen(full_file, "rb");
    fread(&buffer, sizeof(float), 1,hydroPtr); //random stuff about the file from fortran
    fread(&all_index_buffer, sizeof(int)*1, 1,hydroPtr); //min and max indexes for the grid, dont need anymore so just save to dummy variable
    fread(&all_index_buffer, sizeof(int)*1, 1,hydroPtr);
    fread(&all_index_buffer, sizeof(int)*1, 1,hydroPtr);
    fread(&all_index_buffer, sizeof(int)*1, 1,hydroPtr);
    fread(&all_index_buffer, sizeof(int)*1, 1,hydroPtr);
    fread(&all_index_buffer, sizeof(int)*1, 1,hydroPtr);
    fread(&buffer, sizeof(float), 1,hydroPtr); //random stuff about the file from fortran
    fread(&buffer, sizeof(float), 1,hydroPtr); //random stuff about the file from fortran

    fread(vel_r_unprc, sizeof(float),elem, hydroPtr);
    fclose(hydroPtr);

    /*
    for (i=0;i<5;i++)
    {
        fprintf(fPtr,"V_r %d: %e\n", i, *(vel_r_unprc+i));
        fflush(fPtr);
    }
    */

    //v_theta
    snprintf(hydrofile,sizeof(hydrofile),"%s%s%d%s",hydro_prefix,"u0", 3,"-" );
    modifyFlashName(file_num, hydrofile, frame,1);
    snprintf(full_file, sizeof(full_file), "%s%s", file_num, file_end);
    /*
    fprintf(fPtr,"Reading v_theta: %s\n", full_file);
    fflush(fPtr);
    */
    hydroPtr=fopen(full_file, "rb");
    fread(&buffer, sizeof(float), 1,hydroPtr); //random stuff about the file from fortran
    fread(&all_index_buffer, sizeof(int)*1, 1,hydroPtr); //min and max indexes for the grid, dont need anymore so just save to dummy variable
    fread(&all_index_buffer, sizeof(int)*1, 1,hydroPtr);
    fread(&all_index_buffer, sizeof(int)*1, 1,hydroPtr);
    fread(&all_index_buffer, sizeof(int)*1, 1,hydroPtr);
    fread(&all_index_buffer, sizeof(int)*1, 1,hydroPtr);
    fread(&all_index_buffer, sizeof(int)*1, 1,hydroPtr);
    fread(&buffer, sizeof(float), 1,hydroPtr); //random stuff about the file from fortran
    fread(&buffer, sizeof(float), 1,hydroPtr); //random stuff about the file from fortran

    fread(vel_theta_unprc, sizeof(float),elem, hydroPtr);
    fclose(hydroPtr);

    /*
    for (i=0;i<5;i++)
    {
        fprintf(fPtr,"V_theta %d: %e\n", i, *(vel_theta_unprc+i));
        fflush(fPtr);
    }
    */

    //v_phi
    snprintf(hydrofile,sizeof(hydrofile),"%s%s%d%s",hydro_prefix,"u0", 4,"-" );
    modifyFlashName(file_num, hydrofile, frame,1);
    snprintf(full_file, sizeof(full_file), "%s%s", file_num, file_end);
    /*
    fprintf(fPtr,"Reading v_phi: %s\n", full_file);
    fflush(fPtr);
    */
    hydroPtr=fopen(full_file, "rb");
    fread(&buffer, sizeof(float), 1,hydroPtr); //random stuff about the file from fortran
    fread(&all_index_buffer, sizeof(int)*1, 1,hydroPtr); //min and max indexes for the grid, dont need anymore so just save to dummy variable
    fread(&all_index_buffer, sizeof(int)*1, 1,hydroPtr);
    fread(&all_index_buffer, sizeof(int)*1, 1,hydroPtr);
    fread(&all_index_buffer, sizeof(int)*1, 1,hydroPtr);
    fread(&all_index_buffer, sizeof(int)*1, 1,hydroPtr);
    fread(&all_index_buffer, sizeof(int)*1, 1,hydroPtr);
    fread(&buffer, sizeof(float), 1,hydroPtr); //random stuff about the file from fortran
    fread(&buffer, sizeof(float), 1,hydroPtr); //random stuff about the file from fortran

    fread(vel_phi_unprc, sizeof(float),elem, hydroPtr);
    fclose(hydroPtr);

    /*
    for (i=0;i<5;i++)
    {
        fprintf(fPtr,"V_phi %d: %e\n", i, *(vel_phi_unprc+i));
        fflush(fPtr);
    }
    */

    //pressure (divided by c^2)
    snprintf(hydrofile,sizeof(hydrofile),"%s%s%d%s",hydro_prefix,"u0", 8,"-" );
    modifyFlashName(file_num, hydrofile, frame,1);
    snprintf(full_file, sizeof(full_file), "%s%s", file_num, file_end);
    /*
    fprintf(fPtr,"Reading pres: %s\n", full_file);
    fflush(fPtr);
    */
    hydroPtr=fopen(full_file, "rb");
    fread(&buffer, sizeof(float), 1,hydroPtr); //random stuff about the file from fortran
    fread(&all_index_buffer, sizeof(int)*1, 1,hydroPtr); //min and max indexes for the grid, dont need anymore so just save to dummy variable
    fread(&all_index_buffer, sizeof(int)*1, 1,hydroPtr);
    fread(&all_index_buffer, sizeof(int)*1, 1,hydroPtr);
    fread(&all_index_buffer, sizeof(int)*1, 1,hydroPtr);
    fread(&all_index_buffer, sizeof(int)*1, 1,hydroPtr);
    fread(&all_index_buffer, sizeof(int)*1, 1,hydroPtr);
    fread(&buffer, sizeof(float), 1,hydroPtr); //random stuff about the file from fortran
    fread(&buffer, sizeof(float), 1,hydroPtr); //random stuff about the file from fortran

    fread(pres_unprc, sizeof(float),elem, hydroPtr);
    fclose(hydroPtr);
    /*
    for (i=PHI_DIM-1;i<PHI_DIM;i++)
    {
        for (j=THETA_DIM-1;j<THETA_DIM;j++)
        {
            for (k=R_DIM-5;k<R_DIM;k++)
            {

            fprintf(fPtr,"Pres %d: %e\n", (i*R_DIM*THETA_DIM + j*R_DIM + k  ), *(pres_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k  )));
            fflush(fPtr);

            }
        }
    }
    */

    // see how many elements there are to test if reading correctly
    /*
    hydroPtr=fopen(full_file, "rb");
    fread(&buffer, sizeof(float), 1,hydroPtr); //random stuff about the file from fortran
    while (1 == fread(&buffer,sizeof(float),1,hydroPtr))
    {
        elem++;
    }
    //fread(pres_unprc, sizeof(double)*HYDRO_DIM,HYDRO_DIM, hydroPtr);
    fclose(hydroPtr);

    fprintf(fPtr,"Elem %d\n", elem);
    */

    //R
    //remapping_indexes=getIndexesForRadialRemapping(hydro_prefix); //can run this once in debug mode to find out delta index for each remapping and number of total r elements
    //for given set of remappings on July 12th 2017, grid00-x1.data[420]=grid01-x1.data[0], grid01-x1.data[420]=grid02-x1.data[0], etc. total number of r is 3780

    if (frame<=1300)
    {
        snprintf(hydrofile,sizeof(hydrofile),"%s%s%d%s",hydro_prefix,"grid0", 0,"-x1.data" );
        //adjusted_remapping_index=(*(remapping_indexes+0))+r_min_index ; //this if I am not hardcoding the dr index values
        adjusted_remapping_index=(0*420)+r_min_index;
    }
    else if (frame<=2000)
    {
        snprintf(hydrofile,sizeof(hydrofile),"%s%s%d%s",hydro_prefix,"grid0", 1,"-x1.data" );
        //adjusted_remapping_index=(*(remapping_indexes+1))+r_min_index;
        adjusted_remapping_index=(1*420)+r_min_index;
    }
    else if (frame<=10000)
    {
        snprintf(hydrofile,sizeof(hydrofile),"%s%s%d%s",hydro_prefix,"grid0", 2,"-x1.data" );
        //adjusted_remapping_index=(*(remapping_indexes+2))+r_min_index;
        adjusted_remapping_index=(2*420)+r_min_index;
    }
    else if (frame<=20000)
    {
        snprintf(hydrofile,sizeof(hydrofile),"%s%s%d%s",hydro_prefix,"grid0", 3,"-x1.data" );
        //adjusted_remapping_index=(*(remapping_indexes+3))+r_min_index;
        adjusted_remapping_index=(3*420)+r_min_index;
    }
    else if (frame<=35000)
    {
        snprintf(hydrofile,sizeof(hydrofile),"%s%s%d%s",hydro_prefix,"grid0", 4,"-x1.data" );
        //adjusted_remapping_index=(*(remapping_indexes+4))+r_min_index;
        adjusted_remapping_index=(4*420)+r_min_index;
    }
    else if (frame<=50000)
    {
        snprintf(hydrofile,sizeof(hydrofile),"%s%s%d%s",hydro_prefix,"grid0", 5,"-x1.data" );
        //adjusted_remapping_index=(*(remapping_indexes+5))+r_min_index;
        adjusted_remapping_index=(5*420)+r_min_index;
    }
    else if (frame<=60000)
    {
        snprintf(hydrofile,sizeof(hydrofile),"%s%s%d%s",hydro_prefix,"grid0", 6,"-x1.data" );
        //adjusted_remapping_index=(*(remapping_indexes+6))+r_min_index;
        adjusted_remapping_index=(6*420)+r_min_index;
    }

    fprintf(fPtr,"Reading Radius: %s\n", hydrofile);
    fflush(fPtr);

    hydroPtr=fopen(hydrofile, "r");

    i=0;
    while (i<R_DIM)
    {
        fscanf(hydroPtr, "%lf", (r_unprc+i));  //read value
        fgets(buf, 3,hydroPtr); //read comma
        /*
        if (i<5)
        {
            fprintf(fPtr,"R %d: %e\n", i, *(r_unprc+i));
            fflush(fPtr);
        }
        */
        i++;
    }

    fclose(hydroPtr);

    r_edge=malloc(sizeof(double)*(3780+1));
    dr=malloc(sizeof(double)*(3780));

    //calculate radial grid edges
    *(r_edge+0)=r_in;
    i=0;
    for (i=1;i<3780;i++)
    {
        *(r_edge+i)=(*(r_edge+i-1))+((*(r_edge+i-1))*(M_PI/560)/(1+((*(r_edge+i-1))/r_ref))); //r_i = r_(i-1) + Dq r_(i-1) [1 + r_(i-1)/r0]-1
        *(dr+i-1)=(*(r_edge+i))-(*(r_edge+i-1));
        /*
        if (i<5)
        {
            fprintf(fPtr,"R Edge: %d: %e Dr: %e\n", i, *(r_edge+i), *(dr+i-1));
            fflush(fPtr);
        }
        */
    }
    free(r_edge);

    //Theta
    snprintf(hydrofile,sizeof(hydrofile),"%s%s",hydro_prefix,"grid-x2.data" );

    fprintf(fPtr,"Reading Theta: %s\n", hydrofile);
    fflush(fPtr);

    hydroPtr=fopen(hydrofile, "r");

    i=0;
    while (i<THETA_DIM)
    {
        fscanf(hydroPtr, "%lf", (theta_unprc+i));  //read value
        fgets(buf, 3,hydroPtr); //read comma
        /*
        if (i<5)
        {
            fprintf(fPtr,"R %d: %e\n", i, *(theta_unprc+i));
            fflush(fPtr);
        }
        */
        i++;
    }

    fclose(hydroPtr);

    //Phi
    snprintf(hydrofile,sizeof(hydrofile),"%s%s",hydro_prefix,"grid-x3.data" );
    /*
    fprintf(fPtr,"Reading Phi: %s\n", hydrofile);
    fflush(fPtr);
    */
    hydroPtr=fopen(hydrofile, "r");

    i=0;
    while (i<PHI_DIM)
    {
        fscanf(hydroPtr, "%lf", (phi_unprc+i));  //read value
        fgets(buf, 3,hydroPtr); //read comma
        /*
        if (i<5)
        {
            fprintf(fPtr,"R %d: %e\n", i, *(phi_unprc+i));
            fflush(fPtr);
        }
        */
        i++;
    }

    fclose(hydroPtr);

    //limit number of array elements PUT WHILE LOOP TO MAKE SURE NUMBER OF ELEMENTS >0
    elem_factor=0;
    elem=0;
    while (elem==0)
    {
        elem=0;
        elem_factor++;
        for (i=0;i<(phi_max_index+1-phi_min_index);i++)
        {
            for (j=0;j<(theta_max_index+1-theta_min_index);j++)
            {
                for (k=0;k<(r_max_index+1-r_min_index);k++)
                {
                    r_index=r_min_index+k;
                    //if I have photons do selection differently than if injecting photons
                    if (ph_inj_switch==0)
                    {
                        //printf("R's:%d, %e\n", k, *(r_unprc+r_index));
                        //if calling this function when propagating photons, choose blocks based on where the photons are
                        if (((ph_rmin - elem_factor*C_LIGHT/fps)<(*(r_unprc+r_index))) && (*(r_unprc+r_index)  < (ph_rmax + elem_factor*C_LIGHT/fps) ))
                        {
                            // *(pres_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k  )
                            elem++;
                        }
                    }
                    else
                    {
                        //if calling this function to inject photons choose blocks based on injection parameters, r_inj, which is sufficient
                        if (((r_inj - elem_factor*C_LIGHT/fps)<(*(r_unprc+r_index))) && (*(r_unprc+r_index)  < (r_inj + elem_factor*C_LIGHT/fps) ))
                        {
                            // *(pres_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k  )
                            elem++;
                        }

                    }

                }
            }
        }
    }

    fprintf(fPtr,"Number of post restricted Elems: %d %e\n", elem, r_inj);
    //fprintf(fPtr,"Ph_min, Ph_max: %e, %e\n With c: min: %e max: %e \n", ph_rmin, ph_rmax, (ph_rmin - 2*C_LIGHT/fps), (ph_rmax + 2*C_LIGHT/fps));
    fflush(fPtr);

    //allocate space for new set of data
    (*pres)=malloc (elem * sizeof (double ));
    (*velx)=malloc (elem * sizeof (double ));
    (*vely)=malloc (elem * sizeof (double ));
    (*velz)=malloc (elem * sizeof (double ));
    (*dens)=malloc (elem * sizeof (double ));
    (*x)=malloc (elem * sizeof (double ));
    (*y)=malloc (elem * sizeof (double ));
    (*z)=malloc (elem * sizeof (double ));
    (*r)=malloc (elem * sizeof (double ));
    (*theta)=malloc (elem * sizeof (double ));
    (*phi)=malloc (elem * sizeof (double ));
    (*gamma)=malloc (elem * sizeof (double ));
    (*dens_lab)=malloc (elem * sizeof (double ));
    (*szx)=malloc (elem * sizeof (double )); //theta and phi resolution
    (*szy)=malloc (elem * sizeof (double )); //r resolution
    (*temp)=malloc (elem * sizeof (double ));

    //limit number of array elements
    elem=0;
    for (i=0;i<(phi_max_index+1-phi_min_index);i++)
    {
        for (j=0;j<(theta_max_index+1-theta_min_index);j++)
        {
            for (k=0;k<(r_max_index+1-r_min_index);k++)
            {
                r_index=r_min_index+k; //look at indexes of r that are included in small hydro file
                theta_index=theta_min_index+j;
                phi_index=phi_min_index+i;
                dr_index=adjusted_remapping_index+k;
                hydro_index=(i*(r_max_index+1-r_min_index)*(theta_max_index+1-theta_min_index) + j*(r_max_index+1-r_min_index) + k  );

                //if I have photons do selection differently than if injecting photons
                if (ph_inj_switch==0)
                {
                    //if calling this function when propagating photons, choose blocks based on where the photons are
                    if (((ph_rmin - elem_factor*C_LIGHT/fps)<(*(r_unprc+r_index))) && (*(r_unprc+r_index)  < (ph_rmax + elem_factor*C_LIGHT/fps) ))
                    {
                        (*pres)[elem] = *(pres_unprc+hydro_index);
                        (*dens)[elem] = *(dens_unprc+hydro_index);
                        (*temp)[elem] =  pow(3*(*(pres_unprc+hydro_index))*pow(C_LIGHT,2.0)/(A_RAD) ,1.0/4.0);

                        (*gamma)[elem] = pow(pow(1.0-(pow(*(vel_r_unprc+hydro_index),2)+ pow(*(vel_theta_unprc+hydro_index),2)+pow(*(vel_phi_unprc+hydro_index),2)),0.5),-1);
                        (*dens_lab)[elem] = (*(dens_unprc+hydro_index))*pow(pow(1.0-(pow(*(vel_r_unprc+hydro_index),2)+ pow(*(vel_theta_unprc+hydro_index),2)+pow(*(vel_phi_unprc+hydro_index),2)),0.5),-1);
                        (*r)[elem] = *(r_unprc+r_index);
                        (*theta)[elem] = *(theta_unprc+theta_index);
                        (*phi)[elem] = *(phi_unprc+phi_index);
                        (*x)[elem] = (*(r_unprc+r_index))*sin(*(theta_unprc+theta_index))*cos(*(phi_unprc+phi_index));
                        (*y)[elem] = (*(r_unprc+r_index))*sin(*(theta_unprc+theta_index))*sin(*(phi_unprc+phi_index));
                        (*z)[elem] = (*(r_unprc+r_index))*cos(*(theta_unprc+theta_index));
                        (*szx)[elem] =  *(dr+dr_index);
                        (*szy)[elem] =  M_PI/560;
                        (*velx)[elem]=((*(vel_r_unprc+hydro_index))*sin(*(theta_unprc+theta_index))*cos(*(phi_unprc+phi_index))) + ((*(vel_theta_unprc+hydro_index))*cos(*(theta_unprc+theta_index))*cos(*(phi_unprc+phi_index))) - ((*(vel_phi_unprc+hydro_index))*sin(*(phi_unprc+phi_index)));
                        (*vely)[elem]=((*(vel_r_unprc+hydro_index))*sin(*(theta_unprc+theta_index))*sin(*(phi_unprc+phi_index))) + ((*(vel_theta_unprc+hydro_index))*cos(*(theta_unprc+theta_index))*sin(*(phi_unprc+phi_index))) + ((*(vel_phi_unprc+hydro_index))*cos(*(phi_unprc+phi_index)));
                        (*velz)[elem]=((*(vel_r_unprc+hydro_index))*cos(*(theta_unprc+theta_index))) - ((*(vel_theta_unprc+hydro_index))*sin(*(theta_unprc+theta_index)));

                        elem++;

                    }
                }
                else
                {
                    //if calling this function to inject photons choose blocks based on injection parameters, r_inj, which is sufficient
                    if (((r_inj - elem_factor*C_LIGHT/fps)<(*(r_unprc+r_index))) && (*(r_unprc+r_index)  < (r_inj + elem_factor*C_LIGHT/fps) ))
                    {
                        (*pres)[elem] = *(pres_unprc+hydro_index);
                        (*dens)[elem] = *(dens_unprc+hydro_index);
                        (*temp)[elem] =  pow(3*(*(pres_unprc+hydro_index))*pow(C_LIGHT,2.0)/(A_RAD) ,1.0/4.0);
                        (*gamma)[elem] = pow(pow(1.0-(pow(*(vel_r_unprc+hydro_index),2)+ pow(*(vel_theta_unprc+hydro_index),2)+pow(*(vel_phi_unprc+hydro_index),2)),0.5),-1);
                        (*dens_lab)[elem] = (*(dens_unprc+hydro_index))*pow(pow(1.0-(pow(*(vel_r_unprc+hydro_index),2)+ pow(*(vel_theta_unprc+hydro_index),2)+pow(*(vel_phi_unprc+hydro_index),2)),0.5),-1);
                        (*r)[elem] = *(r_unprc+r_index);
                        (*theta)[elem] = *(theta_unprc+theta_index);
                        (*phi)[elem] = *(phi_unprc+phi_index);
                        (*x)[elem] = (*(r_unprc+r_index))*sin(*(theta_unprc+theta_index))*cos(*(phi_unprc+phi_index));
                        (*y)[elem] = (*(r_unprc+r_index))*sin(*(theta_unprc+theta_index))*sin(*(phi_unprc+phi_index));
                        (*z)[elem] = (*(r_unprc+r_index))*cos(*(theta_unprc+theta_index));
                        (*szx)[elem] =  *(dr+dr_index);
                        (*szy)[elem] =  M_PI/560;
                        (*velx)[elem]=((*(vel_r_unprc+hydro_index))*sin(*(theta_unprc+theta_index))*cos(*(phi_unprc+phi_index))) + ((*(vel_theta_unprc+hydro_index))*cos(*(theta_unprc+theta_index))*cos(*(phi_unprc+phi_index))) - ((*(vel_phi_unprc+hydro_index))*sin(*(phi_unprc+phi_index)));
                        (*vely)[elem]=((*(vel_r_unprc+hydro_index))*sin(*(theta_unprc+theta_index))*sin(*(phi_unprc+phi_index))) + ((*(vel_theta_unprc+hydro_index))*cos(*(theta_unprc+theta_index))*sin(*(phi_unprc+phi_index))) + ((*(vel_phi_unprc+hydro_index))*cos(*(phi_unprc+phi_index)));
                        (*velz)[elem]=((*(vel_r_unprc+hydro_index))*cos(*(theta_unprc+theta_index))) - ((*(vel_theta_unprc+hydro_index))*sin(*(theta_unprc+theta_index)));

                        elem++;
                    }

                }

            }
        }
    }

    *number=elem;

    free(pres_unprc); free(dens_unprc); free(r_unprc); free(theta_unprc); free(phi_unprc);free(dr);free(vel_r_unprc); free(vel_theta_unprc); free(vel_phi_unprc);

}

void photonInjection3D( struct photon **ph, int *ph_num, double r_inj, double ph_weight, int min_photons, int max_photons, char spect, int array_length, double fps, double theta_min, double theta_max,\
 double *x, double *y, double *z, double *szx, double *szy, double *r, double *theta, double *phi, double *temps, double *vx, double *vy, double *vz, gsl_rng * rand, FILE *fPtr)
{
    int i=0, block_cnt=0, *ph_dens=NULL, ph_tot=0, j=0,k=0;
    double ph_dens_calc=0.0, fr_dum=0.0, y_dum=0.0, yfr_dum=0.0, fr_max=0, bb_norm=0, position_phi, ph_weight_adjusted, theta_prime=0;
    double com_v_phi, com_v_theta, *p_comv=NULL, *boost=NULL; //comoving phi, theta, comoving 4 momentum for a photon, and boost for photon(to go to lab frame)
    double *l_boost=NULL; //pointer to hold array of lorentz boost, to lab frame, values
    float num_dens_coeff;

    if (spect=='w') //from MCRAT paper, w for wien spectrum
    {
        num_dens_coeff=8.44;
        //printf("in wien spectrum\n");
    }
    else
    {
        num_dens_coeff=20.29; //this is for black body spectrum
        //printf("in BB spectrum");
    }

    //find how many blocks are near the injection radius within the angles defined in mc.par, get temperatures and calculate number of photons to allocate memory for
    //and then rcord which blocks have to have "x" amount of photons injected there
    printf("%e, %e\n",*(phi+i), theta_max);
    for(i=0;i<array_length;i++)
    {
        //look at all boxes in width delta r=c/fps and within angles we are interested in NEED TO modify for RIKEN data- dont need r anymore, just theta and phi? (didnt work), just look at pojection on x-z plane
        theta_prime=acos(*(y+i)/(*(r+i))); //jet axis here is the y axis
            if ( (theta_prime< theta_max) && (theta_prime >= theta_min) ) //(*(r+i) > (r_inj - C_LIGHT/fps))  &&   (*(r+i)  < (r_inj + C_LIGHT/fps)  ) &&
            {
                //printf("%e\n", theta_prime );
                block_cnt++;
            }
    }
    printf("Blocks: %d\n", block_cnt);

    ph_dens=malloc(block_cnt * sizeof(int));

    //calculate the photon density for each block and save it to the array
    j=0;
    ph_tot=0;
    ph_weight_adjusted=ph_weight;
    //printf("%d %d\n", max_photons, min_photons);
    while ((ph_tot>max_photons) || (ph_tot<min_photons) )
    {
        j=0;
        ph_tot=0;
        //allocate memory to record density of photons for each block
        //ph_dens=malloc(block_cnt * sizeof(int));

        for (i=0;i<array_length;i++)
        {
            //printf("%d\n",i);
            //printf("%e, %e, %e, %e, %e, %e\n", *(r+i),(r_inj - C_LIGHT/fps), (r_inj + C_LIGHT/fps), *(theta+i) , theta_max, theta_min);
                //NEED TO modify for RIKEN data - modified
                theta_prime=acos(*(y+i)/(*(r+i)));
                if ( (theta_prime< theta_max) && (theta_prime >= theta_min) )
                {
                    //NEED TO modify for RIKEN data - modified
                    ph_dens_calc=(num_dens_coeff*pow(*(temps+i),3.0)*pow(*(r+i),2)*sin(*(theta+i))* pow(*(szy+i),2.0)*(*(szx+i)) /(ph_weight_adjusted))*pow(pow(1.0-(pow(*(vx+i),2)+pow(*(vy+i),2)+pow(*(vz+i),2)),0.5),-1) ; //a*T^3/(weight) dV, dV=2*PI*x*dx^2,

                     (*(ph_dens+j))=gsl_ran_poisson(rand,ph_dens_calc) ; //choose from poission distribution with mean of ph_dens_calc

                    //printf("%d, %lf \n",*(ph_dens+j), ph_dens_calc);

                     //sum up all the densities to get total number of photons
                     ph_tot+=(*(ph_dens+j));

                     j++;
                }
        }

        if (ph_tot>max_photons)
        {
            //if the number of photons is too big make ph_weight larger
            ph_weight_adjusted*=10;
            //free(ph_dens);
        }
        else if (ph_tot<min_photons)
        {
            ph_weight_adjusted*=0.5;
            //free(ph_dens);
        }
        //printf("dens: %d, photons: %d\n", *(ph_dens+(j-1)), ph_tot);

    }

    printf("%d\n", ph_tot);

    //allocate memory for that many photons and also allocate memory to hold comoving 4 momentum of each photon and the velocity of the fluid
    (*ph)=malloc (ph_tot * sizeof (struct photon ));

    p_comv=malloc(4*sizeof(double));
    boost=malloc(3*sizeof(double));
    l_boost=malloc(4*sizeof(double));

    //go through blocks and assign random energies/locations to proper number of photons
    ph_tot=0;
    k=0;
    for (i=0;i<array_length;i++)
    {
        theta_prime=acos(*(y+i)/(*(r+i)));
        if ( (theta_prime< theta_max) && (theta_prime >= theta_min) )  //NEED TO modify for RIKEN data - modified
        {

            //*(temps+i)=0.76*(*(temps+i));
            for(j=0;j<( *(ph_dens+k) ); j++ )
            {
                    //have to get random frequency for the photon comoving frequency
                    y_dum=1; //initalize loop
                    yfr_dum=0;
                    while (y_dum>yfr_dum)
                    {
                        fr_dum=gsl_rng_uniform_pos(rand)*6.3e11*(*(temps+i)); //in Hz
                        //printf("%lf, %lf ",gsl_rng_uniform_pos(rand), (*(temps+i)));
                        y_dum=gsl_rng_uniform_pos(rand);
                        //printf("%lf ",fr_dum);

                        if (spect=='w')
                        {
                            yfr_dum=(1.0/(1.29e31))*pow((fr_dum/(*(temps+i))),3.0)/(exp((PL_CONST*fr_dum)/(K_B*(*(temps+i)) ))-1); //curve is normalized to maximum
                        }
                        else
                        {
                            fr_max=(5.88e10)*(*(temps+i));//(C_LIGHT*(*(temps+i)))/(0.29); //max frequency of bb
                            bb_norm=(PL_CONST*fr_max * pow((fr_max/C_LIGHT),2.0))/(exp(PL_CONST*fr_max/(K_B*(*(temps+i))))-1); //find value of bb at fr_max
                            yfr_dum=((1.0/bb_norm)*PL_CONST*fr_dum * pow((fr_dum/C_LIGHT),2.0))/(exp(PL_CONST*fr_dum/(K_B*(*(temps+i))))-1); //curve is normalized to vaue of bb @ max frequency
                        }
                        //printf("%lf, %lf,%lf,%e \n",(*(temps+i)),fr_dum, y_dum, yfr_dum);

                    }
                   //printf("%lf\n ",fr_dum);
                   //position_phi= gsl_rng_uniform(rand)*2*M_PI; //NEED TO modify for RIKEN data-modified, dont need anymore
                   com_v_phi=gsl_rng_uniform(rand)*2*M_PI;
                   com_v_theta=acos((gsl_rng_uniform(rand)*2)-1);
                   //printf("%lf, %lf, %lf\n", position_phi, com_v_phi, com_v_theta);

                   //populate 4 momentum comoving array
                   *(p_comv+0)=PL_CONST*fr_dum/C_LIGHT;
                   *(p_comv+1)=(PL_CONST*fr_dum/C_LIGHT)*sin(com_v_theta)*cos(com_v_phi);
                   *(p_comv+2)=(PL_CONST*fr_dum/C_LIGHT)*sin(com_v_theta)*sin(com_v_phi);
                   *(p_comv+3)=(PL_CONST*fr_dum/C_LIGHT)*cos(com_v_theta);

                    //populate boost matrix, not sure why multiplying by -1, seems to give correct answer in old python code...
                    //NEED TO modify for RIKEN data - modified
                    *(boost+0)=-1*(*(vx+i));
                    *(boost+1)=-1*(*(vy+i));
                    *(boost+2)=-1*(*(vz+i));
                    //printf("%lf, %lf, %lf\n", *(boost+0), *(boost+1), *(boost+2));

                    //boost to lab frame
                    lorentzBoost(boost, p_comv, l_boost, 'p', fPtr);
                    //printf("Assignemnt: %e, %e, %e, %e\n", *(l_boost+0), *(l_boost+1), *(l_boost+2),*(l_boost+3));

                (*ph)[ph_tot].p0=(*(l_boost+0));
                (*ph)[ph_tot].p1=(*(l_boost+1));
                (*ph)[ph_tot].p2=(*(l_boost+2));
                (*ph)[ph_tot].p3=(*(l_boost+3));
                //NEED TO modify for RIKEN data-modified
                (*ph)[ph_tot].r0= (*(x+i)); //put photons @ center of box that they are supposed to be in with random phi
                (*ph)[ph_tot].r1=(*(y+i)) ;
                (*ph)[ph_tot].r2=(*(z+i)); //y coordinate in flash becomes z coordinate in MCRaT
                (*ph)[ph_tot].num_scatt=0;
                (*ph)[ph_tot].weight=ph_weight_adjusted;
                //printf("%d\n",ph_tot);
                ph_tot++;
            }
            k++;
        }
    }

    *ph_num=ph_tot; //save number of photons
    //printf(" %d: %d\n", *(ph_dens+(k-1)), *ph_num);
    free(ph_dens); free(p_comv);free(boost); free(l_boost);

}

void phMinMax(struct photon *ph, int ph_num, double *min, double *max, double *min_theta, double *max_theta, FILE *fPtr)
{
    double temp_r_max=0, temp_r_min=DBL_MAX, temp_theta_max=0, temp_theta_min=DBL_MAX;
    int i=0, num_thread=omp_get_num_threads();
    double ph_r=0, ph_theta=0;

#pragma omp parallel for num_threads(num_thread) firstprivate(ph_r, ph_theta) reduction(min:temp_r_min) reduction(max:temp_r_max) reduction(min:temp_theta_min) reduction(max:temp_theta_max)
    for (i=0;i<ph_num;i++)
    {
        if ((ph+i)->weight != 0)
        {
            ph_r=pow(pow( ((ph+i)->r0), 2.0) + pow(((ph+i)->r1),2.0 ) + pow(((ph+i)->r2) , 2.0),0.5);
            ph_theta=acos(((ph+i)->r2) /ph_r); //this is the photons theta psition in the FLASH grid, gives in radians
            if (ph_r > temp_r_max )
            {
                temp_r_max=ph_r;
                //fprintf(fPtr, "The new max is: %e from photon %d with x: %e y: %e z: %e\n", temp_r_max, i, ((ph+i)->r0), (ph+i)->r1, (ph+i)->r2);
            }

            //if ((i==0) || (ph_r<temp_r_min))
            if (ph_r<temp_r_min)
            {
                temp_r_min=ph_r;
                //fprintf(fPtr, "The new min is: %e from photon %d with x: %e y: %e z: %e\n", temp_r_min, i, ((ph+i)->r0), (ph+i)->r1, (ph+i)->r2);
            }

            if (ph_theta > temp_theta_max )
            {
                temp_theta_max=ph_theta;
                //fprintf(fPtr, "The new max is: %e from photon %d with x: %e y: %e z: %e\n", temp_r_max, i, ((ph+i)->r0), (ph+i)->r1, (ph+i)->r2);
            }

            //if ((i==0) || (ph_r<temp_r_min))
            if (ph_theta<temp_theta_min)
            {
                temp_theta_min=ph_theta;
                //fprintf(fPtr, "The new min is: %e from photon %d with x: %e y: %e z: %e\n", temp_r_min, i, ((ph+i)->r0), (ph+i)->r1, (ph+i)->r2);
            }
        }
    }

    *max=temp_r_max;
    *min=temp_r_min;
    *max_theta=temp_theta_max;
    *min_theta=temp_theta_min;
}

int *getIndexesForRadialRemapping(char hydro_prefix[200])
{
    FILE *hydroPtr=NULL;
    char hydrofile[200]="";
    char buf[10]="";
    int i=0, j=0;
    int *remapping_start_index=malloc(sizeof(int)*7); //what index out of total range of r does each remapping begin at
    double r_in=1e10, r_ref=2e13;
    double *r_unprc_0=malloc(sizeof(double)*R_DIM), *r_unprc_1=malloc(sizeof(double)*R_DIM), *r_unprc_2=malloc(sizeof(double)*R_DIM), *r_unprc_3=malloc(sizeof(double)*R_DIM);
    double *r_unprc_4=malloc(sizeof(double)*R_DIM), *r_unprc_5=malloc(sizeof(double)*R_DIM), *r_unprc_6=malloc(sizeof(double)*R_DIM);
    double *r_edge=NULL, *dr=NULL, *rPtr=NULL;

    for (i=0;i<7;i++)
    {
        snprintf(hydrofile,sizeof(hydrofile),"%s%s%d%s",hydro_prefix,"grid0", i,"-x1.data" );

        hydroPtr=fopen(hydrofile, "r");

        j=0;
        while (j<R_DIM)
        {
            switch (i)
            {
                case 0: fscanf(hydroPtr, "%lf", (r_unprc_0+j));  //read value
                case 1: fscanf(hydroPtr, "%lf", (r_unprc_1+j));
                case 2: fscanf(hydroPtr, "%lf", (r_unprc_2+j));
                case 3: fscanf(hydroPtr, "%lf", (r_unprc_3+j));
                case 4: fscanf(hydroPtr, "%lf", (r_unprc_4+j));
                case 5: fscanf(hydroPtr, "%lf", (r_unprc_5+j));
                case 6: fscanf(hydroPtr, "%lf", (r_unprc_6+j));

            }
            fgets(buf, 3,hydroPtr); //read comma
            /*
            if (i<5)
            {
                fprintf(fPtr,"R %d: %e\n", i, *(r_unprc+i));
                fflush(fPtr);
            }
            */
            j++;
        }

        fclose(hydroPtr);
    }

    //calculate the indexes in which each remapping takes over
    j=0; //keeps track of indexes of all of the possible r values
    i=0; //keeps track of index within a certain remapping, when get to R_DIM know that were @ end of last remapping b/c remappings overlap with one another
    rPtr=r_unprc_0; //start off looking at 0th remapping
    *(remapping_start_index+1)=j; //0th remapping starts at index 0
    while (i<R_DIM)
    {
        if (*(rPtr+i)== *(r_unprc_1+0))
        {
            //if the element of the 0th remapping is equal to the 1st element of the 1st remapping, start to look at the 1st remapping
            rPtr=r_unprc_1;
            i=0;
            *(remapping_start_index+1)=j;
        }
        else if (*(rPtr+i)== *(r_unprc_2+0))
        {
            rPtr=r_unprc_2;
            i=0;
            *(remapping_start_index+2)=j;
        }
        else if (*(rPtr+i)== *(r_unprc_3+0))
        {
            rPtr=r_unprc_3;
            i=0;
            *(remapping_start_index+3)=j;
        }
        else if (*(rPtr+i)== *(r_unprc_4+0))
        {
            rPtr=r_unprc_4;
            i=0;
            *(remapping_start_index+4)=j;
        }
        else if (*(rPtr+i)== *(r_unprc_5+0))
        {
            rPtr=r_unprc_5;
            i=0;
            *(remapping_start_index+5)=j;
        }
        else if (*(rPtr+i)== *(r_unprc_6+0))
        {
            rPtr=r_unprc_6;
            i=0;
            *(remapping_start_index+6)=j;
        }

        j++;
        i++;
    }

    printf("Indexes %d, %d, %d, %d, %d, %d, %d\n Elems: %d\n", *(remapping_start_index+0), *(remapping_start_index+1), *(remapping_start_index+2), *(remapping_start_index+3), *(remapping_start_index+4), *(remapping_start_index+5), *(remapping_start_index+6), j);
    //exit(0);
    r_edge=malloc(sizeof(double)*(j+1));
    dr=malloc(sizeof(double)*j);

    //calculate radial grid edges
    *(r_edge+0)=r_in;
    i=0;
    for (i=1;i<j;i++)
    {
        *(r_edge+i)=(*(r_edge+i-1))+((*(r_edge+i-1))*(M_PI/560)/(1+((*(r_edge+i-1))/r_ref))); //r_i = r_(i-1) + Dq r_(i-1) [1 + r_(i-1)/r0]-1
        *(dr+i-1)=(*(r_edge+i))-(*(r_edge+i-1));
        /*
        if (i<5)
        {
            fprintf(fPtr,"R Edge: %d: %e Dr: %e\n", i, *(r_edge+i), *(dr+i-1));
            fflush(fPtr);
        }
        */
    }
    free(r_edge);
    free(r_unprc_0);
    free(r_unprc_1);
    free(r_unprc_2);
    free(r_unprc_3);
    free(r_unprc_4);
    free(r_unprc_5);
    free(r_unprc_6);

    return remapping_start_index;

}

max_stars_repo_name: outflows/MCRaT

max_issues_repo_name: outflows/MCRaT

max_forks_repo_name: outflows/MCRaT


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

#include <cblas.h>
#include <clBlas.h>

void
cblas_dgemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA,
	    const enum CBLAS_TRANSPOSE TransB, const int M, const int N,
	    const int K, const double alpha, const double *A,
	    const int lda, const double *B, const int ldb,
	    const double beta, double *C, const int ldc)
{
	cl_int		err;
	cl_platform_id	platform = 0;
	cl_device_id	device = 0;
	cl_context_properties props[3] = {CL_CONTEXT_PLATFORM, 0, 0};
	cl_context	ctx = 0;
	cl_command_queue queue = 0;
	cl_mem		bufA   , bufB, bufC;
	cl_event	event = NULL;
	cl_uint available = 0;
	int		ret = 0;

	// HACK: ignore order / trans
    clblasOrder order = clblasColumnMajor;
    clblasTranspose transA = clblasNoTrans;
    clblasTranspose transB = clblasNoTrans;

	/* Setup OpenCL environment. */
	err = clGetPlatformIDs(1, &platform, &available);
	if (err != CL_SUCCESS) {
		printf("clGetPlatformIDs() failed with %d and %d\n", err, available);
		exit(1);
	}
	printf("found %d OpenCL platforms\n", available);

	// CL_DEVICE_TYPE_GPU forces GPU use
	err = clGetDeviceIDs(platform, CL_DEVICE_TYPE_DEFAULT, 1, &device, &available);
	if (err != CL_SUCCESS) {
		printf("clGetDeviceIDs() failed with %d\n", err);
		exit(1);
	}
	printf("found %d OpenCL devices\n", available);

	props[1] = (cl_context_properties) platform;
	ctx = clCreateContext(props, 1, &device, NULL, NULL, &err);
	if (err != CL_SUCCESS) {
		printf("clCreateContext() failed with %d\n", err);
		exit(1);
	}
	printf("created context\n");

	queue = clCreateCommandQueue(ctx, device, 0, &err);
	if (err != CL_SUCCESS) {
		printf("clCreateCommandQueue() failed with %d\n", err);
		clReleaseContext(ctx);
		exit(1);
	}
	printf("created command queue\n");

	/* Setup clblas. */
	err = clblasSetup();
	if (err != CL_SUCCESS) {
		printf("clblasSetup() failed with %d\n", err);
		clReleaseCommandQueue(queue);
		clReleaseContext(ctx);
		exit(1);
	}
	printf("setup clblas\n");

	/* Prepare OpenCL memory objects and place matrices inside them. */
	bufA = clCreateBuffer(ctx, CL_MEM_READ_ONLY, M * K * sizeof(*A), NULL, &err);
	bufB = clCreateBuffer(ctx, CL_MEM_READ_ONLY, K * N * sizeof(*B), NULL, &err);
	bufC = clCreateBuffer(ctx, CL_MEM_READ_WRITE, M * N * sizeof(*C), NULL, &err);
	printf("created buffers\n");

	err = clEnqueueWriteBuffer(queue, bufA, CL_TRUE, 0, M * K * sizeof(double), A, 0, NULL, NULL);
	err = clEnqueueWriteBuffer(queue, bufB, CL_TRUE, 0, K * N * sizeof(double), B, 0, NULL, NULL);
	err = clEnqueueWriteBuffer(queue, bufC, CL_TRUE, 0, M * N * sizeof(double), C, 0, NULL, NULL);
	printf("enqueud buffers\n");

	/*
	 * Call clblas extended function. Perform gemm for the lower right
	 * sub-matrices
	 */
	err = clblasDgemm(order, transA, transB, M, N, K,
                         alpha, bufA, 0, lda,
                         bufB, 0, ldb, beta,
                         bufC, 0, ldc,
                         1, &queue, 0, NULL, &event);
	if (err != CL_SUCCESS) {
		printf("clblasSgemmEx() failed with %d\n", err);
		ret = 1;
	} else {
		printf("no errors for calculation\n");
		fflush(stdout);
		/* Wait for calculations to be finished. */
		err = clWaitForEvents(1, &event);

		/* Fetch results of calculations from GPU memory. */
		err = clEnqueueReadBuffer(queue, bufC, CL_TRUE, 0,
					  M * N * sizeof(double),
					  C, 0, NULL, NULL);
		printf("got result\n");
	}

	/* Release OpenCL memory objects. */
	clReleaseMemObject(bufC);
	clReleaseMemObject(bufB);
	clReleaseMemObject(bufA);
	printf("released\n");

	/* Finalize work with clblas. */
	clblasTeardown();
	printf("teardown clblas\n");

	/* Release OpenCL working objects. */
	clReleaseCommandQueue(queue);
	printf("release command queue\n");
	clReleaseContext(ctx);
	printf("release context\n");
}