/*********************************************************** * Copyright Univ. of Texas M.D. Anderson Cancer Center * 1992. * * Monte Carlo simulation of photon distribution in * multi-layered turbid media in ANSI Standard C. **** * Starting Date: 10/1991. * Current Date: 6/1992. * * Lihong Wang, Ph. D. * Steven L. Jacques, Ph. D. * Laser Biology Research Laboratory - 17 * M.D. Anderson Cancer Center * University of Texas * 1515 Holcombe Blvd. * Houston, TX 77030 * USA * * This program was based on: * (1) The Pascal code written by Marleen Keijzer and * Steven L. Jacques in this laboratory in 1989, which * deals with multi-layered turbid media. * * (2) Algorithm for semi-infinite turbid medium by * S.A. Prahl, M. Keijzer, S.L. Jacques, A.J. Welch, * SPIE Institute Series Vol. IS 5 (1989), and by * A.N. Witt, The Astrophysical journal Supplement * Series 35, 1-6 (1977). * * Major modifications include: * . Conform to ANSI Standard C. * . Removal of limit on number of array elements, * because arrays in this program are dynamically * allocated. This means that the program can accept * any number of layers or gridlines as long as the * memory permits. * . Avoiding global variables whenever possible. This * program has not used global variables so far. * . Grouping variables logically using structures. * . Top-down design, keep each subroutine clear & * short. * . Reflectance and transmittance are angularly * resolved. **** * General Naming Conventions: * Preprocessor names: all capital letters, * e.g. #define PREPROCESSORS * Globals: first letter of each word is capital, no * underscores, * e.g. short GlobalVar; * Dummy variables: first letter of each word is capital, * and words are connected by underscores, * e.g. void NiceFunction(char Dummy_Var); * Local variables: all lower cases, words are connected * by underscores, * e.g. short local_var; * Function names or data types: same as Globals. * **** * Dimension of length: cm. * ****/ #include #include #include #include #include #include #include #define PI 3.1415926 #define WEIGHT 1E-4 /* Critical weight for roulette. */ #define CHANCE 0.1 /* Chance of roulette survival. */ #define STRLEN 256 /* String length. */ #define Boolean char #define SIGN(x) ((x)>=0 ? 1:-1) /****************** Stuctures *****************************/ /**** * Structure used to describe a photon packet. ****/ typedef struct { double x, y ,z; /* Cartesian coordinates.[cm] */ double ux, uy, uz;/* directional cosines of a photon. */ double w; /* weight. */ Boolean dead; /* 1 if photon is terminated. */ short layer; /* index to layer where the photon */ /* packet resides. */ double s; /* current step size. [cm]. */ double sleft; /* step size left. dimensionless [-]. */ } PhotonStruct; /**** * Structure used to describe the geometry and optical * properties of a layer. * z0 and z1 are the z coordinates for the upper boundary * and lower boundary respectively. * * cos_crit0 and cos_crit1 are the cosines of the * critical angle of total internal reflection for the * upper boundary and lower boundary respectively. * They are set to zero if no total internal reflection * exists. * They are used for computation speed. ****/ typedef struct { double z0, z1; /* z coordinates of a layer. [cm] */ double n; /* refractive index of a layer. */ double mua; /* absorption coefficient. [1/cm] */ double mus; /* scattering coefficient. [1/cm] */ double g; /* anisotropy. */ double cos_crit0, cos_crit1; } LayerStruct; /**** * Input parameters for each independent run. * * z and r are for the cylindrical coordinate system. [cm] * a is for the angle alpha between the photon exiting * direction and the surface normal. [radian] * * The grid line separations in z, r, and alpha * directions are dz, dr, and da respectively. The numbers * of grid lines in z, r, and alpha directions are * nz, nr, and na respectively. * * The member layerspecs will point to an array of * structures which store parameters of each layer. * This array has (number_layers + 2) elements. One * element is for a layer. * The layers 0 and (num_layers + 1) are for top ambient * medium and the bottom ambient medium respectively. ****/ typedef struct { char out_fname[STRLEN]; /* output file name. */ char out_fformat; /* output file format. */ /* 'A' for ASCII, */ /* 'B' for binary. */ long num_photons; /* to be traced. */ double Wth; /* play roulette if photon */ /* weight < Wth.*/ double dz; /* z grid separation.[cm] */ double dr; /* r grid separation.[cm] */ double da; /* alpha grid separation. */ /* [radian] */ short nz; /* array range 0..nz-1. */ short nr; /* array range 0..nr-1. */ short na; /* array range 0..na-1. */ short num_layers; /* number of layers. */ LayerStruct * layerspecs; /* layer parameters. */ } InputStruct; /**** * Structures for scoring physical quantities. * z and r represent z and r coordinates of the * cylindrical coordinate system. [cm] * a is the angle alpha between the photon exiting * direction and the normal to the surfaces. [radian] * See comments of the InputStruct. * See manual for the physcial quantities. ****/ typedef struct { double Rsp; /* specular reflectance. [-] */ double ** Rd_ra; /* 2D distribution of diffuse */ /* reflectance. [1/(cm2 sr)] */ double * Rd_r; /* 1D radial distribution of diffuse */ /* reflectance. [1/cm2] */ double * Rd_a; /* 1D angular distribution of diffuse */ /* reflectance. [1/sr] */ double Rd; /* total diffuse reflectance. [-] */ double ** A_rz; /* 2D probability density in turbid */ /* media over r & z. [1/cm3] */ double * A_z; /* 1D probability density over z. */ /* [1/cm] */ double * A_l; /* each layer's absorption */ /* probability. [-] */ double A; /* total absorption probability. [-] */ double ** Tt_ra; /* 2D distribution of total */ /* transmittance. [1/(cm2 sr)] */ double * Tt_r; /* 1D radial distribution of */ /* transmittance. [1/cm2] */ double * Tt_a; /* 1D angular distribution of */ /* transmittance. [1/sr] */ double Tt; /* total transmittance. [-] */ } OutStruct; /*********************************************************** * Routine prototypes for dynamic memory allocation and * release of arrays and matrices. * Modified from Numerical Recipes in C. ****/ double *AllocVector(short, short); double **AllocMatrix(short, short,short, short); void FreeVector(double *, short, short); void FreeMatrix(double **, short, short, short, short); void nrerror(char *);