/*********************************************************** * Program Name: conv. * * A program used to process the data from mcml - * A Monte Carlo simulation of photon distribution in * multi-layered turbid media in ANSI Standard C. **** * Creation Date: 11/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 * **** * 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 #include #include #define PI 3.1415926 #define STRLEN 256 /* String length. */ #define Boolean char #define SIGN(x) ((x)>=0 ? 1:-1) #define MAX(x, y) ((x)>(y) ? (x):(y)) #define MIN(x, y) ((x)<(y) ? (x):(y)) #define NULLCONVSTRUCT {0, 0, 0, 0, 0} #define NULLOUTSTRUCT {0.0,\ NULL,NULL,NULL,0.0,\ NULL,NULL,NULL,0.0,\ NULL,NULL,NULL,0.0,\ NULL,NULL,NULL,NULL,NULL,\ 0, NULLCONVSTRUCT} /****************** Stuctures *****************************/ /**** * 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. ****/ 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. */ } LayerStruct; /**** * Parameters to describe a photon beam. * Pencil: infinitely narrow beam. This is default for the * beam from the mcml output. * Flat: Flat beam with radius R. * Gaussian: Gaussian with 1/e2 radius R. * Others: general beam described by points with interpolation. ****/ typedef struct { enum { pencil, flat, gaussian, others } type; /* beam type. */ double P; /* total power. [J] */ double R; /* radius. [cm] */ } BeamStruct; /**** * 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. * * For convolution, the grid line separations in z, and alpha * directions are still dz, and da respectively. The numbers * of grid lines in z, and alpha directions are still * nz, and na respectively. However, the grid line separation * and the number of grid lines in r direction are drc and * nrc respectively. ****/ typedef struct { char in_fname[STRLEN]; /* name of mcml output file. */ char in_fformat; /* format of mcml output file . */ /* '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. */ BeamStruct beam; /* incident beam of finite size. */ double drc; /* convolution r grid separation.[cm] */ short nrc; /* convolution array range 0..nrc-1. */ float eps; /* relative error in convolution. */ } InputStruct; /**** * Structure to keep track of what quantities are convolved. * They are initialized to 0, and are set to 1 when the * represented quantity is convolved. * * When the beam profile is changed, they are all reset to 0. ****/ typedef struct { Boolean Rd_ra; Boolean Rd_r; Boolean A_rz; Boolean Tt_ra; Boolean Tt_r; } ConvStruct; /**** * Structures for scored physical quantities * from mcml and to be convolved for photon * beams of finite size. Therefore, "Out" * here means the output of both mcml and conv. * * The member allocated is used to keep the status * of the arrays. It is set to 1 if all the arrays * are allocated and assigned values. It is set to * 0 otherwise. * * 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. [-] */ double **Rd_rac; /* convolved data. [J/(cm2 sr)] */ double *Rd_rc; /* [J/cm2] */ double **A_rzc; /* [J/cm3] */ double **Tt_rac; /* [J/(cm2 sr)] */ double *Tt_rc; /* [J/cm2] */ char allocated; /* set to 1 when arrays are allocated. */ ConvStruct conved; } 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 *); /*********************************************************** * Other prototypes. ****/ void IsoPlot(double **Z, /* the 2D array Z[i][j]. */ long int IXmax, long int IYmax, /* the 0<=i<=IXmax, 0<=j<=IYmax. */ double Dx, double Dy); /* the gridline separations. */ short GetShort(short, short); float GetFloat(float, float);