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
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
|
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)
{
pos = rts::rtsQuad<ptype, 3>(pMin, pMax, normal);
}
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
109
110
111
112
113
114
|
//for each sphere
for(int i=0; i<sVector.size(); i++)
{
//a = sVector[i].a;
//calculate the required order
sVector[i].calcNl(lambda);
|
3f56f1f9
dmayerich
initial commit
|
115
116
117
118
|
//set the refractive index for the sphere
int imat = sVector[i].iMaterial;
rts::rtsComplex<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
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
|
void nearfieldStruct::calcUs()
{
if(lut_us)
scalarUpLut();
else
scalarUs();
}
void nearfieldStruct::calcUf()
{
if(lut_uf)
scalarUfLut();
else
scalarUf();
}
|
3f56f1f9
dmayerich
initial commit
|
155
156
|
void nearfieldStruct::Simulate()
{
|
51b6469a
dmayerich
added look-up tables
|
157
158
159
160
|
//initialize timings
t_Uf = 0;
t_Us = 0;
|
3f56f1f9
dmayerich
initial commit
|
161
162
163
164
|
//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
|
165
|
calcUf();
|
3f56f1f9
dmayerich
initial commit
|
166
|
calcSpheres();
|
51b6469a
dmayerich
added look-up tables
|
167
|
calcUs();
|
3f56f1f9
dmayerich
initial commit
|
168
|
sumUf();
|
51b6469a
dmayerich
added look-up tables
|
169
|
|
318b88d0
dmayerich
consolidate output
|
170
171
|
if(verbose)
U.Mag().toImage("testU.bmp");
|
51b6469a
dmayerich
added look-up tables
|
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
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
|
}
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
for(int l=0; l<=Nl; l++)
{
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
|
220
|
if(verbose)
|
51b6469a
dmayerich
added look-up tables
|
221
|
{
|
318b88d0
dmayerich
consolidate output
|
222
223
|
ofstream outfile("uf_besselout.txt");
for(int ir = 0; ir < dR; ir++)
|
51b6469a
dmayerich
added look-up tables
|
224
|
{
|
318b88d0
dmayerich
consolidate output
|
225
226
227
228
229
230
|
outfile<<ir*dr + d_min<<endl;
for(int l = 0; l<=Nl; l++)
{
outfile<<j[ir * (Nl+1) + l]<<" -- ";
}
outfile<<endl;
|
51b6469a
dmayerich
added look-up tables
|
231
|
}
|
318b88d0
dmayerich
consolidate output
|
232
|
outfile.close();
|
51b6469a
dmayerich
added look-up tables
|
233
|
}
|
51b6469a
dmayerich
added look-up tables
|
234
|
|
3f56f1f9
dmayerich
initial commit
|
235
|
}
|