/********************************************
 *  trdt.c 
 *  Time-resolved diffusion theory
 *
 *  by Steven L. Jacques, Dec. 28, 1998.
 *                        updated July 8, 1999.
 *
 *  Computes the time-resolved fluence rate F(r,t) [W/cm^2]
 *  in response to an impulse istropic point source Uo(t) [J]
 *  deposited at time zero at the origin
 *  within a uniform unbounded medium with optical properties of 
 *  absorption, reduced scattering, and refractive index.
 *
 *  Creates output file "trdt.out" which records snapshots of
 *  F(r,t) vs r at each of 3 timepoints.
 *
 **********/

#include <math.h>
#include <stdio.h>

#define	PI     3.1415926
#define CC     2.997925e10

int main() {

/* problem variables */
double	mua, mus, g, nt; /* optical properties */
double	D, delta;        /* transport parameters */
double	c;               /* speed of light in medium */
double	r;               /* distance from source to observation point */
double	t[3];            /* times of snapshots [s] */

double  Uo;              /* impulse energy of isotropic source point [J] */
double	T;               /* transport to observation point [1/(cm^2 s)] */
double  F;               /* fluence rate at observation point [W/cm^2] */

/* other variables */
double	N;               /* number of points in array */
double	rmax;            /* maximum position considered [cm] */
double  dr;              /* incremental step size [cm] */
short	i, j;            /* dummy indices */
double	b, u;            /* dummy variables */
FILE*	target;          /* pointer to output file */


/**** USER SPECIFIED VALUES
 *****/

mua         = 0.1;     /* absorption coeff. [cm^-1] */
mus         = 100.0;   /* scattering coeff. [cm^-1] */
g           = 0.90;    /* anistoropy [-] */
nt          = 1.33;    /* refractive index of medium (eg., biological tissue) [-] */

Uo      	= 1;       /* energy of impulse [J] */
rmax        = 2;       /* max distance between source and observation point [cm] */
N           = 100;	   /* set number of points in array */
t[0]        = 5e-12;   /* 1st time of observation [s] */
t[1]        = 100e-12; /* 2nd time of observation [s] */
t[2]        = 1e-9;    /* 3rd time of observation [s] */


/**** RUN
 * Calculate Fss, M, and P.
 *****/

dr = rmax/N;                            /* step size [cm] */
c = CC/nt;                              /* speed of light in medium [cm/s] */
D = 1/( 3*mus*(1 - g) );                /* diffusion length [cm] */
delta = sqrt(D/mua);                    /* 1/e optical penetration depth [cm] */


	/**** PRINTOUT
	 * Create file called "trdt.out"
	 * with labeled columns for plotting by user's graphics software.
	 * Columns =  r [cm]    F @t[0]    F @t[1]     F @t[2]
	 *****/

	target = fopen("trdt.out", "w");      /* open file */
	/* print column labels */
	fprintf(target, " r [cm] \t %5.3e [s] \t %5.3e [s] \t %5.3e [s] \n", t[0], t[1], t[2]); 

	/* left half of Gaussian spread */
	for (i=N; i>=0; i--) {                /* do loop */
       r = -(i + 0.5)*dr;                 /* radial position r */
       fprintf(target, "%5.4f", r);       /* print r for this row */
       for (j=0; j<=2; j++) {
          T = exp( -r*r/(4*c*D*t[j]) )/pow(4*PI*c*D*t[j],1.5);     /* T(r,t[j]), [cm^-3] */
          F = c*Uo*T;                       /* F(r,t[j]), [W/cm2] */
          fprintf(target, "\t %5.3e", F); /* print tab and next F value in row */
          }
       fprintf(target,"\n");               /* print carriage return at end of row */
	   }

	/* right half of Gaussian spread */
	for (i=0; i<=N; i++) {                /* do loop */
       r = (i + 0.5)*dr;                  /* radial position r */
       fprintf(target, "%5.4f", r);       /* print r for this row */
       for (j=0; j<=2; j++) {
          T = exp( -r*r/(4*c*D*t[j]) )/pow(4*PI*c*D*t[j],1.5);     /* T(r,t[j]), [cm^-3] */
          F = c*Uo*T;                       /* F(r,t[j]), [W/cm2] */
          fprintf(target, "\t %5.3e", F); /* print tab and next F value in row */
          }
       fprintf(target,"\n");               /* print carriage return at end of row */
	   }

	fclose(target);                               /* close file */
	
/***** end RUN ****/


} /* end of main */