/********************************************
 *  ssdt.c
 *  Steady-state diffusion theory
 *
 *  by Steven L. Jacques, Dec. 28, 1998.
 *
 *  Computes the steady-state fluence rate F(r) [W/cm^2]
 *  in response to an istropic point source of power Po(t) [W]
 *  deposited at the origin
 *  within a uniform unbounded medium with optical properties of 
 *  absorption, reduced scattering, and refractive index.
 *
 *  Creates output file "ssdt.out" which records F(r).
 *
 **********/

#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  Po;              /* power of isotropic source point [W] */
double	T;               /* transport to observation point [1/cm^2] */
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) [-] */

Po          = 1.0;     /* power of source [W] */
rmax        = 2;       /* max distance between source and observation point [cm] */
N           = 100;	   /* set number of points in array */


/**** 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 "ssdt.out"
	 * with labeled columns for plotting by user's graphics software.
	 * Columns =  r [cm]    F(r)
	 *****/

	target = fopen("ssdt.out", "w");      /* open file */
	/* print column labels */
	fprintf(target, " r [cm] \t F(r) [W/cm^2] \n"); 

	/* 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( -fabs(r)/delta )/(4*PI*mua*delta*delta*fabs(r));     /* T(r) */
          F = Po*T;                       /* F(r) */
          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( -fabs(r)/delta )/(4*PI*mua*delta*delta*fabs(r));     /* T(r) */
          F = Po*T;                       /* F(r) */
          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 */