3f56f1f9
dmayerich
initial commit
|
1
|
#include "nearfield.h"
|
51b6469a
dmayerich
added look-up tables
|
2
3
4
5
6
7
8
9
10
11
|
#include <time.h>
#include <math.h>
#ifdef _WIN32
#define isnan(x) _isnan(x)
#define isinf(x) (!_finite(x))
#endif
int bessjyv_sph(int v, double z, double &vm, double* cjv,
double* cyv, double* cjvp, double* cyvp);
|
3f56f1f9
dmayerich
initial commit
|
12
13
14
15
16
|
nearfieldStruct::nearfieldStruct()
{
scalarSim = true;
planeWave = false;
|
51b6469a
dmayerich
added look-up tables
|
17
18
|
lut_us = true;
lut_uf = false;
|
3f56f1f9
dmayerich
initial commit
|
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
|
nWaves = 0;
}
void nearfieldStruct::init()
{
//set the field parameters
U.scalarField = scalarSim;
Uf.scalarField = scalarSim;
//initialize dynamic memory
U.init_gpu();
Uf.init_gpu();
}
void nearfieldStruct::destroy()
{
U.kill_gpu();
Uf.kill_gpu();
}
void nearfieldStruct::setPos(bsPoint pMin, bsPoint pMax, bsVector normal)
{
|
3f36b18e
David Mayerich
Adding planewave ...
|
42
|
pos = rts::quad<ptype, 3>(pMin, pMax, normal);
|
3f56f1f9
dmayerich
initial commit
|
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
|
}
void nearfieldStruct::setRes(int x_res, int y_res)
{
U.R[0] = Uf.R[0] = x_res;
U.R[1] = Uf.R[1] = y_res;
}
std::string nearfieldStruct::toStr()
{
std::stringstream ss;
ss<<"------Field Parameters-------"<<std::endl;
ss<<"Wavelength: "<<lambda<<"um"<<std::endl;
ss<<"K Vector (r, theta, phi): "<<k.cart2sph()<<std::endl;
ss<<"Condenser NA: "<<condenser[0]<<" to "<<condenser[1]<<std::endl;
ss<<"Focal Point: "<<focus[0]<<", "<<focus[1]<<", "<<focus[2]<<std::endl;
ss<<"Field Slice: "<<std::endl;
|
51b6469a
dmayerich
added look-up tables
|
61
62
|
if(lut_us)
ss<<"LUT Parameters --- min: "<<d_min<<" max: "<<d_max<<std::endl;
|
3f56f1f9
dmayerich
initial commit
|
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
|
ss<<pos<<std::endl;
ss<<std::endl<<"---------Materials-----------"<<std::endl;
ss<<"Number of Materials: "<<mVector.size()<<std::endl;
ss<<"Refractive Indices at lambda = "<<lambda<<"um"<<std::endl;
//output each material
for(unsigned int m=0; m<mVector.size(); m++)
ss<<" "<<m<<": "<<mVector[m](lambda)<<std::endl;
ss<<"---------Spheres-------------"<<std::endl;
ss<<"Number of Spheres: "<<sVector.size()<<std::endl;
//output each sphere
for(unsigned int s=0; s<sVector.size(); s++)
ss<<sVector[s].toStr()<<std::endl;
|
51b6469a
dmayerich
added look-up tables
|
78
79
80
81
|
ss<<"---------Timings-------------"<<std::endl;
ss<<"Uf = "<<t_Uf<<"ms"<<std::endl;
ss<<"Us = "<<t_Us<<"ms"<<std::endl;
|
3f56f1f9
dmayerich
initial commit
|
82
83
84
85
86
87
88
89
90
|
return ss.str();
}
//generate monte-carlo waves
void nearfieldStruct::calcWaves()
{
inWaves.resize(nWaves);
//re-seed the random number generator
|
51b6469a
dmayerich
added look-up tables
|
91
92
|
//srand(time(NULL));
srand(NULL);
|
3f56f1f9
dmayerich
initial commit
|
93
94
95
96
97
98
99
100
101
102
103
104
105
|
//calculate the monte-carlo samples
mcSampleNA(&inWaves[0], nWaves, k, condenser[0], condenser[1]);
}
void nearfieldStruct::calcSpheres()
{
//calculate all of the constants necessary to evaluate the scattered field
//estimate the order required to represent the scattered field for each sphere
|
51b6469a
dmayerich
added look-up tables
|
106
107
|
|
3f56f1f9
dmayerich
initial commit
|
108
|
//for each sphere
|
274c3d2b
dmayerich
fixed visual stud...
|
109
|
for(unsigned int i=0; i<sVector.size(); i++)
|
3f56f1f9
dmayerich
initial commit
|
110
111
112
113
114
|
{
//a = sVector[i].a;
//calculate the required order
sVector[i].calcNl(lambda);
|
3f56f1f9
dmayerich
initial commit
|
115
116
117
|
//set the refractive index for the sphere
int imat = sVector[i].iMaterial;
|
3f36b18e
David Mayerich
Adding planewave ...
|
118
|
rts::complex<ptype> n = mVector[imat](lambda);
|
3f56f1f9
dmayerich
initial commit
|
119
120
121
122
123
124
125
|
//calculate the scattering coefficients
sVector[i].calcCoeff(lambda, n);
//save the refractive index
sVector[i].n = n;
|
51b6469a
dmayerich
added look-up tables
|
126
127
128
129
130
131
132
|
//if the LUT is used, calculate Usp(theta, r)
if(lut_us)
{
sVector[i].calcUp(lambda, n, pos, max(U.R[0], U.R[1]));
}
|
3f56f1f9
dmayerich
initial commit
|
133
134
135
136
|
}
}
|
51b6469a
dmayerich
added look-up tables
|
137
138
139
|
void nearfieldStruct::calcUs()
{
|
3f36b18e
David Mayerich
Adding planewave ...
|
140
141
142
143
144
145
146
147
148
149
150
|
if(scalarSim)
{
if(lut_us)
scalarUpLut();
else
scalarUs();
if(lut_us)
std::cout<<"Using LUT for Us"<<std::endl;
std::cout<<"Calculate Us Scalar Sim."<<std::endl;
}
|
51b6469a
dmayerich
added look-up tables
|
151
|
else
|
3f36b18e
David Mayerich
Adding planewave ...
|
152
153
154
|
{
std::cout<<"Calculate Us Vector Sim."<<std::endl;
}
|
51b6469a
dmayerich
added look-up tables
|
155
156
157
158
|
}
void nearfieldStruct::calcUf()
{
|
3f36b18e
David Mayerich
Adding planewave ...
|
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
|
if(scalarSim)
{
if(lut_uf)
scalarUfLut();
else
scalarUf();
//if(lut_uf)
// std::cout<<"Using LUT for Uf"<<std::endl;
//std::cout<<"Calculate Us Scalar Sim."<<std::endl;
}
else
{
std::cout<<"Calculating Uf Vector Sim..."<<std::endl;
if(lut_uf)
vectorUfLut();
else
vectorUf();
}
|
51b6469a
dmayerich
added look-up tables
|
181
182
|
}
|
3f56f1f9
dmayerich
initial commit
|
183
184
|
void nearfieldStruct::Simulate()
{
|
51b6469a
dmayerich
added look-up tables
|
185
186
187
188
|
//initialize timings
t_Uf = 0;
t_Us = 0;
|
3f56f1f9
dmayerich
initial commit
|
189
190
191
192
|
//compute a set of plane waves for Monte-Carlo simulation
calcWaves();
//the near field has to be simulated no matter what the output rtsPoint is
|
51b6469a
dmayerich
added look-up tables
|
193
|
calcUf();
|
3f56f1f9
dmayerich
initial commit
|
194
|
calcSpheres();
|
51b6469a
dmayerich
added look-up tables
|
195
|
calcUs();
|
3f56f1f9
dmayerich
initial commit
|
196
|
sumUf();
|
51b6469a
dmayerich
added look-up tables
|
197
|
|
318b88d0
dmayerich
consolidate output
|
198
199
|
if(verbose)
U.Mag().toImage("testU.bmp");
|
51b6469a
dmayerich
added look-up tables
|
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
227
228
229
230
231
232
233
|
}
void nearfieldStruct::calcBesselLut(ptype* j, ptype d_min, ptype d_max, int dR)
{
/*Compute the look-up-table for spherical bessel functions used for the incident field
j = (Nl + 1) x aR array of values
aR = resolution of j
*/
//compute the wavenumber
ptype k = 2 * PI / lambda;
unsigned int Nl = m;
//allocate space for the Bessel functions of the first and second kind (and derivatives -- which will be ignored)
int bytes = sizeof(double) * (Nl + 1);
double* cjv_kd = (double*)malloc(bytes);
double* cyv_kd = (double*)malloc(bytes);
double* cjvp_kd = (double*)malloc(bytes);
double* cyvp_kd = (double*)malloc(bytes);
//compute the bessel functions using the CPU-based algorithm
double vm;
//for each sample along r
ptype dr = (d_max - d_min) / (dR - 1);
ptype d;
ptype jv;
for(int id = 0; id < dR; id++)
{
d = id * dr + d_min;
double kd = k*d;
bessjyv_sph(Nl, kd, vm, cjv_kd, cyv_kd, cjvp_kd, cyvp_kd);
//copy the double data to the bsComplex array
|
274c3d2b
dmayerich
fixed visual stud...
|
234
|
for(unsigned int l=0; l<=Nl; l++)
|
51b6469a
dmayerich
added look-up tables
|
235
236
237
238
239
240
241
242
243
244
245
246
247
|
{
jv = cjv_kd[l];
if(isnan(jv) || isinf(jv))
{
if(kd == 0 && l == 0)
jv = 1;
else
jv = 0;
}
j[id * (Nl+1) + l] = jv;
}
}
|
318b88d0
dmayerich
consolidate output
|
248
|
if(verbose)
|
51b6469a
dmayerich
added look-up tables
|
249
|
{
|
318b88d0
dmayerich
consolidate output
|
250
251
|
ofstream outfile("uf_besselout.txt");
for(int ir = 0; ir < dR; ir++)
|
51b6469a
dmayerich
added look-up tables
|
252
|
{
|
318b88d0
dmayerich
consolidate output
|
253
|
outfile<<ir*dr + d_min<<endl;
|
274c3d2b
dmayerich
fixed visual stud...
|
254
|
for(unsigned int l = 0; l<=Nl; l++)
|
318b88d0
dmayerich
consolidate output
|
255
256
257
258
|
{
outfile<<j[ir * (Nl+1) + l]<<" -- ";
}
outfile<<endl;
|
51b6469a
dmayerich
added look-up tables
|
259
|
}
|
318b88d0
dmayerich
consolidate output
|
260
|
outfile.close();
|
51b6469a
dmayerich
added look-up tables
|
261
|
}
|
51b6469a
dmayerich
added look-up tables
|
262
|
|
3f56f1f9
dmayerich
initial commit
|
263
|
}
|