/********************************************
 *  frtvsr.c
 *
 *  by Steven L. Jacques, Dec. 28, 1998.
 *
 *  Computes a snapshot of the spatial resolution of fluence rate F(r,t)
 *  at a specific time  t = 0 +/- multiples of k'r/w
 *  in response to a sinusoidally modulated istropic point source S(t) [W]
 *        S(t) =  So( 1 + Mo sin(wt) )
 *  within a uniform unbounded medium with optical properties of 
 *  absorption, reduced scattering, and refractive index.
 *
 *  Uses equations equivalent to the equations published by 
 *  Fishkin et al., Proc. Photo-Opt. Instrum. Eng. 1431:122-135, 1991,
 *  Schmitt, Knuettel, Knutson, JOSA 9:1832-1843, 1992,
 *  Fishkin and Gratton, JOSA 10:127-140, 1993.
 *
 **********/

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

#define	PI     3.1415926
#define CC     2.997925e10

int main() {

/* 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  dr;              /* incremental step size [cm] */
double	rmax;            /* maximum position considered [cm] */

double	f;               /* modulation frequency [Hz] */
double	w;               /* angular frequency of modulation [radians/s] */
double	th;              /* angle per photon lifetime */

double  So;              /* steady-state power of isotropic source point [W] */
double  Mo;              /* modulation of source power [-] */
double  S;               /* time-resolved power of source  [W] */
double	Tss;             /* steady-state transport to observation point [1/cm^2] */
double	M;               /* magnitude of modulation at observation point [-] */
double	P;               /* phase of modulation at observation point [radians] */
double  F;               /* fluence rate at observation point [W/cm^2] */

double	N;               /* number of points in array */
double	t;               /* time [s] */

/* dummy variables */
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) [-] */

rmax        = 100;     /* max distance between source and observation point [cm] */
So          = 1.0;     /* steady-state power of isotropic source point [W] */
Mo          = 1.0;     /* modulation of source power [-] */
f           = 400e6;   /* modulation frequency [Hz] */
N           = 200;	   /* set number of points in array */
t           = 0;       /* time of observation [s] */


/**** RUN
 * Calculate Tss, M, and P,  save F(r,t), M*sin(w*t - 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] */
w = 2.0*PI*f;                           /* angular frequency [1/s] */

b   = sqrt( sqrt( 1.0 + w*w/(mua*mua*c*c) ) ); /* temporary variable b [-]  */
th  = (w/(mua*c));                      /* reference angle theta [radians] */


	/**** PRINTOUT
	 * Create file called "frtvsr.out"
	 * with labeled columns for plotting by user's graphics software.
	 * Columns =  time t [ns]    F(r,t) [W/cm^2]    M*sin(w*t - P)
	 *****/

	target = fopen("frtvsr.out", "w");                /* open file */
	fprintf(target, " r [cm] \t F [W/cm2] \t M*sin(w*t - P) \n"); /* print column labels */

	for (i=0; i<=N; i++) {                        /* do loop */
       r = (i + 0.5)*dr;                         /* time t */
	   Tss = exp(-r/delta)/(4.0*PI*D*r);         /* steady-state transport [1/cm^2] */
	   M   = exp(-r/delta*( b*cos(th/2) -1) );   /* modulation magnitude M [-] */
	   P   = r/delta*b*sin(th/2);                /* phase lag P [radians] */
       S = So*( 1.0 + Mo*sin(w*t) );             /* S(t) */
       F = So*Tss*( 1 + Mo*M*sin(w*t - P) );     /* F(r,t) */
       fprintf(target, "%5.4f \t%5.4e \t%5.4f \n", r, F, M*sin(w*t - P)); /* print next row */
	   }

	fclose(target);                               /* close file */

/***** end RUN ****/


} /* end of main */