da3d4e0e
dmayerich
Initial commit.
|
1
2
3
4
5
6
7
8
9
10
|
#include <math.h>
#include <complex>
#include <iostream>
#include <fstream>
#include "globals.h"
//#include "cufft.h"
using namespace std;
#define pi 3.14159
|
52a5fe9d
dmayerich
Added double supp...
|
11
|
typedef complex<double> scComplex;
|
da3d4e0e
dmayerich
Initial commit.
|
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
|
int cbessjyva(double v,complex<double> z,double &vm,complex<double>*cjv,
complex<double>*cyv,complex<double>*cjvp,complex<double>*cyvp);
int bessjyv(double v,double x,double &vm,double *jv,double *yv,
double *djv,double *dyv);
complex<double> Jl_neg(complex<double> x)
{
//this function computes the bessel function of the first kind Jl(x) for l = -0.5
return ( sqrt(2.0/pi) * cos(x) )/sqrt(x);
}
double Jl_neg(double x)
{
//this function computes the bessel function of the first kind Jl(x) for l = -0.5
return ( sqrt(2.0/pi) * cos(x) )/sqrt(x);
}
double Yl_neg(double x)
{
//this function computes the bessel function of the second kind Yl(x) for l = -0.5;
return ( sqrt(2.0/pi) * sin(x) )/sqrt(x);
}
|
52a5fe9d
dmayerich
Added double supp...
|
36
|
void computeB(complex<double>* B, double radius, complex<double> refIndex, double lambda, int Nl)
|
da3d4e0e
dmayerich
Initial commit.
|
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
|
{
double k = (2*pi)/lambda;
int b = 2;
//allocate space for the real bessel functions
double* jv = (double*)malloc(sizeof(double)*(Nl+b));
double* yv = (double*)malloc(sizeof(double)*(Nl+b));
double* jvp = (double*)malloc(sizeof(double)*(Nl+b));
double* yvp = (double*)malloc(sizeof(double)*(Nl+b));
//allocate space for the complex bessel functions
complex<double>* cjv = (complex<double>*)malloc(sizeof(complex<double>)*(Nl+b));
complex<double>* cyv = (complex<double>*)malloc(sizeof(complex<double>)*(Nl+b));
complex<double>* cjvp = (complex<double>*)malloc(sizeof(complex<double>)*(Nl+b));
complex<double>* cyvp = (complex<double>*)malloc(sizeof(complex<double>)*(Nl+b));
double kr = k*radius;
complex<double> knr = k*refIndex*(double)radius;
complex<double> n = refIndex;
//compute the bessel functions for k*r
double vm;// = Nl - 1;
bessjyv((Nl)+0.5, kr, vm, jv, yv, jvp, yvp);
//cout<<"Nl: "<<Nl<<" vm: "<<vm<<endl;
//printf("Nl: %f, vm: %f\n", (float)Nl, (float)vm);
//compute the bessel functions for k*n*r
cbessjyva((Nl)+0.5, knr, vm, cjv, cyv, cjvp, cyvp);
//scale factor for spherical bessel functions
double scale_kr = sqrt(pi/(2.0*kr));
complex<double> scale_knr = sqrt(pi/(2.0*knr));
complex<double> numer, denom;
double j_kr;
double y_kr;
complex<double> j_knr;
complex<double> j_d_knr;
double j_d_kr;
complex<double> h_kr;
complex<double> h_d_kr;
complex<double> h_neg;
complex<double> h_pos;
//cout<<"B coefficients:"<<endl;
for(int l=0; l<Nl; l++)
{
//compute the spherical bessel functions
j_kr = jv[l] * scale_kr;
y_kr = yv[l] * scale_kr;
j_knr = cjv[l] * scale_knr;
//compute the Hankel function
h_kr = complex<double>(j_kr, y_kr);
//compute the derivatives
if(l == 0)
{
//spherical bessel functions for l=0
j_d_kr = scale_kr * (Jl_neg(kr) - (jv[l] + kr*jv[l+1])/kr )/2.0;
j_d_knr = scale_knr * ( Jl_neg(knr) - (cjv[l] + knr*cjv[l+1])/knr )/2.0;
h_neg = complex<double>(scale_kr*Jl_neg(kr), scale_kr*Yl_neg(kr));
h_pos = complex<double>(scale_kr*jv[l+1], scale_kr*yv[l+1]);
h_d_kr = (h_neg - (h_kr + kr*h_pos)/kr)/2.0;
}
else
{
//spherical bessel functions
j_d_kr = scale_kr * (jv[l-1] - (jv[l] + kr*jv[l+1])/kr )/2.0;
j_d_knr = scale_knr * ( cjv[l-1] - (cjv[l] + knr*cjv[l+1])/knr )/2.0;
h_neg = complex<double>(scale_kr*jv[l-1], scale_kr*yv[l-1]);
h_pos = complex<double>(scale_kr*jv[l+1], scale_kr*yv[l+1]);
h_d_kr = (h_neg - (h_kr + kr*h_pos)/kr)/2.0;
}
numer = j_kr*j_d_knr*n - j_knr*j_d_kr;
denom = j_knr*h_d_kr - h_kr*j_d_knr*n;
|
52a5fe9d
dmayerich
Added double supp...
|
114
|
B[l] = numer/denom;
|
da3d4e0e
dmayerich
Initial commit.
|
115
|
|
52a5fe9d
dmayerich
Added double supp...
|
116
|
//B[l] = scComplex(temp.real(), temp.imag());
|
da3d4e0e
dmayerich
Initial commit.
|
117
118
119
120
121
122
123
124
125
126
127
128
129
|
//cout<<B[l]<<endl;
}
free(jv);
free(yv);
free(jvp);
free(yvp);
free(cjv);
free(cyv);
free(cjvp);
free(cyvp);
}
|
52a5fe9d
dmayerich
Added double supp...
|
130
|
void Legendre(double* P, double x, int Nl)
|
da3d4e0e
dmayerich
Initial commit.
|
131
132
133
134
135
136
137
138
139
140
141
142
|
{
//computes the legendre polynomials from orders 0 to Nl-1
P[0] = 1;
if(Nl == 1) return;
P[1] = x;
for(int l = 2; l < Nl; l++)
{
P[l] = ((2*l - 1)*x*P[l-1] - (l - 1)*P[l-2])/l;
}
}
|
52a5fe9d
dmayerich
Added double supp...
|
143
|
complex<double> integrateUi(double cAngleI, double cAngleO, double oAngleI, double oAngleO, double M = 2*pi)
|
da3d4e0e
dmayerich
Initial commit.
|
144
145
146
147
148
149
150
151
152
|
{
/*This function integrates the incident field of magnitude M in the far zone
in order to evaluate the field at the central pixel of a detector.
cNAi = condenser inner angle
cNAo = condenser outer angle
oNAi = objective inner angle
oNAo = objective outer angle
M = field magnitude*/
|
52a5fe9d
dmayerich
Added double supp...
|
153
154
|
double alphaIn = max(cAngleI, oAngleI);
double alphaOut = min(cAngleO,oAngleO);
|
da3d4e0e
dmayerich
Initial commit.
|
155
|
|
52a5fe9d
dmayerich
Added double supp...
|
156
|
complex<double> Ui;
|
da3d4e0e
dmayerich
Initial commit.
|
157
|
if(alphaIn > alphaOut)
|
52a5fe9d
dmayerich
Added double supp...
|
158
|
Ui = complex<double>(0.0, 0.0);
|
da3d4e0e
dmayerich
Initial commit.
|
159
|
else
|
52a5fe9d
dmayerich
Added double supp...
|
160
|
Ui = complex<double>(M * 2 * pi * (cos(alphaIn) - cos(alphaOut)), 0.0f);
|
da3d4e0e
dmayerich
Initial commit.
|
161
162
163
164
165
|
return Ui;
}
|
52a5fe9d
dmayerich
Added double supp...
|
166
|
void computeCondenserAlpha(double* alpha, int Nl, double cAngleI, double cAngleO)
|
da3d4e0e
dmayerich
Initial commit.
|
167
168
169
170
171
172
173
|
{
/*This function computes the condenser integral in order to build the field of incident light
alpha = list of Nl floating point values representing the condenser alpha as a function of l
Nl = number of orders in the incident field
cAngleI, cAngleO = inner and outer condenser angles (inner and outer NA)*/
//compute the Legendre polynomials for the condenser aperature
|
52a5fe9d
dmayerich
Added double supp...
|
174
|
double* PcNAo = (double*)malloc(sizeof(double)*(Nl+1));
|
da3d4e0e
dmayerich
Initial commit.
|
175
|
Legendre(PcNAo, cos(cAngleO), Nl+1);
|
52a5fe9d
dmayerich
Added double supp...
|
176
|
double* PcNAi = (double*)malloc(sizeof(double)*(Nl+1));
|
da3d4e0e
dmayerich
Initial commit.
|
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
|
Legendre(PcNAi, cos(cAngleI), Nl+1);
for(int l=0; l<Nl; l++)
{
//integration term
if(l == 0)
alpha[l] = -PcNAo[l+1] + PcNAo[0] + PcNAi[l+1] - PcNAi[0];
else
alpha[l] = -PcNAo[l+1] + PcNAo[l-1] + PcNAi[l+1] - PcNAi[l-1];
alpha[l] *= 2 * pi;
}
}
|
52a5fe9d
dmayerich
Added double supp...
|
192
193
|
complex<double> integrateUs(double r, double lambda, complex<double> eta,
double cAngleI, double cAngleO, double oAngleI, double oAngleO, double M = 2*pi)
|
da3d4e0e
dmayerich
Initial commit.
|
194
195
196
197
198
199
200
201
202
203
204
205
206
|
{
/*This function integrates the incident field of magnitude M in the far zone
in order to evaluate the field at the central pixel of a detector.
r = sphere radius
lambda = wavelength
eta = index of refraction
cNAi = condenser inner NA
cNAo = condenser outer NA
oNAi = objective inner NA
oNAo = objective outer NA
M = field magnitude*/
//compute the required number of orders
|
52a5fe9d
dmayerich
Added double supp...
|
207
208
|
double k = 2*pi/lambda;
int Nl = (int)ceil( k + 4 * exp(log(k*r)/3) + 3 );
|
da3d4e0e
dmayerich
Initial commit.
|
209
210
|
//compute the material coefficients B
|
52a5fe9d
dmayerich
Added double supp...
|
211
|
complex<double>* B = (complex<double>*)malloc(sizeof(complex<double>)*Nl);
|
da3d4e0e
dmayerich
Initial commit.
|
212
|
//compute the Legendre polynomials for the condenser and objective aperatures
|
52a5fe9d
dmayerich
Added double supp...
|
213
|
double* PcNAo = (double*)malloc(sizeof(double)*(Nl+1));
|
da3d4e0e
dmayerich
Initial commit.
|
214
|
Legendre(PcNAo, cos(cAngleO), Nl+1);
|
52a5fe9d
dmayerich
Added double supp...
|
215
|
double* PcNAi = (double*)malloc(sizeof(double)*(Nl+1));
|
da3d4e0e
dmayerich
Initial commit.
|
216
217
|
Legendre(PcNAi, cos(cAngleI), Nl+1);
|
52a5fe9d
dmayerich
Added double supp...
|
218
|
double* PoNAo = (double*)malloc(sizeof(double)*(Nl+1));
|
da3d4e0e
dmayerich
Initial commit.
|
219
|
Legendre(PoNAo, cos(oAngleO), Nl+1);
|
52a5fe9d
dmayerich
Added double supp...
|
220
|
double* PoNAi = (double*)malloc(sizeof(double)*(Nl+1));
|
da3d4e0e
dmayerich
Initial commit.
|
221
222
223
224
225
226
227
228
229
|
Legendre(PoNAi, cos(oAngleI), Nl+1);
//store the index of refraction;
complex<double> IR(eta.real(), eta.imag());
//compute the scattering coefficients
computeB(B, r, IR, lambda, Nl);
//aperature terms for the condenser (alpha) and objective (beta)
|
52a5fe9d
dmayerich
Added double supp...
|
230
231
232
233
|
double alpha;
double beta;
double c;
complex<double> Us(0.0, 0.0);
|
da3d4e0e
dmayerich
Initial commit.
|
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
|
for(int l=0; l<Nl; l++)
{
//integration term
if(l == 0)
{
alpha = -PcNAo[l+1] + PcNAo[0] + PcNAi[l+1] - PcNAi[0];
beta = -PoNAo[l+1] + PoNAo[0] + PoNAi[l+1] - PoNAi[0];
}
else
{
alpha = -PcNAo[l+1] + PcNAo[l-1] + PcNAi[l+1] - PcNAi[l-1];
beta = -PoNAo[l+1] + PoNAo[l-1] + PoNAi[l+1] - PoNAi[l-1];
}
c = (2*pi)/(2.0 * l + 1.0);
Us += c * alpha * beta * B[l] * M;
}
free(PcNAo);
free(PcNAi);
free(PoNAo);
free(PoNAi);
free(B);
return Us;
}
void pointSpectrum()
{
PD.StartTimer(SIMULATE_SPECTRUM);
//clear the previous spectrum
SimSpectrum.clear();
|
52a5fe9d
dmayerich
Added double supp...
|
269
270
|
double dNu = 2.0f;
double lambda;
|
da3d4e0e
dmayerich
Initial commit.
|
271
272
|
//compute the angles based on NA
|
52a5fe9d
dmayerich
Added double supp...
|
273
274
275
276
|
double cAngleI = asin(cNAi);
double cAngleO = asin(cNAo);
double oAngleI = asin(oNAi);
double oAngleO = asin(oNAo);
|
da3d4e0e
dmayerich
Initial commit.
|
277
278
279
280
281
282
283
284
285
286
287
288
289
290
|
//implement a reflection-mode system if necessary
if(opticsMode == ReflectionOpticsType){
//set the condenser to match the objective
cAngleI = oAngleI;
cAngleO = oAngleO;
//invert the objective
oAngleO = pi - cAngleI;
oAngleI = pi - cAngleO;
}
//integrate the incident field at the detector position
|
52a5fe9d
dmayerich
Added double supp...
|
291
292
|
complex<double> Ui = integrateUi(cAngleI, cAngleO, oAngleI, oAngleO, 2*pi);
double I0 = Ui.real() * Ui.real() + Ui.imag() * Ui.imag();
|
da3d4e0e
dmayerich
Initial commit.
|
293
294
295
296
|
I0 *= scaleI0;
|
52a5fe9d
dmayerich
Added double supp...
|
297
|
//double I;
|
da3d4e0e
dmayerich
Initial commit.
|
298
|
SpecPair temp;
|
52a5fe9d
dmayerich
Added double supp...
|
299
300
301
|
double nu;
complex<double> eta;
complex<double> Us, U;
|
da3d4e0e
dmayerich
Initial commit.
|
302
|
|
52a5fe9d
dmayerich
Added double supp...
|
303
304
|
double vecLen = 0.0;
for(unsigned int i=0; i<EtaK.size(); i++)
|
da3d4e0e
dmayerich
Initial commit.
|
305
306
307
308
|
{
nu = EtaK[i].nu;
lambda = 10000.0f/nu;
if(applyMaterial)
|
52a5fe9d
dmayerich
Added double supp...
|
309
|
eta = complex<double>(EtaN[i].A, EtaK[i].A);
|
da3d4e0e
dmayerich
Initial commit.
|
310
|
else
|
52a5fe9d
dmayerich
Added double supp...
|
311
|
eta = complex<double>(baseIR, 0.0);
|
da3d4e0e
dmayerich
Initial commit.
|
312
313
314
315
316
|
//integrate the scattered field at the detector position
Us = integrateUs(radius, lambda, eta, cAngleI, cAngleO, oAngleI, oAngleO, 2*pi);
U = Us + Ui;
|
52a5fe9d
dmayerich
Added double supp...
|
317
|
double I = U.real() * U.real() + U.imag() * U.imag();
|
da3d4e0e
dmayerich
Initial commit.
|
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
|
temp.nu = nu;
//set the spectrum value based on the current display type
if(dispSimType == AbsorbanceSpecType)
temp.A = -log10(I/I0);
else
temp.A = I;
if(dispNormalize)
vecLen += temp.A * temp.A;
SimSpectrum.push_back(temp);
}
vecLen = sqrt(vecLen);
if(dispNormalize)
|
52a5fe9d
dmayerich
Added double supp...
|
335
|
for(unsigned int i=0; i<SimSpectrum.size(); i++)
|
da3d4e0e
dmayerich
Initial commit.
|
336
337
338
339
340
|
SimSpectrum[i].A = (SimSpectrum[i].A / vecLen) * dispNormFactor;
PD.EndTimer(SIMULATE_SPECTRUM);
}
|
52a5fe9d
dmayerich
Added double supp...
|
341
|
void updateSpectrum(double* I, double I0, int n)
|
da3d4e0e
dmayerich
Initial commit.
|
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
|
{
SimSpectrum.clear();
SpecPair temp;
//update the displayed spectrum based on the computed intensity I
for(int i=0; i<n; i++)
{
temp.nu = EtaK[i].nu;
//set the spectrum value based on the current display type
if(dispSimType == AbsorbanceSpecType)
temp.A = -log10(I[i]/I0);
else
temp.A = I[i];
SimSpectrum.push_back(temp);
}
}
|
52a5fe9d
dmayerich
Added double supp...
|
361
|
void computeCassegrainAngles(double& cAngleI, double& cAngleO, double& oAngleI, double& oAngleO)
|
da3d4e0e
dmayerich
Initial commit.
|
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
|
{
//compute the angles based on NA
cAngleI = asin(cNAi);
cAngleO = asin(cNAo);
oAngleI = asin(oNAi);
oAngleO = asin(oNAo);
//implement a reflection-mode system if necessary
if(opticsMode == ReflectionOpticsType){
//set the condenser to match the objective
cAngleI = oAngleI;
cAngleO = oAngleO;
//invert the objective
oAngleO = pi - cAngleI;
oAngleI = pi - cAngleO;
}
}
int computeNl()
{
|
52a5fe9d
dmayerich
Added double supp...
|
386
387
388
389
|
double maxNu = EtaK.back().nu;
double maxLambda = 10000.0f/maxNu;
double k = 2*pi/maxLambda;
int Nl = (int)ceil( k + 4 * exp(log(k*radius)/3) + 3 );
|
da3d4e0e
dmayerich
Initial commit.
|
390
391
392
393
|
return Nl;
}
|
52a5fe9d
dmayerich
Added double supp...
|
394
|
void computeBArray(complex<double>* B, int Nl, int nLambda)
|
da3d4e0e
dmayerich
Initial commit.
|
395
|
{
|
52a5fe9d
dmayerich
Added double supp...
|
396
397
398
|
double nu;
complex<double> eta;
double* Lambda = (double*)malloc(sizeof(double) * nLambda);
|
da3d4e0e
dmayerich
Initial commit.
|
399
400
|
//for each wavenumber nu
|
52a5fe9d
dmayerich
Added double supp...
|
401
|
for(unsigned int i=0; i<EtaK.size(); i++)
|
da3d4e0e
dmayerich
Initial commit.
|
402
403
404
405
406
|
{
//compute information based on wavelength and material
nu = EtaK[i].nu;
Lambda[i] = 10000.0f/nu;
if(applyMaterial)
|
52a5fe9d
dmayerich
Added double supp...
|
407
|
eta = complex<double>(EtaN[i].A, EtaK[i].A);
|
da3d4e0e
dmayerich
Initial commit.
|
408
|
else
|
52a5fe9d
dmayerich
Added double supp...
|
409
|
eta = complex<double>(baseIR, 0.0);
|
da3d4e0e
dmayerich
Initial commit.
|
410
411
412
413
414
415
416
417
418
|
//allocate memory for the scattering coefficients
//complex<float>* B = (complex<float>*)malloc(sizeof(complex<float>)*Nl);
complex<double> IR(eta.real(), eta.imag());
computeB(&B[i * Nl], radius, IR, Lambda[i], Nl);
}
}
|
52a5fe9d
dmayerich
Added double supp...
|
419
|
void computeOpticalParameters(double& cAngleI, double& cAngleO, double& oAngleI, double& oAngleO, double& I0, double* alpha, int Nl)
|
da3d4e0e
dmayerich
Initial commit.
|
420
421
422
423
424
|
{
computeCassegrainAngles(cAngleI, cAngleO, oAngleI, oAngleO);
//evaluate the incident field intensity
I0 = 0.0;
|
52a5fe9d
dmayerich
Added double supp...
|
425
|
complex<double> Ui;
|
da3d4e0e
dmayerich
Initial commit.
|
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
|
Ui = integrateUi(cAngleI, cAngleO, oAngleI, oAngleO, 2*pi);
I0 = Ui.real()*2*pi;
//compute alpha (condenser integral)
computeCondenserAlpha(alpha, Nl, cAngleI, cAngleO);
}
void gpuDetectorSpectrum(int numSamples)
{
//integrate across the objective aperature and calculate the resulting intensity on a detector
PD.StartTimer(SIMULATE_SPECTRUM);
//clear the previous spectrum
SimSpectrum.clear();
//compute Nl (maximum order of the spectrum)
int Nl = computeNl();
|
52a5fe9d
dmayerich
Added double supp...
|
444
445
|
double* alpha = (double*)malloc(sizeof(double)*(Nl + 1));
double cAngleI, cAngleO, oAngleI, oAngleO, I0;
|
da3d4e0e
dmayerich
Initial commit.
|
446
447
448
449
450
451
|
computeOpticalParameters(cAngleI, cAngleO, oAngleI, oAngleO, I0, alpha, Nl);
//allocate space for a list of wavelengths
int nLambda = EtaK.size();
//allocate space for the 2D array (Nl x nu) of scattering coefficients
|
52a5fe9d
dmayerich
Added double supp...
|
452
|
complex<double>* B = (complex<double>*)malloc(sizeof(complex<double>) * Nl * nLambda);
|
da3d4e0e
dmayerich
Initial commit.
|
453
454
455
456
|
computeBArray(B, Nl, nLambda);
//allocate temporary space for the spectrum
|
52a5fe9d
dmayerich
Added double supp...
|
457
|
double* I = (double*)malloc(sizeof(double) * EtaK.size());
|
da3d4e0e
dmayerich
Initial commit.
|
458
459
460
|
//compute the spectrum on the GPU
PD.StartTimer(SIMULATE_GPU);
|
52a5fe9d
dmayerich
Added double supp...
|
461
|
cudaComputeSpectrum(I, (double*)B, alpha, Nl, nLambda, oAngleI, oAngleO, cAngleI, cAngleO, numSamples);
|
da3d4e0e
dmayerich
Initial commit.
|
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
|
PD.EndTimer(SIMULATE_GPU);
updateSpectrum(I, I0, nLambda);
PD.EndTimer(SIMULATE_SPECTRUM);
}
void SimulateSpectrum()
{
if(pointDetector)
pointSpectrum();
else
gpuDetectorSpectrum(objectiveSamples);
//detectorSpectrum(objectiveSamples);
}
|
52a5fe9d
dmayerich
Added double supp...
|
479
|
double absorbanceDistortion(){
|
da3d4e0e
dmayerich
Initial commit.
|
480
481
|
//compute the mean of the spectrum
|
52a5fe9d
dmayerich
Added double supp...
|
482
483
|
double sumSim = 0.0;
for(unsigned int i=0; i<SimSpectrum.size(); i++)
|
da3d4e0e
dmayerich
Initial commit.
|
484
485
486
|
{
sumSim += SimSpectrum[i].A;
}
|
52a5fe9d
dmayerich
Added double supp...
|
487
|
double meanSim = sumSim/SimSpectrum.size();
|
da3d4e0e
dmayerich
Initial commit.
|
488
489
|
//compute the distortion (MSE from the mean)
|
52a5fe9d
dmayerich
Added double supp...
|
490
491
|
double sumSE = 0.0;
for(unsigned int i=0; i<SimSpectrum.size(); i++)
|
da3d4e0e
dmayerich
Initial commit.
|
492
493
494
|
{
sumSE += pow(SimSpectrum[i].A - meanSim, 2);
}
|
52a5fe9d
dmayerich
Added double supp...
|
495
|
double MSE = sumSE / SimSpectrum.size();
|
da3d4e0e
dmayerich
Initial commit.
|
496
497
498
499
|
return MSE;
}
|
52a5fe9d
dmayerich
Added double supp...
|
500
|
double intensityDistortion(){
|
da3d4e0e
dmayerich
Initial commit.
|
501
502
|
//compute the magnitude of the spectrum
|
52a5fe9d
dmayerich
Added double supp...
|
503
504
|
double sumSim = 0.0;
for(unsigned int i=0; i<SimSpectrum.size(); i++)
|
da3d4e0e
dmayerich
Initial commit.
|
505
506
507
|
{
sumSim += SimSpectrum[i].A * SimSpectrum[i].A;
}
|
52a5fe9d
dmayerich
Added double supp...
|
508
|
double magSim = sqrt(sumSim);
|
da3d4e0e
dmayerich
Initial commit.
|
509
510
|
//compute the distortion (MSE from the mean)
|
52a5fe9d
dmayerich
Added double supp...
|
511
512
|
double sumSE = 0.0;
for(unsigned int i=0; i<SimSpectrum.size(); i++)
|
da3d4e0e
dmayerich
Initial commit.
|
513
514
515
|
{
sumSE += (SimSpectrum[i].A/magSim) * (1.0/SimSpectrum.size());
}
|
52a5fe9d
dmayerich
Added double supp...
|
516
|
double MSE = sumSE;
|
da3d4e0e
dmayerich
Initial commit.
|
517
518
519
520
521
522
523
524
|
return MSE;
}
void MinimizeDistortion(){
ofstream outFile("distortion.txt");
//set the parameters for the distortion simulation
|
52a5fe9d
dmayerich
Added double supp...
|
525
|
double step = 0.001;
|
da3d4e0e
dmayerich
Initial commit.
|
526
527
528
529
530
531
532
533
|
oNAi = 0.2;
oNAo = 0.5;
//compute the optical parameters
//compute Nl (maximum order of the spectrum)
int Nl = computeNl();
|
52a5fe9d
dmayerich
Added double supp...
|
534
535
|
double* alpha = (double*)malloc(sizeof(double)*(Nl + 1));
double cAngleI, cAngleO, oAngleI, oAngleO, I0;
|
da3d4e0e
dmayerich
Initial commit.
|
536
537
538
539
540
|
//allocate space for a list of wavelengths
int nLambda = EtaK.size();
//allocate temporary space for the spectrum
|
52a5fe9d
dmayerich
Added double supp...
|
541
|
double* I = (double*)malloc(sizeof(double) * EtaK.size());
|
da3d4e0e
dmayerich
Initial commit.
|
542
543
544
|
//calculate the material parameters
//allocate space for the 2D array (Nl x nu) of scattering coefficients
|
52a5fe9d
dmayerich
Added double supp...
|
545
|
complex<double>* B = (complex<double>*)malloc(sizeof(complex<double>) * Nl * nLambda);
|
da3d4e0e
dmayerich
Initial commit.
|
546
547
548
549
|
computeBArray(B, Nl, nLambda);
|
52a5fe9d
dmayerich
Added double supp...
|
550
551
552
|
double D;
double e = 0.001;
for(double i=0.0; i<=oNAo-step; i+=step)
|
da3d4e0e
dmayerich
Initial commit.
|
553
554
|
{
|
52a5fe9d
dmayerich
Added double supp...
|
555
|
for(double o=oNAi+step; o<=1.0; o+=step)
|
da3d4e0e
dmayerich
Initial commit.
|
556
557
558
559
560
561
562
563
564
565
566
|
{
//set the current optical parameters
cNAi = i;
cNAo = o;
//compute the optical parameters
computeOpticalParameters(cAngleI, cAngleO, oAngleI, oAngleO, I0, alpha, Nl);
//simulate the spectrum
|
52a5fe9d
dmayerich
Added double supp...
|
567
|
cudaComputeSpectrum(I, (double*)B, alpha, Nl, nLambda, oAngleI, oAngleO, cAngleI, cAngleO, objectiveSamples);
|
da3d4e0e
dmayerich
Initial commit.
|
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
|
updateSpectrum(I, I0, nLambda);
if(dispSimType == AbsorbanceSpecType)
{
if(i + e >= o || i + e >= oNAo || oNAi + e >= o || oNAi + e >= oNAo)
D = 0.0;
else
D = absorbanceDistortion();
}
else
{
if(i >= o || oNAi >= oNAo)
D=0;
else
D = intensityDistortion();
}
outFile<<D<<" ";
}
outFile<<endl;
cout<<i<<endl;
}
outFile.close();
}
|