a9275be5
David Mayerich
added vector fiel...
|
1
2
3
4
5
6
|
#ifndef RTS_PLANEWAVE
#define RTS_PLANEWAVE
#include <string>
#include <sstream>
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
7
8
9
|
#include "../math/vector.h"
#include "../math/quaternion.h"
#include "../math/constants.h"
|
559d0fcb
David Mayerich
wave interactions...
|
10
|
#include "../math/plane.h"
|
8a86bd56
David Mayerich
changed rts names...
|
11
|
#include "../cuda/callable.h"
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
12
13
14
15
|
/*Basic conversions used here (assuming a vacuum)
lambda =
*/
|
a9275be5
David Mayerich
added vector fiel...
|
16
|
|
8a86bd56
David Mayerich
changed rts names...
|
17
|
namespace stim{
|
a9275be5
David Mayerich
added vector fiel...
|
18
|
|
559d0fcb
David Mayerich
wave interactions...
|
19
|
template<typename T>
|
a9275be5
David Mayerich
added vector fiel...
|
20
21
|
class planewave{
|
559d0fcb
David Mayerich
wave interactions...
|
22
23
24
|
protected:
vec<T> k; //k = tau / lambda
|
ecfd14df
David Mayerich
implemented D fie...
|
25
26
|
vec< complex<T> > E0; //amplitude
//T phi;
|
559d0fcb
David Mayerich
wave interactions...
|
27
|
|
ecfd14df
David Mayerich
implemented D fie...
|
28
|
CUDA_CALLABLE planewave<T> bend(rts::vec<T> kn) const{
|
559d0fcb
David Mayerich
wave interactions...
|
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
|
vec<T> kn_hat = kn.norm(); //normalize the new k
vec<T> k_hat = k.norm(); //normalize the current k
//std::cout<<"PLANE WAVE BENDING------------------"<<std::endl;
//std::cout<<"kn_hat: "<<kn_hat<<" k_hat: "<<k_hat<<std::endl;
planewave<T> new_p; //create a new plane wave
//if kn is equal to k or -k, handle the degenerate case
T k_dot_kn = k_hat.dot(kn_hat);
//if k . n < 0, then the bend is a reflection
//flip k_hat
if(k_dot_kn < 0) k_hat = -k_hat;
//std::cout<<"k dot kn: "<<k_dot_kn<<std::endl;
//std::cout<<"k_dot_kn: "<<k_dot_kn<<std::endl;
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;
}
vec<T> r = k_hat.cross(kn_hat); //compute the rotation vector
//std::cout<<"r: "<<r<<std::endl;
T theta = asin(r.len()); //compute the angle of the rotation about r
//deal with a zero vector (both k and kn point in the same direction)
//if(theta == (T)0)
//{
// new_p = *this;
// return new_p;
//}
//create a quaternion to capture the rotation
quaternion<T> q;
q.CreateRotation(theta, r.norm());
//apply the rotation to E0
|
ecfd14df
David Mayerich
implemented D fie...
|
79
|
vec< complex<T> > E0n = q.toMatrix3() * E0;
|
559d0fcb
David Mayerich
wave interactions...
|
80
81
82
83
84
85
86
|
new_p.k = kn_hat * kmag();
new_p.E0 = E0n;
return new_p;
}
|
a9275be5
David Mayerich
added vector fiel...
|
87
|
public:
|
a9275be5
David Mayerich
added vector fiel...
|
88
89
|
|
a9275be5
David Mayerich
added vector fiel...
|
90
|
///constructor: create a plane wave propagating along z, polarized along x
|
559d0fcb
David Mayerich
wave interactions...
|
91
|
/*planewave(T lambda = (T)1)
|
a9275be5
David Mayerich
added vector fiel...
|
92
|
{
|
559d0fcb
David Mayerich
wave interactions...
|
93
94
|
k = rts::vec<T>(0, 0, 1) * (TAU/lambda);
E0 = rts::vec<T>(1, 0, 0);
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
95
|
}*/
|
559d0fcb
David Mayerich
wave interactions...
|
96
97
|
///constructor: create a plane wave propagating along k, polarized along _E0, at frequency _omega
CUDA_CALLABLE planewave(vec<T> kvec = rts::vec<T>(0, 0, rtsTAU),
|
ecfd14df
David Mayerich
implemented D fie...
|
98
|
vec< complex<T> > E = rts::vec<T>(1, 0, 0), T phase = 0)
|
a9275be5
David Mayerich
added vector fiel...
|
99
|
{
|
ecfd14df
David Mayerich
implemented D fie...
|
100
101
|
//phi = phase;
|
559d0fcb
David Mayerich
wave interactions...
|
102
|
k = kvec;
|
ecfd14df
David Mayerich
implemented D fie...
|
103
|
vec< complex<T> > k_hat = k.norm();
|
559d0fcb
David Mayerich
wave interactions...
|
104
105
106
107
|
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...
|
108
109
|
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
|
559d0fcb
David Mayerich
wave interactions...
|
110
111
|
E0 = E_hat * E_hat.dot(E); //compute the projection of _E0 onto E0_hat
}
|
ecfd14df
David Mayerich
implemented D fie...
|
112
113
|
E0 = E0 * exp( complex<T>(0, phase) );
|
a9275be5
David Mayerich
added vector fiel...
|
114
115
116
|
}
///multiplication operator: scale E0
|
559d0fcb
David Mayerich
wave interactions...
|
117
|
CUDA_CALLABLE planewave<T> & operator* (const T & rhs)
|
a9275be5
David Mayerich
added vector fiel...
|
118
119
120
121
122
123
|
{
E0 = E0 * rhs;
return *this;
}
|
559d0fcb
David Mayerich
wave interactions...
|
124
|
CUDA_CALLABLE T lambda() const
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
125
|
{
|
559d0fcb
David Mayerich
wave interactions...
|
126
|
return rtsTAU / k.len();
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
127
128
|
}
|
559d0fcb
David Mayerich
wave interactions...
|
129
|
CUDA_CALLABLE T kmag() const
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
130
131
132
133
|
{
return k.len();
}
|
ecfd14df
David Mayerich
implemented D fie...
|
134
|
CUDA_CALLABLE vec< complex<T> > E(){
|
559d0fcb
David Mayerich
wave interactions...
|
135
136
137
138
139
140
141
|
return E0;
}
CUDA_CALLABLE vec<T> kvec(){
return k;
}
|
ecfd14df
David Mayerich
implemented D fie...
|
142
|
/*CUDA_CALLABLE T phase(){
|
559d0fcb
David Mayerich
wave interactions...
|
143
144
145
146
147
|
return phi;
}
CUDA_CALLABLE void phase(T p){
phi = p;
|
ecfd14df
David Mayerich
implemented D fie...
|
148
|
}*/
|
559d0fcb
David Mayerich
wave interactions...
|
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
|
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...
|
168
|
|
559d0fcb
David Mayerich
wave interactions...
|
169
|
CUDA_CALLABLE planewave<T> refract(rts::vec<T> kn) const
|
a9275be5
David Mayerich
added vector fiel...
|
170
|
{
|
559d0fcb
David Mayerich
wave interactions...
|
171
172
|
return bend(kn);
}
|
d609550e
David Mayerich
fixed bug in plan...
|
173
|
|
559d0fcb
David Mayerich
wave interactions...
|
174
|
void scatter(rts::plane<T> P, T nr, planewave<T> &r, planewave<T> &t){
|
d609550e
David Mayerich
fixed bug in plan...
|
175
|
|
559d0fcb
David Mayerich
wave interactions...
|
176
|
int facing = P.face(k); //determine which direction the plane wave is coming in
|
d609550e
David Mayerich
fixed bug in plan...
|
177
|
|
559d0fcb
David Mayerich
wave interactions...
|
178
179
180
181
182
183
184
|
//if(facing == 0) //if the wave is tangent to the plane, return an identical wave
// return *this;
//else
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...
|
185
|
|
559d0fcb
David Mayerich
wave interactions...
|
186
187
188
189
190
191
192
193
194
195
|
//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;
theta_t = rtsPI / (T)2;
|
d609550e
David Mayerich
fixed bug in plan...
|
196
197
|
}
|
559d0fcb
David Mayerich
wave interactions...
|
198
199
200
201
202
203
|
//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...
|
204
205
|
vec< complex<T> > Er = E0 * rp; //compute the E vectors
vec< complex<T> > Et = E0 * tp;
|
559d0fcb
David Mayerich
wave interactions...
|
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
|
T phase_t = P.p().dot(k - kt); //compute the phase offset
T phase_r = P.p().dot(k - kr);
//std::cout<<"Degeneracy: Head-On"<<std::endl;
//std::cout<<"rs: "<<rp<<" rp: "<<rp<<" ts: "<<tp<<" tp: "<<tp<<std::endl;
//std::cout<<"phase r: "<<phase_r<<" phase t: "<<phase_t<<std::endl;
//create the plane waves
r = planewave<T>(kr, Er, phase_r);
t = planewave<T>(kt, Et, phase_t);
//std::cout<<"i + r: "<<pos()[0] + r.pos()[0]<<pos()[1] + r.pos()[1]<<pos()[2] + r.pos()[2]<<std::endl;
//std::cout<<"t: "<<t.pos()[0]<<t.pos()[1]<<t.pos()[2]<<std::endl;
//std::cout<<"--------------------------------"<<std::endl;
return;
}
|
d609550e
David Mayerich
fixed bug in plan...
|
221
|
|
d609550e
David Mayerich
fixed bug in plan...
|
222
|
|
559d0fcb
David Mayerich
wave interactions...
|
223
224
225
226
227
228
229
230
231
232
233
234
|
//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...
|
235
|
|
559d0fcb
David Mayerich
wave interactions...
|
236
237
238
239
240
241
242
243
244
245
246
|
//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...
|
247
248
249
250
251
|
complex<T> Ei_s = E0.dot(x_hat);
//int sgn = (0 < E0.dot(y_hat)) - (E0.dot(y_hat) < 0);
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...
|
252
253
|
//T Ei_p = ( E0 - x_hat * Ei_s ).len();
//compute the magnitude of the p- and s-polarized components of the reflected E vector
|
ecfd14df
David Mayerich
implemented D fie...
|
254
255
|
complex<T> Er_s = Ei_s * rs;
complex<T> Er_p = Ei_p * rp;
|
559d0fcb
David Mayerich
wave interactions...
|
256
|
//compute the magnitude of the p- and s-polarized components of the transmitted E vector
|
ecfd14df
David Mayerich
implemented D fie...
|
257
258
|
complex<T> Et_s = Ei_s * ts;
complex<T> Et_p = Ei_p * tp;
|
559d0fcb
David Mayerich
wave interactions...
|
259
260
261
262
263
264
265
266
267
268
|
//std::cout<<"E0: "<<E0<<std::endl;
//std::cout<<"E0 dot y_hat: "<<E0.dot(y_hat)<<std::endl;
//std::cout<<"theta i: "<<theta_i<<" theta t: "<<theta_t<<std::endl;
//std::cout<<"x_hat: "<<x_hat<<" y_hat: "<<y_hat<<" z_hat: "<<z_hat<<std::endl;
//std::cout<<"Ei_s: "<<Ei_s<<" Ei_p: "<<Ei_p<<" Er_s: "<<Er_s<<" Er_p: "<<Er_p<<" Et_s: "<<Et_s<<" Et_p: "<<Et_p<<std::endl;
//std::cout<<"rs: "<<rs<<" rp: "<<rp<<" ts: "<<ts<<" tp: "<<tp<<std::endl;
//compute the reflected E vector
|
ecfd14df
David Mayerich
implemented D fie...
|
269
|
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...
|
270
|
//compute the transmitted E vector
|
ecfd14df
David Mayerich
implemented D fie...
|
271
|
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...
|
272
273
274
275
276
277
278
279
280
281
|
T phase_t = P.p().dot(k - kt);
T phase_r = P.p().dot(k - kr);
//std::cout<<"phase r: "<<phase_r<<" phase t: "<<phase_t<<std::endl;
//std::cout<<"phase: "<<phase<<std::endl;
//create the plane waves
r.k = kr;
|
ecfd14df
David Mayerich
implemented D fie...
|
282
283
|
r.E0 = Er * exp( complex<T>(0, phase_r) );
//r.phi = phase_r;
|
559d0fcb
David Mayerich
wave interactions...
|
284
285
286
287
288
|
//t = bend(kt);
//t.k = t.k * nr;
t.k = kt;
|
ecfd14df
David Mayerich
implemented D fie...
|
289
290
|
t.E0 = Et * exp( complex<T>(0, phase_t) );
//t.phi = phase_t;
|
559d0fcb
David Mayerich
wave interactions...
|
291
292
293
294
295
296
297
|
//std::cout<<"i: "<<str()<<std::endl;
//std::cout<<"r: "<<r.str()<<std::endl;
//std::cout<<"t: "<<t.str()<<std::endl;
//std::cout<<"i + r: "<<pos()[0] + r.pos()[0]<<pos()[1] + r.pos()[1]<<pos()[2] + r.pos()[2]<<std::endl;
//std::cout<<"t: "<<t.pos()[0]<<t.pos()[1]<<t.pos()[2]<<std::endl;
//std::cout<<"--------------------------------"<<std::endl;
|
d609550e
David Mayerich
fixed bug in plan...
|
298
|
|
a9275be5
David Mayerich
added vector fiel...
|
299
300
301
302
303
|
}
std::string str()
{
std::stringstream ss;
|
559d0fcb
David Mayerich
wave interactions...
|
304
305
|
ss<<"Plane Wave:"<<std::endl;
ss<<" "<<E0<<" e^i ( "<<k<<" . r )";
|
a9275be5
David Mayerich
added vector fiel...
|
306
307
308
309
310
|
return ss.str();
}
};
}
|
559d0fcb
David Mayerich
wave interactions...
|
311
312
313
314
315
316
317
|
template <typename T>
std::ostream& operator<<(std::ostream& os, rts::planewave<T> p)
{
os<<p.str();
return os;
}
|
8a86bd56
David Mayerich
changed rts names...
|
318
|
#endif
|