1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116
|
ir = (long)(r2 / dr) + 1; /* round to 1 <= ir */
iz = (long)(fabs(z) / dz) + 1; /* round to 1 <= iz */
if (ir >= NR) ir = NR; /* last bin is for overflow */
if (iz >= NZ) iz = NZ; /* last bin is for overflow */
F[iz][ir] += absorb; /* DROP absorbed weight into bin */
/**** SPIN
* Scatter photon into new trajectory defined by theta and psi.
* Theta is specified by cos(theta), which is determined
* based on the Henyey-Greenstein scattering function
* Convert theta and psi into cosines ux, uy, uz
*****/
/* Sample for costheta */
if (skinAnis == 0.0)
costheta = 2.0*rnd - 1.0;
else if (skinAnis == 1.0)
costheta = 1.0;
else {
temp = (1.0 - skinAnis*skinAnis) / (1.0 - skinAnis + 2 * skinAnis*rnd);
costheta = (1.0 + skinAnis*skinAnis - temp*temp) / (2.0*skinAnis);
}
sintheta = sqrt(1.0 - costheta*costheta);/*sqrt faster than sin()*/
/* Sample psi. */
psi = 2.0*PI*printRandomDouble(0.0, 1.01);
cospsi = cos(psi);
if (psi < PI)
sinpsi = sqrt(1.0 - cospsi*cospsi); /*sqrt faster */
else
sinpsi = -sqrt(1.0 - cospsi*cospsi);
/* New trajectory. */
if (1 - fabs(uz) <= 1.0e-12) { /* close to perpendicular. */
uxx = sintheta*cospsi;
uyy = sintheta*sinpsi;
uzz = costheta*((uz) >= 0 ? 1 : -1);
}
else { /* usually use this option */
temp = sqrt(1.0 - uz*uz);
uxx = sintheta*(ux*uz*cospsi - uy*sinpsi) / temp + ux*costheta;
uyy = sintheta*(uy*uz*cospsi + ux*sinpsi) / temp + uy*costheta;
uzz = -sintheta*cospsi*temp + uz*costheta;
}
/* Update trajectory */
ux = uxx;
uy = uyy;
uz = uzz;
}
if (rouletteChoice == 1) {
if (photon_weight < THRESHOLD) {
if (rnd <= CHANCE)
photon_weight /= CHANCE;
else photon_condition = DEAD;
}
}
else { //if person chose roulette option 2
int indexSelector1 = printRandomInt(0, monowaveNum);
double percentkeLoss1 = (energyLeft[indexSelector1] / photonE[indexSelector1])*100.0;
if (percentkeLoss1 > 65.0) {
photon_condition = DEAD;
}
else {
if(randomNum > PROB_SURVIVAL) {
photon_condition = DEAD;
}
else {
photon_weight = photon_weight / PROB_SURVIVAL;
}
}
}//rouletteChoice 2
} while (photon_condition == ALIVE);
int a = 0;
rhoSpherVals[a] = distsofarSph;
rhoCylinVals[a] = distsofarCyl;
rhoPlanarVals[a] = distsofarPla;
rhoCounter += 1;
xVal[a] = x;
yVal[a] = y;
zVal[a] = z;
xPlaneDir[a] = ux;
yPlaneDir[a] = uy;
zPlaneDir[a] = uz;
++a;
current_Phot += 1;
p1.totalDistTravelled(distsofarSph, "continuous", "spherical", current_Phot + 1);
p1.totalDistTravelled(distsofarCyl, "continuous", "cylindrical", current_Phot + 1);
p1.totalDistTravelled(distsofarPla, "continuous", "planar", current_Phot + 1);
}while (current_Phot < cr1.photonNum);
cout << "Continuous IR radiation simulation complete." << endl;
free(mieAlbedo);
free(wavelengths);
free(energyLeft);
free(KEMedium);
/************************
* NORMALIZE
* J[ir] escaping flux density [W/cm^2 per W incident]
* where bin = 2.0*PI*r[ir]*dr [cm^2].
* F[iz][ir] fluence rate [W/cm^2 per W incident]
* where bin = 2.0*PI*r[ir]*dr*dz [cm^3].
************************/
temp = 0.0;
for (ir = 0; ir < NR; ir++) {
r1 = (ir - 0.5)*dr;
temp += J[ir]; /* accumulate total escaped photon weight */
J[ir] /= 2.0*PI*r1*dr*cr1.photonNum; /* flux density */
for (iz = 0; iz < NZ; iz++)
F[iz][ir] /= 2.0*PI*r1*dr*dz*cr1.photonNum*totalSkinMua;
/* fluence rate */
}
specularReflec = Rsptot / cr1.photonNum;
tabsorbedFrac = Atot / cr1.photonNum;
tescapingFrac = temp / cr1.photonNum;
Sptr = &specularReflec;
Aptr = &tabsorbedFrac;
Eptr = &tescapingFrac;
/*SaveFile(int Nfile, double *J, double **F, double specularreflec, double absorbingFrac, double escapingFrac,double absorbCoeff, double scattCoeff, double anisG, double n1, double n2,int mcflag, double radius, double waist, double xs, double ys, double zs,int NR, int NZ, double dr, double dz, const int Nphotons);*/
p1.SaveFile(NFile, J,(double**)F, specularReflec, tabsorbedFrac, tescapingFrac, totalSkinMua, avgSkinScattCoeff, skinAnis, cr1.ecm_refrac, cr1.air_refracIND, mcflag, radius, waist, xs, ys, zs, NR, NZ, dr, dz, cr1.photonNum);
| |