/*************
 * mcmain.c
 *************/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>

/*************************/
/**** USER CHOICES ******/
/*************************/
#define USER_LABEL  "mcmain.c for class"
#define BINS        101
/*************************/
/*************************/

/**********************
 * DECLARE SUBROUTINES
 *********************/
/* The Monte Carlo subroutine */
void mcsub(double mua, double mus, double g, double n1, double n2, 
           long NR, long NZ, double dr, double dz, double Nphotons,
           int mcflag, double xs, double ys, double zs,
           double radius, double waist, double zfocus,
           double *J, double **F);

/* Computes internal reflectance at tissue/air interface */
double RFresnel(double n1, double n2, double ca1, double *ca2_Ptr);

/* Saves surface escape R(r) and fluence rate distribution F(z,r) */
void SaveFile(double *J, double **F, int Nfile,
              long NR, long NZ, double dr, double dz);

/* Random number generator 
   Initiate by RandomGen(0,1,NULL)
   Use as rnd = RandomGen(1,0,NULL) */
double RandomGen(char Type, long Seed, long *Status); 

/* Memory allocation routines 
 * from MCML ver. 1.0, 1992 L. V. Wang, S. L. Jacques,
 * which are modified versions from Numerical Recipes in C. */
void   nrerror(char error_text[]);
double *AllocVector(short nl, short nh);
double **AllocMatrix(short nrl,short nrh,short ncl,short nch);
void   FreeVector(double *v,short nl,short nh);
void   FreeMatrix(double **m,short nrl,short nrh,short ncl,short nch);


/**********************
 * MAIN PROGRAM
 *********************/
int main() {

/**** Declare variables ******/
double   mua;   /* excitation absorption coeff. [cm^-1] */
double   mus;   /* excitation scattering coeff. [cm^-1] */
double   g;     /* excitation anisotropy [dimensionless] */
double   n1;    /* refractive index of medium */
double   n2;    /* refractive index outside medium */
short    mcflag; /* 0 = collimated, 1 = focused Gaussian, 
                           2 = isotropic pt */
double   radius; /* used if mcflag = 0 or 1 */
double   waist;  /* used if mcflag = 1 */
double   zfocus; /* used if mcflag = 1 */
double   xs;     /* used if mcflag = 2 */
double   ys;     /* used if mcflag = 2 */
double   zs;     /* used if mcflag = 2 */
/* other parameters */
double   Nruns;  /* number photons launched = Nruns x 1e6 */
double   dr;     /* radial bin size [cm] */
double   dz;     /* depth bin size [cm] */

char     label[1];
double   Nphotons;
long     ir, iz;
double   temp; /* dummy variables */
double   start_time, finish_time1; /* for clock() */
double   timeA, timeB;
time_t   now;
double   *Jx;
double   **Fx;
long     NR   = BINS;    /* number of radial bins */
long     NZ   = BINS;    /* number of depth bins */
Jx       = AllocVector(1,BINS);
Fx       = AllocMatrix(1, BINS, 1, BINS); /* for absorbed excitation */

strcpy(label, USER_LABEL);
printf("\n||-------------------------\n");
printf(  "|| %s\n",label);
printf(  "||-------------------------\n\n");

start_time = clock();
now = time(NULL);
printf("%s\n", ctime(&now));

/*************************/
/**** USER CHOICES ******/
/*************************/
mua    = 1.0;  /* excitation absorption coeff. [cm^-1] */
mus    = 100.0;/* excitation scattering coeff. [cm^-1] */
g      = 0.90; /* excitation anisotropy [dimensionless] */
n1     = 1.33; /* refractive index of medium */
n2     = 1.00; /* refractive index outside medium */
mcflag = 1;    /* 0 = collimated, 1 = focused Gaussian, 
                           2 = isotropic pt */
radius = 0.0300;  /* used if mcflag = 0 or 1 */
waist  = 0.0010;  /* used if mcflag = 1 */
zfocus = 0.0300;  /* used if mcflag = 1 */
xs     = 0.0;  /* used if mcflag = 2 */
ys     = 0.0;  /* used if mcflag = 2 */
zs     = 0.090909;  /* used if mcflag = 2 */
/* other parameters */
Nruns  = 0.01;    /* number photons launched = Nruns x 1e6 */
dr     = 0.0100; /* radial bin size [cm] */
dz     = 0.0100; /* depth bin size [cm] */


if (1) { /* Switch printout ON=1 or OFF=0 */
	/* print out summary of parameters to user */
	printf("----- USER CHOICES -----\n");
	printf("mua     = %1.3f\n",mua);
	printf("mus     = %1.3f\n",mus);
	printf("g       = %1.3f\n",g);
	printf("n1      = %1.3f\n",n1);
	printf("n2      = %1.3f\n",n2);
	printf("mcflag  = %d\n"   ,mcflag);
	printf("radius  = %1.4f\n",radius);
	printf("waist   = %1.4f\n",waist);
	printf("zfocus  = %1.4f\n",radius);
	printf("xs      = %1.4f\n",xs);
	printf("ys      = %1.4f\n",ys);
	printf("zs      = %1.4f\n",zs);
	printf("OTHER\n");
	printf("Nruns   = %1.1e\n",Nruns);
	printf("dr      = %1.4e\n",dr);
	printf("dz      = %1.4e\n",dz);
	printf("---------------\n\n");
	}

/* Initialize arrays */
for (ir=1; ir<=NR; ir++) {
	Jx[ir] = 0.0;
	for (iz=1; iz<=NR; iz++) {
		Fx[iz][ir] = 0.0;
		}
	}

/* Time completion estimate */
timeA = clock();
mcsub(	mua, mus, g, n1, n2,   /* CALL THE MONTE CARLO SUBROUTINE */
		NR, NZ, dr, dz, 1000,
		mcflag, xs, ys, zs,
		radius, waist, zfocus,
		Jx, Fx);               /* returns J, F */
timeB = clock();
temp = (timeB - timeA)/CLOCKS_PER_SEC/60/1000; /* min per photon */
printf("%5.3e min/photon \n", temp);
printf("estimated completion times = %5.2f min\n", temp*1e6*Nruns);


/*********************
 * CALL mcsub()
 *********************/
Nphotons = 1e6*Nruns;

mcsub(	mua, mus, g, n1, n2,   /* CALL THE MONTE CARLO SUBROUTINE */
		NR, NZ, dr, dz, Nphotons,
		mcflag, xs, ys, zs,
		radius, waist, zfocus,
		Jx, Fx);               /* returns J, F */

/* Save results to files Ji.dat and Fi.dat, where i = mcflag. */
SaveFile(Jx, Fx, mcflag, NR, NZ, dr, dz); 
printf("------------------------------------------------------\n");

finish_time1 = clock();
printf("--------------------------------------------------------------------\n");
printf("Elapsed Time for excitation = %5.2f min\n", (double)(finish_time1-start_time)/CLOCKS_PER_SEC/60);
now = time(NULL);
printf("%s\n", ctime(&now));


FreeVector(Jx, 1, BINS);
FreeMatrix(Fx, 1, BINS, 1, BINS);
return(0);
}

/**********************************************************
 *                     SUBROUTINES
 **********************************************************/

/**********************************************************
 * The Monte Carlo SUBROUTINE
 **********************************************************/

void mcsub(double mua, double mus, double g, double n1, double n2, 
           long NR, long NZ, double dr, double dz, double Nphotons,
           int mcflag, double xs, double ys, double zs, 
           double radius, double waist, double zfocus,
           double *J, double **F) 
{
/* Constants */
double	PI          = 3.1415926;
short	ALIVE       = 1;           /* if photon not yet terminated */
short	DEAD        = 0;           /* if photon is to be terminated */
double	THRESHOLD   = 0.0001;        /* used in roulette */
double	CHANCE      = 0.1;           /* used in roulette */

/* Variable parameters */
double	mut, albedo, absorb, rsp, Rsptot, Atot;
double	rnd, xfocus;
double	x,y,z, ux,uy,uz,uz1, uxx,uyy,uzz, s,r,W,temp;
double	psi,costheta,sintheta,cospsi,sinpsi;
long   	iphoton, ir, iz, CNT;
short	photon_status;

/**** INITIALIZATIONS *****/
RandomGen(0,1,NULL); /* initiate with seed = 1, or any long integer. */
CNT = 0;
mut    = mua + mus;
albedo = mus/mut;
Rsptot = 0.0; /* accumulate absorbed photon weight */
Atot   = 0.0; /* accumulate specular reflectance per photon */

/* initialize arrays to zero */
for (ir=1; ir<=NR; ir++) {
	J[ir] = 0.0;
	for (iz=1; iz<=NZ; iz++)
  		F[iz][ir] = 0.0;
  	}

/*============================================================
======================= RUN N photons =====================
 * Launch N photons, initializing each one before progation.
============================================================*/
for (iphoton=1; iphoton<=Nphotons; iphoton++) {

	/* Print out progress for user if mcflag < 3 */
	temp = (double)iphoton;
	if ((mcflag < 3) & (temp >= 1000)) {
		if (temp<10000) {
			if (fmod(temp,1000)==0) printf("%1.0f     photons\n",temp);
			}
		else if (temp<100000) {
			if (fmod(temp,10000)==0) printf("%1.0f    photons\n",temp);
			}
		else if (temp<1000000) {
			if (fmod(temp,100000)==0) printf("%1.0f   photons\n",temp);
			}
		else if (temp<10000000) {
			if (fmod(temp,1000000)==0) printf("%1.0f  photons\n",temp);
			}
		else if (temp<100000000) {
			if (fmod(temp,10000000)==0) printf("%1.0f photons\n",temp);
			}
		}
	
/**** LAUNCH 
   Initialize photon position and trajectory.
   Implements an isotropic point source.
*****/

if (mcflag == 0) {
	/* UNIFORM COLLIMATED BEAM INCIDENT AT SURFACE */
	/* Launch at (r,z) = (radius*sqrt(rnd), 0).
	 * Due to cylindrical symmetry, radial launch position is
	 * assigned to x while y = 0. 
	 * radius = radius of uniform beam. */
	/* Initial position */
	rnd = RandomGen(1,0,NULL); 
	x = radius*sqrt(rnd); 
	y = 0;
	z = 0;
	/* Initial trajectory as cosines */
	ux = 0;
	uy = 0;
	uz = 1;  
	/* specular reflectance */
	temp   = n1/n2; /* refractive index mismatch, internal/external */
	temp   = (1.0 - temp)/(1.0 + temp);
	rsp    = temp*temp; /* specular reflectance at boundary */
	}
else if (mcflag == 1) {
	/* GAUSSIAN BEAM AT SURFACE */
	/* Launch at (r,z) = (radius*sqrt(-log(rnd)), 0).
	 * Due to cylindrical symmetry, radial launch position is
	 * assigned to x while y = 0. 
	 * radius = 1/e radius of Gaussian beam at surface. 
	 * waist  = 1/e radius of Gaussian focus.
	 * zfocus = depth of focal point. */
	/* Initial position */
	while ((rnd = RandomGen(1,0,NULL)) <= 0.0); /* avoids rnd = 0 */
	x = radius*sqrt(-log(rnd));
	y = 0.0;
	z = 0.0;
	/* Initial trajectory as cosines */
	/* Due to cylindrical symmetry, radial launch trajectory is
	 * assigned to ux and uz while uy = 0. */
	xfocus = waist/radius*x;
	temp = sqrt((x - xfocus)*(x - xfocus) + zfocus*zfocus);
	sintheta = -(x - xfocus)/temp;
	costheta = zfocus/temp;
	ux = sintheta;
	uy = 0.0;
	uz = costheta;
	/* specular reflectance and diffraction */
	rsp = RFresnel(n2, n1, costheta, &uz); /* new uz */
	ux  = -sqrt(1.0 - uz*uz);         /* new ux */
	}
else  {
	/* ISOTROPIC POINT SOURCE AT POSITION xs,ys,zs */
	/* Initial position */
	x = xs;
	y = ys;
	z = zs;
	/* Initial trajectory as cosines */
	costheta = 1.0 - 2.0*RandomGen(1,0,NULL);
	sintheta = sqrt(1.0 - costheta*costheta);
	psi = 2.0*PI*RandomGen(1,0,NULL);
	cospsi = cos(psi);
	if (psi < PI)
		sinpsi = sqrt(1.0 - cospsi*cospsi); 
	else
		sinpsi = -sqrt(1.0 - cospsi*cospsi);
	ux = sintheta*cospsi;
	uy = sintheta*sinpsi;
	uz = costheta;
	/* specular reflectance */
	rsp = 0.0;
	}

W             = 1.0 - rsp;  /* set photon initial weight */
Rsptot        += rsp; /* accumulate specular reflectance per photon */
photon_status = ALIVE;

/******************************************
****** HOP_ESCAPE_SPINCYCLE **************
* Propagate one photon until it dies by ESCAPE or ROULETTE. 
*******************************************/
do {

/**** HOP
 * Take step to new position
 * s = stepsize
 * ux, uy, uz are cosines of current photon trajectory
 *****/
	while ((rnd = RandomGen(1,0,NULL)) <= 0.0);   /* avoids rnd = 0 */
	s = -log(rnd)/mut;   /* Step size.  Note: log() is base e */
	x += s*ux;           /* Update positions. */
	y += s*uy;
	z += s*uz;

	/* Does photon ESCAPE at surface? ... z <= 0? */
 	if (z <= 0) {
		rnd = RandomGen(1,0,NULL); 
		/* Check Fresnel reflectance at surface boundary */
		if (rnd > RFresnel(n1, n2, -uz, &uz1)) {  
			/* Photon escapes at external angle, uz1 = cos(angle) */
			x -= s*ux;       /* return to original position */
			y -= s*uy;
			z -= s*uz;
			s  = fabs(z/uz); /* calculate stepsize to reach surface */
			x += s*ux;       /* partial step to reach surface */
			y += s*uy;
			r = sqrt(x*x + y*y);   /* find radial position r */
			ir = (long)(r/dr) + 1; /* round to 1 <= ir */
			if (ir > NR) ir = NR;  /* ir = NR is overflow bin */
			J[ir] += W;      /* increment escaping flux */
			photon_status = DEAD;
			}
		else z = -z;  /* Total internal reflection. */
		}

if (photon_status  == ALIVE) {
	/*********************************************
	 ****** SPINCYCLE = DROP_SPIN_ROULETTE ******
	 *********************************************/

	/**** DROP
	 * Drop photon weight (W) into local bin.
	 *****/
	absorb = W*(1 - albedo);  /* photon weight absorbed at this step */
	W -= absorb;              /* decrement WEIGHT by amount absorbed */
	Atot += absorb;           /* accumulate absorbed photon weight */
	/* deposit power in cylindrical coordinates z,r */
	r  = sqrt(x*x + y*y);     /* current cylindrical radial position */
	ir = (long)(r/dr) + 1;        /* round to 1 <= ir */
	iz = (long)(fabs(z)/dz) + 1;  /* round to 1 <= iz */
	if (ir >= NR) ir = NR;        /* last bin is for overflow */
	if (iz >= NZ) iz = NZ;        /* last bin is for overflow */
	F[iz][ir] += absorb;          /* DROP absorbed weight into bin */

	/**** SPIN 
	 * Scatter photon into new trajectory defined by theta and psi.
	 * Theta is specified by cos(theta), which is determined 
	 * based on the Henyey-Greenstein scattering function.
	 * Convert theta and psi into cosines ux, uy, uz. 
	 *****/
	/* Sample for costheta */
	rnd = RandomGen(1,0,NULL);
	if (g == 0.0)
		costheta = 2.0*rnd - 1.0;
	else if (g == 1.0)
		costheta = 1.0;
	else {
		temp = (1.0 - g*g)/(1.0 - g + 2*g*rnd);
		costheta = (1.0 + g*g - temp*temp)/(2.0*g);
		}
	sintheta = sqrt(1.0 - costheta*costheta);/*sqrt faster than sin()*/

	/* Sample psi. */
	psi = 2.0*PI*RandomGen(1,0,NULL);
	cospsi = cos(psi);
	if (psi < PI)
		sinpsi = sqrt(1.0 - cospsi*cospsi); /*sqrt faster */
	else
		sinpsi = -sqrt(1.0 - cospsi*cospsi);

	/* New trajectory. */
	if (1 - fabs(uz) <= 1.0e-12) {  /* close to perpendicular. */
		uxx = sintheta*cospsi;
		uyy = sintheta*sinpsi;
		uzz = costheta*((uz)>=0 ? 1:-1);
		} 
	else {   /* usually use this option */
		temp = sqrt(1.0 - uz*uz);
		uxx = sintheta*(ux*uz*cospsi - uy*sinpsi)/temp + ux*costheta;
		uyy = sintheta*(uy*uz*cospsi + ux*sinpsi)/temp + uy*costheta;
		uzz = -sintheta*cospsi*temp + uz*costheta;
		}

	/* Update trajectory */
	ux = uxx;
	uy = uyy;
	uz = uzz;

	/**** CHECK ROULETTE 
	 * If photon weight below THRESHOLD, then terminate photon using
	 * Roulette technique. Photon has CHANCE probability of having its
	 * weight increased by factor of 1/CHANCE,
	 * and 1-CHANCE probability of terminating.
	 *****/
	if (W < THRESHOLD) {
		rnd = RandomGen(1,0,NULL);
		if (rnd <= CHANCE)
			W /= CHANCE;
		else photon_status = DEAD;
		}

	}/**********************************************
	  **** END of SPINCYCLE = DROP_SPIN_ROULETTE *
	  **********************************************/

} 
while (photon_status == ALIVE);
/******************************************
 ****** END of HOP_ESCAPE_SPINCYCLE ******
 ****** when photon_status == DEAD) ******
 ******************************************/

/* If photon dead, then launch new photon. */
} /*======================= End RUN N photons =====================
=====================================================================*/

/************************
 * NORMALIZE 
 *   J[ir]      escaping flux density [W/cm^2 per W incident] 
 *              where bin = 2.0*PI*r[ir]*dr [cm^2].
 *	  F[iz][ir]  fluence rate [W/cm^2 per W incident] 
 *              where bin = 2.0*PI*r[ir]*dr*dz [cm^3].
 ************************/
temp = 0.0;
for (ir=1; ir<=NR; ir++) {
	r = (ir - 0.5)*dr;
	temp += J[ir];  /* accumulate total escaped photon weight */
	J[ir] /= 2.0*PI*r*dr*Nphotons;                /* flux density */
	for (iz=1; iz<=NZ; iz++)
		F[iz][ir] /= 2.0*PI*r*dr*dz*Nphotons*mua; /* fluence rate */
	}

if (mcflag < 2) {
	printf("Specular = %5.6f\n", Rsptot/Nphotons);
	printf("Absorbed = %5.6f\n", Atot/Nphotons);
	printf("Escaped  = %5.6f\n", temp/Nphotons);
	printf("total    = %5.6f\n", (Rsptot + Atot + temp)/Nphotons);
	}

} /******** END SUBROUTINE **********/


/***********************************************************
 *	FRESNEL REFLECTANCE
 * Computes reflectance as photon passes from medium 1 to 
 * medium 2 with refractive indices n1,n2. Incident
 * angle a1 is specified by cosine value ca1 = cos(a1).
 * Program returns value of transmitted angle a1 as
 * value in *ca2_Ptr = cos(a2).
 ****/
double RFresnel(	double n1,		/* incident refractive index.*/
					double n2,		/* transmit refractive index.*/
					double ca1,		/* cosine of the incident */
									/* angle a1, 0<a1<90 degrees. */
					double *ca2_Ptr) 	/* pointer to the cosine */
									/* of the transmission */
									/* angle a2, a2>0. */
{
double r;

if(n1==n2) { /** matched boundary. **/
	*ca2_Ptr = ca1;
	r = 0.0;
	}
else if(ca1>(1.0 - 1.0e-12)) { /** normal incidence. **/
	*ca2_Ptr = ca1;
	r = (n2-n1)/(n2+n1);
	r *= r;
	}
else if(ca1< 1.0e-6)  {	/** very slanted. **/
	*ca2_Ptr = 0.0;
	r = 1.0;
	}
else  {			  		/** general. **/
	double sa1, sa2; /* sine of incident and transmission angles. */
	double ca2;      /* cosine of transmission angle. */
	sa1 = sqrt(1-ca1*ca1);
	sa2 = n1*sa1/n2;
	if(sa2>=1.0) {	
		/* double check for total internal reflection. */
		*ca2_Ptr = 0.0;
		r = 1.0;
		}
	else {
		double cap, cam;	/* cosines of sum ap or diff am of the two */
							/* angles: ap = a1 + a2, am = a1 - a2. */
		double sap, sam;	/* sines. */
		*ca2_Ptr = ca2 = sqrt(1-sa2*sa2);
		cap = ca1*ca2 - sa1*sa2; /* c+ = cc - ss. */
		cam = ca1*ca2 + sa1*sa2; /* c- = cc + ss. */
		sap = sa1*ca2 + ca1*sa2; /* s+ = sc + cs. */
		sam = sa1*ca2 - ca1*sa2; /* s- = sc - cs. */
		r = 0.5*sam*sam*(cam*cam+cap*cap)/(sap*sap*cam*cam); 
		/* rearranged for speed. */
		}
	}
return(r);
} /******** END SUBROUTINE **********/


/***********************************************************
 * SAVE RESULTS TO FILES 
 *   to files named „Ji.dat¾ and „Fi.dat¾ where i = mcflag.
 * Saves „Ji.dat¾ in following format:
 *   Saves r[ir]  values in first  column, (2:NR,1) = (rows,cols).
 *   Saves Ji[ir] values in second column, (2:NR,2) = (rows,cols).
 * Saves „Fi.dat¾ in following format:
 *   The upper element (1,1) is filled with zero, and ignored.
 *   Saves z[iz] values in first column, (2:NZ,1) = (rows,cols).
 *   Saves r[ir] values in first row,    (1,2:NZ) = (rows,cols).
 *   Saves Fi[iz][ir] in (2:NZ, 2:NR).
***********************************************************/
void SaveFile(	double *R, double **F, int Nfile,
					long NR, long NZ, double dr, double dz)
{
double	r, z, r1, r2;
long  	ir, iz;
char	name[20];
FILE*	target;

/* SAVE flux density J(r) */
sprintf(name, "mcJ%d.dat",Nfile);
target = fopen(name, "w");
for (ir=1; ir<=NR; ir++) {
	r2 = ir*dr;
	r1 = (ir-1)*dr;
	r = 2.0/3*(r2*r2 + r2*r1 + r1*r1)/(r1 + r2);
	fprintf(target, "%5.5f \t%5.12f\n", r, R[ir]);
	}
fclose(target);

/* SAVE fluence rate F(z,r) */
sprintf(name, "mcF%d.dat",Nfile);
target = fopen(name, "w");
fprintf(target, "%5.5f", 0.0); /* ignore upperleft element of matrix */
for (ir=1; ir<=NR; ir++) { 
	r2 = ir*dr;
	r1 = (ir-1)*dr;
	r = 2.0/3*(r2*r2 + r2*r1 + r1*r1)/(r1 + r2);
	fprintf(target, "\t %5.5f", r);
	}
fprintf(target, "\n");
for (iz=1; iz<=NZ; iz++) {
	z = (iz - 0.5)*dz; /* z values for depth position in 1st column */
	fprintf(target, "%5.5f", z);
	for (ir=1; ir<=NR; ir++)
		fprintf(target, "\t %5.12f", F[iz][ir]);
	fprintf(target, "\n");
	}
fclose(target);

} /******** END SUBROUTINE **********/


/*********************************************************************
 *      RANDOM NUMBER GENERATOR
 *      A random number generator that generates uniformly
 *      distributed random numbers between 0 and 1 inclusive.
 *      The algorithm is based on:
 *      W.H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P.
 *      Flannery, "Numerical Recipes in C," Cambridge University
 *      Press, 2nd edition, (1992).
 *      and
 *      D.E. Knuth, "Seminumerical Algorithms," 2nd edition, vol. 2
 *      of "The Art of Computer Programming", Addison-Wesley, (1981).
 *
 *      When Type is 0, sets Seed as the seed. Make sure 0<Seed<32000.
 *      When Type is 1, returns a random number.
 *      When Type is 2, gets the status of the generator.
 *      When Type is 3, restores the status of the generator.
 *
 *      The status of the generator is represented by Status[0..56].
 *      Make sure you initialize the seed before you get random
 *      numbers.
 ****/
#define MBIG 1000000000
#define MSEED 161803398
#define MZ 0
#define FAC 1.0E-9

double RandomGen(char Type, long Seed, long *Status)
{
static long i1, i2, ma[56];   /* ma[0] is not used. */
long        mj, mk;
short       i, ii;

if (Type == 0) { /* set seed. */
	mj = MSEED - (Seed < 0 ? -Seed : Seed);
	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 < MZ)
			mk += MBIG;
		 mj = ma[ii];
		}
	for (ii = 1; ii <= 4; ii++)
		for (i = 1; i <= 55; i++) {
			ma[i] -= ma[1 + (i + 30) % 55];
			if (ma[i] < MZ)
				ma[i] += MBIG;
			}
	i1 = 0;
	i2 = 31;
	} 
else if (Type == 1) {       /* get a number. */
	if (++i1 == 56)
		i1 = 1;
	if (++i2 == 56)
		i2 = 1;
	mj = ma[i1] - ma[i2];
	if (mj < MZ)
		mj += MBIG;
	ma[i1] = mj;
	return (mj * FAC);
	}
else if (Type == 2) {       /* get status. */
	for (i = 0; i < 55; i++)
		Status[i] = ma[i + 1];
	Status[55] = i1;
	Status[56] = i2;
	} 
else if (Type == 3) {       /* restore status. */
	for (i = 0; i < 55; i++)
		ma[i + 1] = Status[i];
	i1 = Status[55];
	i2 = Status[56];
	} 
else
	puts("Wrong parameter to RandomGen().");
return (0);
} 
#undef MBIG
#undef MSEED
#undef MZ
#undef FAC
/*******  end subroutine  ******/


/***********************************************************
 * MEMORY ALLOCATION
 * REPORT ERROR MESSAGE to stderr, then exit the program
 * with signal 1.
 ****/
void nrerror(char error_text[])
{
  fprintf(stderr,"%s\n",error_text);
  fprintf(stderr,"...now exiting to system...\n");
  exit(1);
}

/***********************************************************
 * MEMORY ALLOCATION
 * by Lihong Wang for MCML version 1.0 code, 1992.
 *	ALLOCATE A 1D ARRAY with index from nl to nh inclusive.
 *	Original matrix and vector from Numerical Recipes in C
 *	don't initialize the elements to zero. This will
 *	be accomplished by the following functions. 
 ****/
double *AllocVector(short nl, short nh)
{
double *v;
short i;
v=(double *)malloc((unsigned) (nh-nl+1)*sizeof(double));
if (!v) nrerror("allocation failure in vector()");
v -= nl;
for(i=nl;i<=nh;i++) v[i] = 0.0;	/* init. */
return v;
}

/***********************************************************
 * MEMORY ALLOCATION
 *	ALLOCATE A 2D ARRAY with row index from nrl to nrh 
 *	inclusive, and column index from ncl to nch inclusive.
****/
double **AllocMatrix(short nrl,short nrh,
		short ncl,short nch)
{
short i,j;
double **m;
m=(double **) malloc((unsigned) (nrh-nrl+1) *sizeof(double*));
if (!m) nrerror("allocation failure 1 in matrix()");
	m -= nrl;
for(i=nrl;i<=nrh;i++) {
	m[i]=(double *) malloc((unsigned) (nch-ncl+1) *sizeof(double));
	if (!m[i]) nrerror("allocation failure 2 in matrix()");
	m[i] -= ncl;
	}
for(i=nrl;i<=nrh;i++)
	for(j=ncl;j<=nch;j++) m[i][j] = 0.0;
return m;
}

/***********************************************************
 * MEMORY ALLOCATION
 *	RELEASE MEMORY FOR 1D ARRAY.
 ****/
void FreeVector(double *v,short nl,short nh)
{
free((char*) (v+nl));
}

/***********************************************************
 * MEMORY ALLOCATION
 *	RELEASE MEMORY FOR 2D ARRAY.
 ****/
void FreeMatrix(double **m,short nrl,short nrh, short ncl,short nch)
{
short i;
for(i=nrh;i>=nrl;i--) free((char*) (m[i]+ncl));
free((char*) (m+nrl));
}