@@ -14,13 +14,15 @@ void calculate::DoSimulation(Ui_MainWindow *ui, parameters *para,
1414 double * mus1000, double *g1000)
1515{
1616 MieSimulation sim;
17- double sumMus, sumMusG, tempMus, tempG ;
18- double tempForwardHalf, tempBackwardHalf ;
17+ double curMus ;
18+ double sumMus, sumMusG ;
1919 double sumForward, sumBackward;
20- double rad;
21- double singleScatCross;
20+ double sumCsca, sumCext, sumCback;
21+ double tempPhase;
22+ double piRadiusSquared;
2223 double refRelRe = 1.0 ;
2324 double refRelIm = 0.0 ;
25+ double xPara;
2426
2527 std::complex <double > *curS1 = new std::complex <double > [para->nTheta ];
2628 std::complex <double > *curS2 = new std::complex <double > [para->nTheta ];
@@ -46,6 +48,9 @@ void calculate::DoSimulation(Ui_MainWindow *ui, parameters *para,
4648 sumMusG = 0.0 ;
4749 sumForward = 0.0 ;
4850 sumBackward = 0.0 ;
51+ sumCsca = 0.0 ;
52+ sumCext = 0.0 ;
53+ sumCback = 0.0 ;
4954 for (int t = 0 ; t < para->nTheta ; t++)
5055 {
5156 sumS1[t] = std::complex <double >(0.0 ,0.0 );
@@ -56,8 +61,8 @@ void calculate::DoSimulation(Ui_MainWindow *ui, parameters *para,
5661 }
5762 for (int r = 0 ; r < para->nRadius ; r++)
5863 {
59- rad = para->radArray [r];
60- xPara = k * rad ;
64+ xPara = k * para->radArray [r];
65+ piRadiusSquared = M_PI * para-> radArray [r] * para-> radArray [r] ;
6166
6267 refRelRe = para->scatRefRealArray [r] / para->medRef ;
6368 refRelIm = para->scatRefImagArray [r] / para->medRef ;
@@ -67,7 +72,7 @@ void calculate::DoSimulation(Ui_MainWindow *ui, parameters *para,
6772 for (int t = 0 ; t < para->nTheta ; t++)
6873 {
6974 mu = cos (para->minTheta + t * para->stepTheta );
70- sim.FarFieldSolutionForRealRefIndex (&cS1, &cS2, &qSca, xPara, refRelRe, mu);
75+ sim.FarFieldSolutionForRealRefIndex (&cS1, &cS2, &qSca, &qExt, &qBack, xPara, refRelRe, mu);
7176 curS1[t] = cS1;
7277 curS2[t] = cS2;
7378 }
@@ -77,40 +82,50 @@ void calculate::DoSimulation(Ui_MainWindow *ui, parameters *para,
7782 for (int t = 0 ; t < para->nTheta ; t++)
7883 {
7984 mu = cos (para->minTheta + t * para->stepTheta );
80- sim.FarFieldSolutionForComplexRefIndex (&cS1, &cS2, &qSca, xPara,
85+ sim.FarFieldSolutionForComplexRefIndex (&cS1, &cS2, &qSca, &qExt, &qBack, xPara,
8186 std::complex <double >(refRelRe,-refRelIm), mu); // multiply by -1 to use "n-ik" convention
8287 curS1[t] = cS1;
8388 curS2[t] = cS2;
8489 }
8590 }
86-
87- singleScatCross = 2.0 * M_PI * rad * rad * qSca;
91+ // Ref: Schmitt and Kumar, Applied Optics 37(13) 1998
8892 // Mus calculation
89- tempMus = singleScatCross * para->numDensityArray [r] * 1e-6 ; // 1e-6--> 1micron2 to 1mm2
90- sumMus += tempMus ; // Σμs
93+ curMus = piRadiusSquared * qSca * para->numDensityArray [r];
94+ sumMus += curMus ; // Σμs
9195 // G calculation
92- tempG = CalculateG (curS1, curS2, para);
93- sumMusG += tempG * tempMus;
96+ sumMusG += CalculateG (curS1, curS2, para) * curMus;
9497 // Forward Scattering
95- tempForwardHalf = CalculateForwardBackward (curS1, curS2, para, 0 , (para->nTheta -1 )/2 );
96- sumForward += tempForwardHalf * tempMus;
98+ sumForward += CalculateForwardBackward (curS1, curS2, para, 0 , (para->nTheta -1 )/2 )* curMus;
9799 // Backward Scattering
98- tempBackwardHalf = CalculateForwardBackward (curS1, curS2, para, ((para->nTheta -1 )/2 ), para->nTheta );
99- sumBackward += tempBackwardHalf * tempMus;
100+ sumBackward += CalculateForwardBackward (curS1, curS2, para, ((para->nTheta -1 )/2 ), para->nTheta )* curMus;
101+ // Coefficients
102+ sumCsca += piRadiusSquared * qSca * para->numDensityArray [r];
103+ sumCext += piRadiusSquared * qExt * para->numDensityArray [r];
104+ sumCback += piRadiusSquared * qBack * para->numDensityArray [r];
100105
101106 // S1 and S2
102- double factor = k*k*singleScatCross ;
107+ double factor = 1.0 / (M_PI * xPara * xPara * qSca); // 1/(pi *X^2 * qSca) ;
103108 for (int t = 0 ; t < para->nTheta ; t++)
104109 {
105- sumS1[t] += curS1[t] * tempMus;
106- sumS2[t] += curS2[t] * tempMus;
107- sumPhaseFuncAve[t] += 0.5 * (util.ComplexAbsSquared (curS1[t])+util.ComplexAbsSquared (curS2[t])) * tempMus / factor;
108- sumPhaseFuncPara[t] += util.ComplexAbsSquared (curS2[t]) * tempMus / factor;
109- sumPhaseFuncPerp[t] += util.ComplexAbsSquared (curS1[t]) * tempMus / factor;
110+ sumS1[t] += curS1[t] * curMus;
111+ sumS2[t] += curS2[t] * curMus;
112+ // Phase function - Parallel
113+ tempPhase = util.ComplexAbsSquared (curS2[t])*factor; // |S2|^2/(pi X^2 Qsca)
114+ sumPhaseFuncPara[t] += tempPhase*curMus;
115+ // Phase function - Perpendicular
116+ tempPhase = util.ComplexAbsSquared (curS1[t])*factor; // |S1|^2/(pi X^2 Qsca)
117+ sumPhaseFuncPerp[t] += tempPhase*curMus;
118+ // Phase function - Average
119+ tempPhase = 0.5 *(util.ComplexAbsSquared (curS1[t])+util.ComplexAbsSquared (curS2[t]))*factor; // (|S1|^2+|S2|^2)/(2 pi X^2 Qsca)
120+ sumPhaseFuncAve[t] += tempPhase*curMus;
110121 }
111122 }
112- para->scatCross [w] = sumMus * 1e6 / sumNumDen ; // 1e6--> 1mm2 to 1micron2
113- para->mus [w] = sumMus;
123+ // Normalize
124+ para->cSca [w] = sumCsca / sumNumDen ;
125+ para->cExt [w] = sumCext / sumNumDen;
126+ para->cBack [w] = sumCback / sumNumDen;
127+ para->SizePara [w] = xPara;
128+ para->mus [w] = sumMus * 1e-6 ; // 1e-6--> 1micron2 to 1mm2
114129 para->g [w] = sumMusG /sumMus;
115130 para->forward [w] = sumForward*100.0 /(sumForward+sumBackward); // Not necessary to divide by sumMus in ratio calculation
116131 para->backward [w] = sumBackward*100.0 /(sumForward+sumBackward); // Not necessary to divide by sumMus in ratio calculation
@@ -130,8 +145,8 @@ void calculate::DoSimulation(Ui_MainWindow *ui, parameters *para,
130145 sumMusG = 0.0 ;
131146 for (int r = 0 ; r < para->nRadius ; r++)
132147 {
133- rad = para->radArray [r];
134- xPara = k * rad ;
148+ xPara = k * para->radArray [r];
149+ piRadiusSquared = M_PI * para-> radArray [r] * para-> radArray [r] ;
135150
136151 refRelRe = para->scatRefRealArray [r] / para->medRef ;
137152 refRelIm = para->scatRefImagArray [r] / para->medRef ;
@@ -141,7 +156,7 @@ void calculate::DoSimulation(Ui_MainWindow *ui, parameters *para,
141156 for (int t = 0 ; t < para->nTheta ; t++)
142157 {
143158 mu = cos (para->minTheta + t * para->stepTheta );
144- sim.FarFieldSolutionForRealRefIndex (&cS1, &cS2, &qSca, xPara, refRelRe, mu);
159+ sim.FarFieldSolutionForRealRefIndex (&cS1, &cS2, &qSca, &qExt, &qBack, xPara, refRelRe, mu);
145160 curS1[t] = cS1;
146161 curS2[t] = cS2;
147162 }
@@ -151,20 +166,17 @@ void calculate::DoSimulation(Ui_MainWindow *ui, parameters *para,
151166 for (int t = 0 ; t < para->nTheta ; t++)
152167 {
153168 mu = cos (para->minTheta + t * para->stepTheta );
154- sim.FarFieldSolutionForComplexRefIndex (&cS1, &cS2, &qSca, xPara,
169+ sim.FarFieldSolutionForComplexRefIndex (&cS1, &cS2, &qSca, &qExt, &qBack, xPara,
155170 std::complex <double >(refRelRe,-refRelIm), mu); // multiply by -1 to use "n-ik" convention
156171 curS1[t] = cS1;
157172 curS2[t] = cS2;
158173 }
159174 }
160-
161- singleScatCross = 2.0 * M_PI * rad * rad * qSca;
162175 // Mus calculation
163- tempMus = singleScatCross * para->numDensityArray [r] * 1e-6 ; // 1e-6--> 1micron2 to 1mm2
164- sumMus += tempMus ; // Σμs
176+ curMus = piRadiusSquared * qSca * para->numDensityArray [r] * 1e-6 ; // 1e-6--> 1micron2 to 1mm2
177+ sumMus += curMus ; // Σμs
165178 // G calculation
166- tempG = CalculateG (curS1, curS2, para);
167- sumMusG += tempG*tempMus;
179+ sumMusG += CalculateG (curS1, curS2, para)*curMus;
168180 }
169181 *mus1000 = sumMus;
170182 *g1000 = sumMusG /sumMus;
@@ -262,7 +274,6 @@ double calculate::CalculateG(std::complex<double> *S1, std::complex<double> *S2,
262274 return num/den;
263275}
264276
265-
266277// Set sphere parameters
267278void calculate::SetSphereRadiusAndRefIndex (parameters *para, int index, bool flagVolOrConc)
268279{
0 commit comments