8e4f8364
David Mayerich
started a new opt...
|
1
2
|
#ifndef STIM_PLANEWAVE_H
#define STIM_PLANEWAVE_H
|
a9275be5
David Mayerich
added vector fiel...
|
3
4
5
|
#include <string>
#include <sstream>
|
8e4f8364
David Mayerich
started a new opt...
|
6
|
#include <cmath>
|
a9275be5
David Mayerich
added vector fiel...
|
7
|
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
8
9
10
|
#include "../math/vector.h"
#include "../math/quaternion.h"
#include "../math/constants.h"
|
559d0fcb
David Mayerich
wave interactions...
|
11
|
#include "../math/plane.h"
|
8e4f8364
David Mayerich
started a new opt...
|
12
|
#include "../math/complex.h"
|
a9275be5
David Mayerich
added vector fiel...
|
13
|
|
8a86bd56
David Mayerich
changed rts names...
|
14
|
namespace stim{
|
8e4f8364
David Mayerich
started a new opt...
|
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
|
namespace optics{
/// evaluate the scalar field produced by a plane wave at a point (x, y, z)
/// @param x is the x-coordinate of the point
/// @param y is the y-coordinate of the point
/// @param z is the z-coordinate of the point
/// @param A is the amplitude of the plane wave, specifically the field at (0, 0, 0)
/// @param kx is the k-vector component in the x direction
/// @param ky is the k-vector component in the y direction
/// @param kz is the k-vector component in the z direction
template<typename T>
stim::complex<T> planewave_scalar(T x, T y, T z, stim::complex<T> A, T kx, T ky, T kz){
T d = x * kx + y * ky + z * kz; //calculate the dot product between k and p = (x, y, z) to find the distance p is along the propagation direction
stim::complex<T> di = stim::complex<T>(0, d); //calculate the phase shift that will have to be applied to propagate the wave distance d
return A * exp(di); //multiply the phase term by the amplitude at (0, 0, 0) to propagate the wave to p
}
/// evaluate the scalar field produced by a plane wave at several positions
/// @param field is a pre-allocated block of memory that will store the complex field at all points
/// @param N is the number of field values to be evaluated
/// @param x is a set of x coordinates defining positions within the field (NULL implies that all values are zero)
/// @param y is a set of y coordinates defining positions within the field (NULL implies that all values are zero)
/// @param z is a set of z coordinates defining positions within the field (NULL implies that all values are zero)
/// @param A is the amplitude of the plane wave, specifically the field at (0, 0, 0)
/// @param kx is the k-vector component in the x direction
/// @param ky is the k-vector component in the y direction
/// @param kz is the k-vector component in the z direction
template<typename T>
void cpu_planewave_scalar(stim::complex<T>* field, size_t N, T* x, T* y = NULL, T* z = NULL, stim::complex<T> A = 1.0, T kx = 0.0, T ky = 0.0, T kz = 0.0){
T px, py, pz;
for(size_t i = 0; i < N; i++){ // for each element in the array
(x == NULL) ? px = 0 : px = x[i]; // test for NULL values
(y == NULL) ? py = 0 : py = y[i];
(z == NULL) ? pz = 0 : pz = z[i];
field[i] = planewave_scalar(px, py, pz, A, kx, ky, kz); // call the single-value plane wave function
}
}
|
a9275be5
David Mayerich
added vector fiel...
|
55
|
|
559d0fcb
David Mayerich
wave interactions...
|
56
|
template<typename T>
|
a9275be5
David Mayerich
added vector fiel...
|
57
58
|
class planewave{
|
559d0fcb
David Mayerich
wave interactions...
|
59
60
|
protected:
|
8e4f8364
David Mayerich
started a new opt...
|
61
62
|
stim::vec<T> k; //k-vector, pointed in propagation direction with magnitude |k| = tau / lambda = 2pi / lambda
stim::vec< stim::complex<T> > E0; //amplitude (for a scalar plane wave, only E0[0] is used)
|
559d0fcb
David Mayerich
wave interactions...
|
63
|
|
8e4f8364
David Mayerich
started a new opt...
|
64
65
|
/// Bend a plane wave via refraction, given that the new propagation direction is known
CUDA_CALLABLE planewave<T> bend(stim::vec<T> kn) const{
|
559d0fcb
David Mayerich
wave interactions...
|
66
|
|
8e4f8364
David Mayerich
started a new opt...
|
67
68
|
stim::vec<T> kn_hat = kn.norm(); //normalize the new k
stim::vec<T> k_hat = k.norm(); //normalize the current k
|
559d0fcb
David Mayerich
wave interactions...
|
69
|
|
8e4f8364
David Mayerich
started a new opt...
|
70
|
planewave<T> new_p; //create a new plane wave
|
559d0fcb
David Mayerich
wave interactions...
|
71
|
|
8e4f8364
David Mayerich
started a new opt...
|
72
|
T k_dot_kn = k_hat.dot(kn_hat); //if kn is equal to k or -k, handle the degenerate case
|
559d0fcb
David Mayerich
wave interactions...
|
73
74
|
//if k . n < 0, then the bend is a reflection
|
8e4f8364
David Mayerich
started a new opt...
|
75
|
if(k_dot_kn < 0) k_hat = -k_hat; //flip k_hat
|
559d0fcb
David Mayerich
wave interactions...
|
76
|
|
559d0fcb
David Mayerich
wave interactions...
|
77
78
79
80
81
82
83
84
85
86
87
|
if(k_dot_kn == -1){
new_p.k = -k;
new_p.E0 = E0;
return new_p;
}
else if(k_dot_kn == 1){
new_p.k = k;
new_p.E0 = E0;
return new_p;
}
|
8e4f8364
David Mayerich
started a new opt...
|
88
89
90
91
92
|
vec<T> r = k_hat.cross(kn_hat); //compute the rotation vector
T theta = asin(r.len()); //compute the angle of the rotation about r
quaternion<T> q; //create a quaternion to capture the rotation
q.CreateRotation(theta, r.norm());
vec< complex<T> > E0n = q.toMatrix3() * E0; //apply the rotation to E0
|
559d0fcb
David Mayerich
wave interactions...
|
93
94
95
96
97
98
|
new_p.k = kn_hat * kmag();
new_p.E0 = E0n;
return new_p;
}
|
a9275be5
David Mayerich
added vector fiel...
|
99
|
public:
|
a9275be5
David Mayerich
added vector fiel...
|
100
|
|
8e4f8364
David Mayerich
started a new opt...
|
101
102
103
|
///constructor: create a plane wave propagating along k
CUDA_CALLABLE planewave(vec<T> kvec = stim::vec<T>(0, 0, stim::TAU),
vec< complex<T> > E = stim::vec<T>(1, 0, 0))
|
a9275be5
David Mayerich
added vector fiel...
|
104
|
{
|
ecfd14df
David Mayerich
implemented D fie...
|
105
106
|
//phi = phase;
|
559d0fcb
David Mayerich
wave interactions...
|
107
|
k = kvec;
|
ecfd14df
David Mayerich
implemented D fie...
|
108
|
vec< complex<T> > k_hat = k.norm();
|
559d0fcb
David Mayerich
wave interactions...
|
109
110
111
112
|
if(E.len() == 0) //if the plane wave has an amplitude of 0
E0 = vec<T>(0); //just return it
else{
|
ecfd14df
David Mayerich
implemented D fie...
|
113
114
|
vec< complex<T> > s = (k_hat.cross(E)).norm(); //compute an orthogonal side vector
vec< complex<T> > E_hat = (s.cross(k)).norm(); //compute a normalized E0 direction vector
|
8e4f8364
David Mayerich
started a new opt...
|
115
|
E0 = E_hat;// * E_hat.dot(E); //compute the projection of _E0 onto E0_hat
|
559d0fcb
David Mayerich
wave interactions...
|
116
|
}
|
ecfd14df
David Mayerich
implemented D fie...
|
117
118
|
E0 = E0 * exp( complex<T>(0, phase) );
|
a9275be5
David Mayerich
added vector fiel...
|
119
120
121
|
}
///multiplication operator: scale E0
|
8e4f8364
David Mayerich
started a new opt...
|
122
|
CUDA_CALLABLE planewave<T> & operator* (const T & rhs){
|
a9275be5
David Mayerich
added vector fiel...
|
123
124
125
126
|
E0 = E0 * rhs;
return *this;
}
|
8e4f8364
David Mayerich
started a new opt...
|
127
128
|
CUDA_CALLABLE T lambda() const{
return stim::TAU / k.len();
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
129
130
|
}
|
8e4f8364
David Mayerich
started a new opt...
|
131
|
CUDA_CALLABLE T kmag() const{
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
132
133
134
|
return k.len();
}
|
ecfd14df
David Mayerich
implemented D fie...
|
135
|
CUDA_CALLABLE vec< complex<T> > E(){
|
559d0fcb
David Mayerich
wave interactions...
|
136
137
138
139
140
141
142
|
return E0;
}
CUDA_CALLABLE vec<T> kvec(){
return k;
}
|
8e4f8364
David Mayerich
started a new opt...
|
143
144
145
|
/// calculate the value of the field produced by the plane wave given a three-dimensional position
CUDA_CALLABLE vec< complex<T> > pos(T x, T y, T z){
return pos( stim::vec<T>(x, y, z) );
|
559d0fcb
David Mayerich
wave interactions...
|
146
147
|
}
|
559d0fcb
David Mayerich
wave interactions...
|
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
|
CUDA_CALLABLE vec< complex<T> > pos(vec<T> p = vec<T>(0, 0, 0)){
vec< complex<T> > result;
T kdp = k.dot(p);
complex<T> x = complex<T>(0, kdp);
complex<T> expx = exp(x);
result[0] = E0[0] * expx;
result[1] = E0[1] * expx;
result[2] = E0[2] * expx;
return result;
}
//scales k based on a transition from material ni to material nt
CUDA_CALLABLE planewave<T> n(T ni, T nt){
return planewave<T>(k * (nt / ni), E0);
}
|
a9275be5
David Mayerich
added vector fiel...
|
166
|
|
8e4f8364
David Mayerich
started a new opt...
|
167
|
CUDA_CALLABLE planewave<T> refract(stim::vec<T> kn) const{
|
559d0fcb
David Mayerich
wave interactions...
|
168
169
|
return bend(kn);
}
|
d609550e
David Mayerich
fixed bug in plan...
|
170
|
|
8e4f8364
David Mayerich
started a new opt...
|
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
|
/// Calculate the result of a plane wave hitting an interface between two refractive indices
/// @param P is a plane representing the position and orientation of the surface
/// @param n0 is the refractive index outside of the surface (in the direction of the normal)
/// @param n1 is the refractive index inside the surface (in the direction away from the normal)
/// @param r is the reflected component of the plane wave
/// @param t is the transmitted component of the plane wave
void scatter(stim::plane<T> P, T n0, T n1, planewave<T> &r, planewave<T> &t){
scatter(P, n1/n0, r, t);
}
/// Calculate the scattering result when nr = n1/n0
/// @param P is a plane representing the position and orientation of the surface
/// @param r is the ration n1/n0
/// @param n1 is the refractive index inside the surface (in the direction away from the normal)
/// @param r is the reflected component of the plane wave
/// @param t is the transmitted component of the plane wave
void scatter(stim::plane<T> P, T nr, planewave<T> &r, planewave<T> &t){
|
d609550e
David Mayerich
fixed bug in plan...
|
190
|
|
559d0fcb
David Mayerich
wave interactions...
|
191
|
int facing = P.face(k); //determine which direction the plane wave is coming in
|
d609550e
David Mayerich
fixed bug in plan...
|
192
|
|
559d0fcb
David Mayerich
wave interactions...
|
193
194
195
196
|
if(facing == -1){ //if the wave hits the back of the plane, invert the plane and nr
P = P.flip(); //flip the plane
nr = 1/nr; //invert the refractive index (now nr = n0/n1)
}
|
d609550e
David Mayerich
fixed bug in plan...
|
197
|
|
559d0fcb
David Mayerich
wave interactions...
|
198
199
200
201
202
203
204
205
206
|
//use Snell's Law to calculate the transmitted angle
T cos_theta_i = k.norm().dot(-P.norm()); //compute the cosine of theta_i
T theta_i = acos(cos_theta_i); //compute theta_i
T sin_theta_t = (1/nr) * sin(theta_i); //compute the sine of theta_t using Snell's law
T theta_t = asin(sin_theta_t); //compute the cosine of theta_t
bool tir = false; //flag for total internal reflection
if(theta_t != theta_t){
tir = true;
|
8e4f8364
David Mayerich
started a new opt...
|
207
|
theta_t = stim::PI / (T)2;
|
d609550e
David Mayerich
fixed bug in plan...
|
208
209
|
}
|
559d0fcb
David Mayerich
wave interactions...
|
210
211
212
213
214
215
|
//handle the degenerate case where theta_i is 0 (the plane wave hits head-on)
if(theta_i == 0){
T rp = (1 - nr) / (1 + nr); //compute the Fresnel coefficients
T tp = 2 / (1 + nr);
vec<T> kr = -k;
vec<T> kt = k * nr; //set the k vectors for theta_i = 0
|
ecfd14df
David Mayerich
implemented D fie...
|
216
217
|
vec< complex<T> > Er = E0 * rp; //compute the E vectors
vec< complex<T> > Et = E0 * tp;
|
559d0fcb
David Mayerich
wave interactions...
|
218
219
|
T phase_t = P.p().dot(k - kt); //compute the phase offset
T phase_r = P.p().dot(k - kr);
|
559d0fcb
David Mayerich
wave interactions...
|
220
221
222
223
|
//create the plane waves
r = planewave<T>(kr, Er, phase_r);
t = planewave<T>(kt, Et, phase_t);
|
559d0fcb
David Mayerich
wave interactions...
|
224
225
|
return;
}
|
d609550e
David Mayerich
fixed bug in plan...
|
226
|
|
d609550e
David Mayerich
fixed bug in plan...
|
227
|
|
559d0fcb
David Mayerich
wave interactions...
|
228
229
230
231
232
233
234
235
236
237
238
239
|
//compute the Fresnel coefficients
T rp, rs, tp, ts;
rp = tan(theta_t - theta_i) / tan(theta_t + theta_i);
rs = sin(theta_t - theta_i) / sin(theta_t + theta_i);
if(tir){
tp = ts = 0;
}
else{
tp = ( 2 * sin(theta_t) * cos(theta_i) ) / ( sin(theta_t + theta_i) * cos(theta_t - theta_i) );
ts = ( 2 * sin(theta_t) * cos(theta_i) ) / sin(theta_t + theta_i);
}
|
d609550e
David Mayerich
fixed bug in plan...
|
240
|
|
559d0fcb
David Mayerich
wave interactions...
|
241
242
243
244
245
246
247
248
249
250
251
|
//compute the coordinate space for the plane of incidence
vec<T> z_hat = -P.norm();
vec<T> y_hat = P.parallel(k).norm();
vec<T> x_hat = y_hat.cross(z_hat).norm();
//compute the k vectors for r and t
vec<T> kr, kt;
kr = ( y_hat * sin(theta_i) - z_hat * cos(theta_i) ) * kmag();
kt = ( y_hat * sin(theta_t) + z_hat * cos(theta_t) ) * kmag() * nr;
//compute the magnitude of the p- and s-polarized components of the incident E vector
|
ecfd14df
David Mayerich
implemented D fie...
|
252
|
complex<T> Ei_s = E0.dot(x_hat);
|
ecfd14df
David Mayerich
implemented D fie...
|
253
254
255
|
int sgn = E0.dot(y_hat).sgn();
vec< complex<T> > cx_hat = x_hat;
complex<T> Ei_p = ( E0 - cx_hat * Ei_s ).len() * sgn;
|
559d0fcb
David Mayerich
wave interactions...
|
256
|
//compute the magnitude of the p- and s-polarized components of the reflected E vector
|
ecfd14df
David Mayerich
implemented D fie...
|
257
258
|
complex<T> Er_s = Ei_s * rs;
complex<T> Er_p = Ei_p * rp;
|
559d0fcb
David Mayerich
wave interactions...
|
259
|
//compute the magnitude of the p- and s-polarized components of the transmitted E vector
|
ecfd14df
David Mayerich
implemented D fie...
|
260
261
|
complex<T> Et_s = Ei_s * ts;
complex<T> Et_p = Ei_p * tp;
|
559d0fcb
David Mayerich
wave interactions...
|
262
|
|
559d0fcb
David Mayerich
wave interactions...
|
263
|
//compute the reflected E vector
|
ecfd14df
David Mayerich
implemented D fie...
|
264
|
vec< complex<T> > Er = vec< complex<T> >(y_hat * cos(theta_i) + z_hat * sin(theta_i)) * Er_p + cx_hat * Er_s;
|
559d0fcb
David Mayerich
wave interactions...
|
265
|
//compute the transmitted E vector
|
ecfd14df
David Mayerich
implemented D fie...
|
266
|
vec< complex<T> > Et = vec< complex<T> >(y_hat * cos(theta_t) - z_hat * sin(theta_t)) * Et_p + cx_hat * Et_s;
|
559d0fcb
David Mayerich
wave interactions...
|
267
268
269
270
|
T phase_t = P.p().dot(k - kt);
T phase_r = P.p().dot(k - kr);
|
559d0fcb
David Mayerich
wave interactions...
|
271
272
|
//create the plane waves
r.k = kr;
|
ecfd14df
David Mayerich
implemented D fie...
|
273
|
r.E0 = Er * exp( complex<T>(0, phase_r) );
|
559d0fcb
David Mayerich
wave interactions...
|
274
275
|
t.k = kt;
|
ecfd14df
David Mayerich
implemented D fie...
|
276
|
t.E0 = Et * exp( complex<T>(0, phase_t) );
|
a9275be5
David Mayerich
added vector fiel...
|
277
278
279
280
281
|
}
std::string str()
{
std::stringstream ss;
|
559d0fcb
David Mayerich
wave interactions...
|
282
283
|
ss<<"Plane Wave:"<<std::endl;
ss<<" "<<E0<<" e^i ( "<<k<<" . r )";
|
a9275be5
David Mayerich
added vector fiel...
|
284
285
|
return ss.str();
}
|
8e4f8364
David Mayerich
started a new opt...
|
286
287
288
|
}; //end planewave class
} //end namespace optics
} //end namespace stim
|
a9275be5
David Mayerich
added vector fiel...
|
289
|
|
559d0fcb
David Mayerich
wave interactions...
|
290
|
template <typename T>
|
8e4f8364
David Mayerich
started a new opt...
|
291
|
std::ostream& operator<<(std::ostream& os, stim::optics::planewave<T> p)
|
559d0fcb
David Mayerich
wave interactions...
|
292
293
294
295
296
|
{
os<<p.str();
return os;
}
|
8e4f8364
David Mayerich
started a new opt...
|
297
|
#endif
|