planewave.h
4.45 KB
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
#ifndef RTS_PLANEWAVE
#define RTS_PLANEWAVE
#include <string>
#include <sstream>
#include "../math/vector.h"
#include "../math/quaternion.h"
#include "../math/constants.h"
/*Basic conversions used here (assuming a vacuum)
lambda =
*/
namespace rts{
template<typename P>
class planewave{
public:
vec<P> k; //k = tau / lambda
vec<P> E0; //amplitude
///constructor: create a plane wave propagating along z, polarized along x
/*planewave(P lambda = (P)1)
{
k = rts::vec<P>(0, 0, 1) * (TAU/lambda);
E0 = rts::vec<P>(1, 0, 0);
}*/
///constructor: create a plane wave propagating along _k, polarized along _E0, at frequency _omega
planewave(vec<P> _k = rts::vec<P>(0, 0, TAU), vec<P> _E0 = rts::vec<P>(1, 0, 0))
{
k = _k;
vec<P> k_hat = _k.norm();
vec<P> s = (k.cross(_E0)).norm(); //compute an orthogonal side vector
vec<P> E0_hat = (s.cross(k)).norm(); //compute a normalized E0 direction vector
E0 = E0_hat * E0_hat.dot(_E0); //compute the projection of _E0 onto E0_hat
}
///multiplication operator: scale E0
planewave<P> & operator* (const P & rhs)
{
E0 = E0 * rhs;
return *this;
}
P lambda() const
{
return TAU / k.len();
}
P kmag() const
{
return k.len();
}
planewave<P> refract(rts::vec<P> kn) const
{
vec<P> kn_hat = kn.norm(); //normalize the new k
vec<P> k_hat = k.norm(); //normalize the current k
vec<P> r = k_hat.cross(kn_hat); //compute the rotation vector
P theta = asin(r.len()); //compute the angle of the rotation about r
planewave<P> new_p; //create a new plane wave
//deal with a zero vector (both k and kn point in the same direction)
if(theta == (P)0)
{
new_p = *this;
return new_p;
}
//create a quaternion to capture the rotation
quaternion<P> q;
q.CreateRotation(theta, r.norm());
//apply the rotation to E0
vec<P> E0n = q.toMatrix3() * E0;
new_p.k = kn_hat * k.len();
new_p.E0 = E0n;
return new_p;
/*vec<P> kn_hat = kn.norm(); //normalize new_k
vec<P> k_hat = k.norm();
//compute the side vector (around which we will be rotating)
vec<P> E0_hat = E0.norm();
rts::vec<P> s = k_hat.cross(E0_hat);
//compute the projection of k' onto the k-E plane
rts::vec<P> s_prime = s * (kn_hat.dot(s));
//compute the angle between
vec<P> kp_hat = (kn_hat - s_prime).norm();
P theta = acos(k_hat.dot( kp_hat ));
if(kn_hat.dot(E0_hat) < 0)
theta = -theta;
//rotate E0 around s by theta
quaternion<P> q;
q.CreateRotation(theta, s);
rts::vec<P> E0_prime = q.toMatrix3() * E0;
//create the refracted plane wave
planewave<P> new_p;
new_p.E0 = E0_prime;
new_p.k = kn_hat * k.len();
return new_p;*/
/*vec<P> kn_hat = kn.norm(); //normalize kn
vec<P> k_hat = k.norm();
vec<P> E0_hat = E0.norm();
vec<P> B = k_hat.cross(E0_hat);
planewave<P> newp;
newp.k = kn_hat * k.len();
newp.E0 = B.cross(kn_hat).norm();
std::cout<<newp.k.norm().dot(newp.E0.norm())<<std::endl;
return newp;*/
/*
//compute the side vector (around which we will be rotating)
rts::vec<P> s_hat = k_hat.cross(E0_hat);
//std::cout<<s.len()<<std::endl;
//project kn_hat into the k-E0 plane
rts::vec<P> sp = s_hat * (kn_hat.dot(s_hat)); //project k_new onto s
rts::vec<P> kp = (kn_hat - sp); //correct k_new so it lies on the E0-k plane
rts::vec<P> kp_hat = kp.norm();
//compute the angle and direction between k_prime and k
P theta = acos(k_hat.dot(kp_hat));
if(kp_hat.dot(E0_hat) < 0)
theta = -theta;
//rotate E0 around s by theta
quaternion<P> q;
q.CreateRotation(theta, s_hat);
rts::vec<P> E0n = q.toMatrix3() * E0;
rts::vec<P> E0n_hat = E0n.norm();
//std::cout<<s_hat.dot(kp_hat)<<" "<<s_hat.dot(E0n_hat)<<" "<<s_hat.dot(E0_hat)<<" "<<s_hat.dot(k_hat)<<" "<<
// E0_hat.dot(k_hat)<<" "<<k_hat.dot(kp_hat)<<" "<<E0_hat.dot(E0n_hat)<<" "<<E0n_hat.dot(kp_hat)<<std::endl;
//create the refracted plane wave
//std::cout<<"cos: "<<cos(theta)<<" k*kp: "<<k_hat.dot(kp_hat)<<" E0*E0p: "<<E0_hat.dot(E0n_hat)<<" E0p*kp: "<<E0n_hat.dot(kp_hat)<<std::endl;
//std::cout<<"kp*s: "<<kp.dot(s_hat)<<" E0n*s: "<<E0n.dot(s_hat)<<" |E0n|: "<<E0n.len()<<" E0n*kp: "<<E0n.dot(kp_hat)<<" E0n*kn: "<<E0n.dot(kn_hat)<<std::endl;
planewave<P> new_p(kn_hat * k.len(), E0n); //create the plane wave
std::cout<<"|E0n|: "<<new_p.E0.len()<<" E0n*kn: "<<(new_p.E0.norm()).dot(new_p.k.norm())<<std::endl;
return new_p;*/
}
std::string str()
{
std::stringstream ss;
ss<<E0<<" e^i ( "<<k<<" . r )";
return ss.str();
}
};
}
#endif