3f56f1f9
 dmayerich
initial commit
|
1
2
3
4
5
6
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
141
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
194
195
196
197
198
199
|
#include "nearfield.h"
#include "rts/sbessel.h"
#include "rts/legendre.h"
#include <stdlib.h>
#include "rts/cuda_handle_error.h"
#include "rts/cuda_timer.h"
__global__ void gpuScalarUfp(bsComplex* Uf, bsVector k, ptype kmag, bsPoint f, ptype A, bsRect ABCD, int uR, int vR)
{
/*Compute the scalar focused field using Debye focusing
k = direction of focused light, where |k| = 2*pi/lambda
P = rect struct describing the field slice
rX, rY = resolution of the field slice
cNAin = inner NA of the condenser
cNAout = outer NA of the condenser
*/
//get the current coordinate in the plane slice
int iu = blockIdx.x * blockDim.x + threadIdx.x;
int iv = blockIdx.y * blockDim.y + threadIdx.y;
//make sure that the thread indices are in-bounds
if(iu >= uR || iv >= vR) return;
//compute the index (easier access to the scalar field array)
int i = iv*uR + iu;
//compute the parameters for u and v
ptype u = (ptype)iu / uR;
ptype v = (ptype)iv / vR;
//get the rtsPoint in world space and then the r vector
bsPoint p = ABCD(u, v);
bsVector r = p - f;
//ptype d = r.len();
ptype k_dot_r = kmag * k.dot(r);
bsComplex d(0, k_dot_r);
Uf[i] = exp(d) * A;
}
__global__ void gpuScalarUf(bsComplex* Uf, bsVector k, ptype kmag, bsPoint f, ptype A, bsRect ABCD, int uR, int vR, ptype cosAlpha, ptype cosBeta, int nl, ptype j_conv = 1.4)
{
/*Compute the scalar focused field using Debye focusing
k = direction of focused light, where |k| = 2*pi/lambda
P = rect struct describing the field slice
rX, rY = resolution of the field slice
cNAin = inner NA of the condenser
cNAout = outer NA of the condenser
*/
//get the current coordinate in the plane slice
int iu = blockIdx.x * blockDim.x + threadIdx.x;
int iv = blockIdx.y * blockDim.y + threadIdx.y;
//make sure that the thread indices are in-bounds
if(iu >= uR || iv >= vR) return;
//compute the index (easier access to the scalar field array)
int i = iv*uR + iu;
//compute the parameters for u and v
ptype u = (ptype)iu / (uR);
ptype v = (ptype)iv / (vR);
//get the rtsPoint in world space and then the r vector
bsPoint p = ABCD(u, v);
bsVector r = p - f;
ptype d = r.len();
if(d < EPSILON_FLOAT)
{
Uf[i] = A * 2 * PI * (cosAlpha - cosBeta);
return;
}
//get info for the light direction and frequency
//k = k.norm();
r = r.norm();
//compute the imaginary factor i^l
bsComplex im = bsComplex(0, 1);
bsComplex il = bsComplex(1, 0);
//Bessel and Legendre functions are computed dynamically to save memory
//initialize the Bessel and Legendre functions
ptype j[2];
ptype kd = kmag * d;
rts::init_sbesselj<ptype>(kd, j);
ptype P[2];
//get the angle between k and r (light direction and position vector)
ptype cosTheta;
cosTheta = k.dot(r);
//deal with the degenerate case where r == 0
//if(isnan(cosTheta))
// cosTheta = 0;
rts::init_legendre<ptype>(cosTheta, P[0], P[1]);
//initialize legendre functions for the cassegrain angles
ptype Palpha[3];
//ptype cosAlpha = cos(asin(cNAin));
rts::init_legendre<ptype>(cosAlpha, Palpha[0], Palpha[1]);
Palpha[2] = 1;
ptype Pbeta[3];
//ptype cosBeta = cos(asin(cNAout));
rts::init_legendre<ptype>(cosBeta, Pbeta[0], Pbeta[1]);
Pbeta[2] = 1;
//for each order l
bsComplex sumUf(0.0, 0.0);
ptype jl = 0.0;
ptype Pl;
for(int l = 0; l<=nl; l++)
{
if(l==0)
{
jl = j[0];
Pl = P[0];
}
else if(l==1)
{
jl = j[1];
Pl = P[1];
//adjust the cassegrain Legendre function
Palpha[2] = Palpha[0];
rts::shift_legendre<ptype>(l+1, cosAlpha, Palpha[0], Palpha[1]);
Pbeta[2] = Pbeta[0];
rts::shift_legendre<ptype>(l+1, cosBeta, Pbeta[0], Pbeta[1]);
}
else
{
rts::shift_sbesselj<ptype>(l, kd, j);//, j_conv);
rts::shift_legendre<ptype>(l, cosTheta, P[0], P[1]);
jl = j[1];
Pl = P[1];
//adjust the cassegrain outer Legendre function
Palpha[2] = Palpha[0];
rts::shift_legendre<ptype>(l+1, cosAlpha, Palpha[0], Palpha[1]);
Pbeta[2] = Pbeta[0];
rts::shift_legendre<ptype>(l+1, cosBeta, Pbeta[0], Pbeta[1]);
}
sumUf += il * jl * Pl * (Palpha[1] - Palpha[2] - Pbeta[1] + Pbeta[2]);
il *= im;
}
Uf[i] = sumUf * 2 * PI * A;
}
void nearfieldStruct::scalarUf()
{
//Compute the incident field via a scalar simulation
//This method uses Debye focusing to approximate the field analytically
//time the calculation of the focused field
//gpuStartTimer();
//set the field slice to a scalar field
//Uf.scalarField = true;
//initialize the GPU arrays
//Uf.init_gpu();
//create one thread for each pixel of the field slice
dim3 dimBlock(SQRT_BLOCK, SQRT_BLOCK);
dim3 dimGrid((Uf.R[0] + SQRT_BLOCK -1)/SQRT_BLOCK, (Uf.R[1] + SQRT_BLOCK - 1)/SQRT_BLOCK);
//if we are computing a plane wave, call the gpuScalarUfp function
if(planeWave)
{
gpuScalarUfp<<<dimGrid, dimBlock>>>(Uf.x_hat, k, 2 * PI / lambda, focus, A, pos, Uf.R[0], Uf.R[1]);
}
//otherwise compute the condenser info and create a focused field
else
{
//pre-compute the cosine of the obscuration and objective angles
ptype cosAlpha = cos(asin(condenser[0]));
ptype cosBeta = cos(asin(condenser[1]));
//compute the scalar Uf field (this will be in the x_hat channel of Uf)
gpuScalarUf<<<dimGrid, dimBlock>>>(Uf.x_hat, k, 2 * PI / lambda, focus, A, pos, Uf.R[0], Uf.R[1], cosAlpha, cosBeta, m);
}
//float t = gpuStopTimer();
//std::cout<<"Scalar Uf Time: "<<t<<"ms"<<std::endl;
//std::cout<<focus<<std::endl;
}
|
