/************* * callmcsubG.c * A calling program that * 1. defines parameters for Monte Carlo run * 2. calls the mcsubG() routine * 3. saves the results into an output file using SaveFile() *************/ #include #include #include #include #include #include "mcsubLIBG.h" /*************************/ /**** USER CHOICES *******/ /*************************/ #define BINS 201 /* number of bins, NZ and NR, for z and r */ /*************************/ /*************************/ /********************** * MAIN PROGRAM *********************/ int main(){ //int argc, const char * argv[]) { /*************************/ /**** USER CHOICES *******/ /*************************/ /* number of file for saving */ int Nfile = 1; // thread#, saves as (mcOUT#_#.dat, jrun,Nfile) int jrun = 1; // jth run int i,GKflag; double af,gf,ab,gb,c; /* tissue parameters */ double mua = 1.0; /* excitation absorption coeff. [cm^-1] */ double mus = 100; /* excitation scattering coeff. [cm^-1] */ double g = 0.90; /* excitation anisotropy [dimensionless] */ double n1 = 1.40; /* refractive index of medium */ double n2 = 1.00; /* refractive index outside medium */ /* beam parameters */ short mcflag = 0; /* 0 = collimated, 1 = focused Gaussian, 2 >= isotropic pt. */ double radius = 0; /* used if mcflag = 0 or 1 */ double waist = 0.0; /* used if mcflag = 1 */ double zfocus = 1.0e3; /* used if mcflag = 1 */ double xs = 0.0; /* used if mcflag = 2, isotropic pt source */ double ys = 0.0; /* used if mcflag = 2, isotropic pt source */ double zs = 0.0; /* used if mcflag = 2, isotropic pt source, or mcflag = 0 collimated*/ int boundaryflag = 1; /* 0 = infinite medium, 1 = air/tissue surface boundary */ /* Run parameters */ double Nphotons = 10000; /* number of photons to be launched */ double dr = 0.0100; /* radial bin size [cm] */ double dz = 0.0100; /* depth bin size [cm] */ short PRINTOUT = 1; /*************************/ /****** Setup output parameters, vectors, arrays **********/ double S; /* specular reflectance at air/tissue boundary */ double A; /* total fraction of light absorbed by tissue */ double E; /* total fraction of light escaping tissue */ double E1; // avg mu of escape double* J; /* escaping flux, J[ir], [W/cm2 per W incident] */ double* J1; /* projecting onto z-axis, -uz1*J[ir] */ double** F; /* fluence rate, F[iz][ir], [W/cm2 per W incident] */ short NR = BINS-1; /* number of radial bins */ short NZ = BINS-1; /* number of depth bins */ double albedo, Nsteps, THRESH, t, tperstep; /* Input/Output */ char myname[32]; // Holds the user's choice of myname, used in input and output files. char filename[32]; // temporary filename for writing output. FILE* fid=NULL; // file ID pointer char buf[32]; // buffer for reading header.dat /*************************/ J = AllocVector(1, BINS); /* for escaping flux */ J1 = AllocVector(1, BINS); /* for escaping flux */ F = AllocMatrix(1, BINS, 1, BINS); /* for absorbed fluence rate */ /*************************/ strcpy(filename, "params.dat"); // acquire name from argument of function call by user. fid = fopen(filename,"r"); fgets(buf, 32, fid); sscanf(buf, "%lf", &af); fgets(buf, 32, fid); sscanf(buf, "%lf", &gf); fgets(buf, 32,fid); sscanf(buf, "%lf", &ab); fgets(buf, 32,fid); sscanf(buf, "%lf", &gb); fgets(buf, 32,fid); sscanf(buf, "%lf", &c); fgets(buf, 32,fid); sscanf(buf, "%d", &Nfile); // Nfile: thread # fgets(buf, 32,fid); sscanf(buf, "%d", &jrun); // jth run: =CH fgets(buf, 32,fid); sscanf(buf, "%lf", &mua); fgets(buf, 32,fid); sscanf(buf, "%lf", &mus); fclose(fid); n1 = 1.0; //tissue n2 = 1.0; // air5 dr = 0.0020; dz = 0.0020; g = 0.90; // for HG, ignored if GKflag==1 GKflag = 1; // 1=Gegenbauer, 0=HG. <----- choose ---- printf("af = %0.1f\n",af); printf("gf = %0.4f\n",gf); printf("ab=%0.2f\n",ab); printf("gb=%0.2f\n",gb); printf("c=%0.6f\n",c); printf("thread = %d\n",Nfile); printf("RUN = %d\n",jrun); printf("mua = %0.4f\n",mua); printf("mus = %0.1f\n",mus); // choose Nphotons THRESH = 1e-4; albedo = mus/(mua + mus); Nsteps = log(THRESH)/log(albedo); tperstep = 24.9e-9; t = 60;//2*60*60; // s <----------------------------------- runtime -------------- printf("runtime = %0.2f min\n",t/60); Nphotons = (double)( (long)( t/(tperstep*Nsteps) ) ); printf("Nphotons = %5.4e\n", Nphotons); Nphotons = 1e6; // over-ride previous estimate for Nphotons mcsubG( mua, mus, g, n1, n2, /* CALL THE MONTE CARLO SUBROUTINE */ NR, NZ, dr, dz, Nphotons, mcflag, xs, ys, zs, boundaryflag, radius, waist, zfocus, J, J1, F, &S, &A, &E, &E1, af,gf,ab,gb,c,GKflag,Nfile,jrun, PRINTOUT); /* returns J, F, S, A, E */ if (1==1) { // save to ("mcOUT#_#.dat",jrun,Nfile) SaveFile(Nfile,jrun,J,J1,F, S, A, E, E1, mua, mus, g, n1, n2, mcflag, radius, waist, xs, ys, zs, NR, NZ, dr, dz, Nphotons); } printf("Nphotons = %5.1e\n", Nphotons); printf("Specular = %5.6f\n", S); printf("Absorbed = %5.6f\n", A); printf("Escaped = %5.6f\n", E); printf("total = %5.6f\n", S+A+E); printf("E1 = %5.6f\n", E1); /*************************/ /*************************/ FreeVector(J, 1, BINS); FreeVector(J1, 1, BINS); FreeMatrix(F, 1, BINS, 1, BINS); return(1); }