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