Skip to content

Commit bceb4d2

Browse files
janakaranajanakarana
committed
mus and musp values in V1.0.6 and V1.0.7 are off by a factor of 2 for real refractive index input. This fix correct that bug.
1 parent be44d9d commit bceb4d2

4 files changed

Lines changed: 75 additions & 57 deletions

File tree

‎src/calculate.cpp‎

Lines changed: 47 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -18,7 +18,7 @@ void calculate::DoSimulation(Ui_MainWindow *ui, parameters *para,
1818
double tempForwardHalf, tempBackwardHalf;
1919
double sumForward, sumBackward;
2020
double rad;
21-
double factor;
21+
double singleScatCross;
2222
double refRelRe = 1.0;
2323
double refRelIm = 0.0;
2424

@@ -57,23 +57,36 @@ void calculate::DoSimulation(Ui_MainWindow *ui, parameters *para,
5757
for (int r = 0; r < para->nRadius; r++)
5858
{
5959
rad = para->radArray[r];
60+
xPara = k * rad;
61+
6062
refRelRe = para->scatRefRealArray[r] / para->medRef;
6163
refRelIm = para->scatRefImagArray[r] / para->medRef;
6264

63-
xPara = k * rad;
64-
for (int t = 0; t < para->nTheta; t++)
65+
if (refRelIm == 0.0) //FarFieldSolutionForRealRefIndex is ~2x faster than FarFieldSolutionForComplexRefIndex
6566
{
66-
mu = cos(para->minTheta + t * para->stepTheta);
67-
if (refRelIm == 0.0)
67+
for (int t = 0; t < para->nTheta; t++)
68+
{
69+
mu = cos(para->minTheta + t * para->stepTheta);
6870
sim.FarFieldSolutionForRealRefIndex(&cS1, &cS2, &qSca, xPara, refRelRe, mu);
69-
else
71+
curS1[t] = cS1;
72+
curS2[t] = cS2;
73+
}
74+
}
75+
else
76+
{
77+
for (int t = 0; t < para->nTheta; t++)
78+
{
79+
mu = cos(para->minTheta + t * para->stepTheta);
7080
sim.FarFieldSolutionForComplexRefIndex(&cS1, &cS2, &qSca, xPara,
7181
std::complex<double>(refRelRe,-refRelIm), mu); //multiply by -1 to use "n-ik" convention
72-
curS1[t] = cS1;
73-
curS2[t] = cS2;
82+
curS1[t] = cS1;
83+
curS2[t] = cS2;
84+
}
7485
}
86+
87+
singleScatCross = 2.0 * M_PI * rad * rad * qSca;
7588
//Mus calculation
76-
tempMus = 1e-6 *(M_PI * rad * rad * qSca) * para->numDensityArray[r]; //1e-6--> 1micron2 to 1mm2
89+
tempMus = singleScatCross * para->numDensityArray[r] * 1e-6; //1e-6--> 1micron2 to 1mm2
7790
sumMus += tempMus; //Σμs
7891
//G calculation
7992
tempG = CalculateG(curS1, curS2, para);
@@ -86,17 +99,17 @@ void calculate::DoSimulation(Ui_MainWindow *ui, parameters *para,
8699
sumBackward += tempBackwardHalf * tempMus;
87100

88101
//S1 and S2
89-
factor = M_PI * xPara * xPara * qSca;
102+
double factor = k*k*singleScatCross;
90103
for (int t = 0; t < para->nTheta; t++)
91104
{
92105
sumS1[t] += curS1[t] * tempMus;
93106
sumS2[t] += curS2[t] * tempMus;
94-
sumPhaseFuncAve[t] += 0.5 * (util.ComplexAbsSquared(curS1[t])+util.ComplexAbsSquared(curS2[t])) * tempMus / factor ;
107+
sumPhaseFuncAve[t] += 0.5 * (util.ComplexAbsSquared(curS1[t])+util.ComplexAbsSquared(curS2[t])) * tempMus / factor;
95108
sumPhaseFuncPara[t] += util.ComplexAbsSquared(curS2[t]) * tempMus / factor;
96109
sumPhaseFuncPerp[t] += util.ComplexAbsSquared(curS1[t]) * tempMus / factor;
97110
}
98111
}
99-
para->scatCross[w] = 1e6 * sumMus / sumNumDen ; //1e6--> 1mm2 to 1micron2
112+
para->scatCross[w] = sumMus * 1e6 / sumNumDen ; //1e6--> 1mm2 to 1micron2
100113
para->mus[w] = sumMus;
101114
para->g[w] = sumMusG /sumMus;
102115
para->forward[w] = sumForward*100.0/(sumForward+sumBackward); //Not necessary to divide by sumMus in ratio calculation
@@ -118,23 +131,36 @@ void calculate::DoSimulation(Ui_MainWindow *ui, parameters *para,
118131
for (int r = 0; r < para->nRadius; r++)
119132
{
120133
rad = para->radArray[r];
134+
xPara = k * rad;
135+
121136
refRelRe = para->scatRefRealArray[r] / para->medRef;
122137
refRelIm = para->scatRefImagArray[r] / para->medRef;
123138

124-
xPara = k * rad;
125-
for (int t = 0; t < para->nTheta; t++)
139+
if (refRelIm == 0.0) //FarFieldSolutionForRealRefIndex is ~2x faster than FarFieldSolutionForComplexRefIndex
126140
{
127-
mu = cos(para->minTheta + t * para->stepTheta);
128-
if (refRelIm == 0.0)
141+
for (int t = 0; t < para->nTheta; t++)
142+
{
143+
mu = cos(para->minTheta + t * para->stepTheta);
129144
sim.FarFieldSolutionForRealRefIndex(&cS1, &cS2, &qSca, xPara, refRelRe, mu);
130-
else
145+
curS1[t] = cS1;
146+
curS2[t] = cS2;
147+
}
148+
}
149+
else
150+
{
151+
for (int t = 0; t < para->nTheta; t++)
152+
{
153+
mu = cos(para->minTheta + t * para->stepTheta);
131154
sim.FarFieldSolutionForComplexRefIndex(&cS1, &cS2, &qSca, xPara,
132-
std::complex<double>(refRelRe,-refRelIm), mu); //multiply by -1 to use "n-ik" convention
133-
curS1[t] = cS1;
134-
curS2[t] = cS2;
155+
std::complex<double>(refRelRe,-refRelIm), mu); //multiply by -1 to use "n-ik" convention
156+
curS1[t] = cS1;
157+
curS2[t] = cS2;
158+
}
135159
}
160+
161+
singleScatCross = 2.0 * M_PI * rad * rad * qSca;
136162
//Mus calculation
137-
tempMus = 1e-6 *(M_PI * rad * rad * qSca)*para->numDensityArray[r]; //1e-6--> 1micron2 to 1mm2
163+
tempMus = singleScatCross * para->numDensityArray[r] * 1e-6; //1e-6--> 1micron2 to 1mm2
138164
sumMus += tempMus; //Σμs
139165
//G calculation
140166
tempG = CalculateG(curS1, curS2, para);

‎src/mainwindow.ui‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -23,7 +23,7 @@
2323
</size>
2424
</property>
2525
<property name="windowTitle">
26-
<string>Mie Simulator v1.0.7: Virtual Photonics Technology Initiative - Beckman Laser Institute </string>
26+
<string>Mie Simulator v1.0.8: Virtual Photonics Technology Initiative - Beckman Laser Institute </string>
2727
</property>
2828
<property name="windowIcon">
2929
<iconset resource="MieRes.qrc">
@@ -2384,7 +2384,7 @@
23842384
</connection>
23852385
</connections>
23862386
<buttongroups>
2387-
<buttongroup name="buttonGroup_ConcVolFrac"/>
23882387
<buttongroup name="buttonGroup_MonoPoly"/>
2388+
<buttongroup name="buttonGroup_ConcVolFrac"/>
23892389
</buttongroups>
23902390
</ui>

‎src/miesimulation.cpp‎

Lines changed: 25 additions & 33 deletions
Original file line numberDiff line numberDiff line change
@@ -29,11 +29,11 @@ MieSimulation::~MieSimulation(void)
2929
void MieSimulation::FarFieldSolutionForRealRefIndex(std::complex<double> *cS1,std::complex<double> *cS2,
3030
double *qSca,double xPara,double relRef,double mu)
3131
{
32-
double xstop;
33-
double mx;
32+
double xstop;
3433
double dPCost0, dPCost1;
3534
double fac0, fac1, fac2, fac3, fac4;
36-
double j_x0, j_x1, y_x0, y_x1;
35+
double j_x0, j_x1, y_x0, y_x1;
36+
double mx;
3737
double dervDn1, dervDn2;
3838
double sumQsca;
3939
std::complex<double> an, bn;
@@ -61,7 +61,7 @@ void MieSimulation::FarFieldSolutionForRealRefIndex(std::complex<double> *cS1,st
6161

6262
Dn_mx[nmx-1]=0;
6363
for (int N = nmx-1; N>0; N--)
64-
Dn_mx[N-1] = (N/mx)-(1.0/(Dn_mx[N]+N/mx));
64+
Dn_mx[N-1] = (double(N)/mx)-(1.0/(Dn_mx[N]+double(N)/mx));
6565

6666
// Legendre Polynomials
6767
dPCost0 = 0.0; //pi0
@@ -109,8 +109,8 @@ void MieSimulation::FarFieldSolutionForRealRefIndex(std::complex<double> *cS1,st
109109

110110
// Calculate an and bn (According to Bohren and Huffman book)
111111
// Remark: GouGouesbet uses size parameter as "ka" instead of "kx"
112-
dervDn1 = Dn_mx[n]/m +n/x;
113-
dervDn2 = m*Dn_mx[n] +n/x;
112+
dervDn1 = (Dn_mx[n]/m) + (double(n)/x);
113+
dervDn2 = (m*Dn_mx[n]) + (double(n)/x);
114114

115115
an = (dervDn1*j_x[n]-j_x[n-1])/ (dervDn1*xi_x[n]-xi_x[n-1]);
116116
bn = (dervDn2*j_x[n]-j_x[n-1])/ (dervDn2*xi_x[n]-xi_x[n-1]);
@@ -145,22 +145,14 @@ void MieSimulation::FarFieldSolutionForRealRefIndex(std::complex<double> *cS1,st
145145
void MieSimulation::FarFieldSolutionForComplexRefIndex(std::complex<double> *cS1, std::complex<double> *cS2,
146146
double *qSca, double xPara, std::complex<double> cRelRef,
147147
double mu)
148-
{
149-
int n;
150-
int nstop, ymod, nmx, arraySize;
151-
double xstop;
152-
std::complex<double> mx;
148+
{
149+
double xstop;
153150
double dPCost0, dPCost1;
154151
double fac0, fac1, fac2, fac3, fac4;
155152
double j_x0, j_x1, y_x0, y_x1;
156-
double *j_x;
157-
double *y_x;
158-
double *piCost;
159-
double *tauCost;
160-
std::complex<double> *Dn_mx;
153+
std::complex<double> mx;
161154
std::complex<double> dervDn1, dervDn2;
162155
double sumQSca;
163-
std::complex<double> *xi_x;
164156
std::complex<double> an, bn;
165157
std::complex<double> tempS1, tempS2;
166158

@@ -172,21 +164,21 @@ void MieSimulation::FarFieldSolutionForComplexRefIndex(std::complex<double> *cS1
172164

173165
mx = m*x;
174166
xstop = x+4.0*(pow(x,(1.0/3.0)))+2.0;
175-
nstop = ceil(xstop);
176-
ymod = ceil(util.ComplexAbs(mx));
177-
nmx = max(xstop,ymod)+15;
178-
arraySize = nstop+1;
179-
180-
Dn_mx = new std::complex<double> [nmx];
181-
j_x = new double [arraySize];
182-
y_x = new double [arraySize];
183-
piCost = new double [arraySize];
184-
tauCost = new double [arraySize];
185-
xi_x = new std::complex<double> [arraySize];
167+
int nstop = ceil(xstop);
168+
int ymod = ceil(util.ComplexAbs(mx));
169+
int nmx = max(xstop,ymod)+15;
170+
int arraySize = nstop+1;
171+
172+
std::complex<double> *Dn_mx = new std::complex<double> [nmx];
173+
double *j_x = new double [arraySize];
174+
double *y_x = new double [arraySize];
175+
double *piCost = new double [arraySize];
176+
double *tauCost = new double [arraySize];
177+
std::complex<double> *xi_x = new std::complex<double> [arraySize];
186178

187179
Dn_mx[nmx-1] = 0;
188-
for (n = nmx-1; n>0; n--)
189-
Dn_mx[n-1] = (double(n)/mx)-(1.0/(Dn_mx[n]+(double(n)/mx)));
180+
for (int N = nmx-1; N>0; N--)
181+
Dn_mx[N-1] = (double(N)/mx)-(1.0/(Dn_mx[N]+double(N)/mx));
190182

191183
// Legendre Polynomials
192184
dPCost0 = 0.0; //pi0
@@ -206,7 +198,7 @@ void MieSimulation::FarFieldSolutionForComplexRefIndex(std::complex<double> *cS1
206198
tempS2 = std::complex<double> (0.0, 0.0);
207199
sumQSca = 0.0;
208200

209-
n=1;
201+
int n=1;
210202
while ((n-1-nstop) < 0)
211203
{
212204
fac0 = double(n); // n
@@ -239,7 +231,7 @@ void MieSimulation::FarFieldSolutionForComplexRefIndex(std::complex<double> *cS1
239231

240232
an = (dervDn1*j_x[n]-j_x[n-1])/ (dervDn1*xi_x[n]-xi_x[n-1]);
241233
bn = (dervDn2*j_x[n]-j_x[n-1])/ (dervDn2*xi_x[n]-xi_x[n-1]);
242-
sumQSca += fac2*(util.ComplexAbs(an)*util.ComplexAbs(an) + util.ComplexAbs(bn)*util.ComplexAbs(bn));
234+
sumQSca += (fac2/(x*x))*(util.ComplexAbs(an)*util.ComplexAbs(an) + util.ComplexAbs(bn)*util.ComplexAbs(bn));
243235

244236
// Calculate cS1 and cS2
245237
tempS1 += fac4*(an*piCost[n-1]+bn*tauCost[n-1]);
@@ -248,7 +240,7 @@ void MieSimulation::FarFieldSolutionForComplexRefIndex(std::complex<double> *cS1
248240
n = n+1;
249241
}
250242

251-
*qSca= (2.0/(x*x))*sumQSca;
243+
*qSca= sumQSca;
252244
*cS1 = tempS1;
253245
*cS2 = tempS2;
254246

‎src/utilities.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -17,7 +17,7 @@ void utilities::Delay()
1717
{
1818
QTime dieTime= QTime::currentTime().addMSecs(1);
1919
while( QTime::currentTime() < dieTime )
20-
QCoreApplication::processEvents(QEventLoop::AllEvents, 100);
20+
QCoreApplication::processEvents(QEventLoop::AllEvents, 50);
2121
}
2222

2323
//Intensity calculation (Amplitude^2)

0 commit comments

Comments
 (0)