codebase / stimlib

Commit b28f312da451d7b7756437ac60ed7ee79db390c2

Authored by David Mayerich
1 parent 024756d5

added code for evaluating waves and beams

Showing 7 changed files with 183 additions and 354 deletions   Show diff stats
Inline Side-by-side
stim/math/matrix_sq.h
  Show/Hide comments   View file @b28f312
... ... @@ -8,7 +8,7 @@
8 8 #include <stim/math/vec3.h>
9 9 #include <stim/cuda/cudatools/callable.h>
10 10  
11   -namespace stim{
  11 +namespace tira {
12 12  
13 13 template <class T, int N>
14 14 struct matrix_sq
... ... @@ -100,7 +100,7 @@ struct matrix_sq
100 100  
101 101 template<typename Y>
102 102 CUDA_CALLABLE vec3<Y> operator*(vec3<Y> rhs){
103   - vec3<Y> result;
  103 + vec3<Y> result(0, 0, 0);
104 104 for(int r=0; r<3; r++)
105 105 for(int c=0; c<3; c++)
106 106 result[r] += (*this)(r, c) * rhs[c];
... ... @@ -152,7 +152,7 @@ struct matrix_sq
152 152 } //end namespace rts
153 153  
154 154 template <typename T, int N>
155   -std::ostream& operator<<(std::ostream& os, stim::matrix_sq<T, N> M)
  155 +std::ostream& operator<<(std::ostream& os, tira::matrix_sq<T, N> M)
156 156 {
157 157 os<<M.toStr();
158 158 return os;
... ...
stim/math/plane.h
  Show/Hide comments   View file @b28f312
1   -#ifndef STIM_PLANE_H
2   -#define STIM_PLANE_H
  1 +#ifndef TIRA_PLANE_H
  2 +#define TIRA_PLANE_H
3 3  
4 4 #include <iostream>
5 5 #include <stim/math/vec3.h>
... ... @@ -7,32 +7,32 @@
7 7 #include <stim/math/quaternion.h>
8 8  
9 9  
10   -namespace stim
  10 +namespace tira
11 11 {
12 12 template<typename T> class plane;
13 13 }
14 14  
15 15 template<typename T>
16   -CUDA_CALLABLE stim::plane<T> operator-(stim::plane<T> v);
  16 +CUDA_CALLABLE tira::plane<T> operator-(tira::plane<T> v);
17 17  
18   -namespace stim
  18 +namespace tira
19 19 {
20 20  
21 21 template <typename T>
22 22 class plane
23 23 {
24 24 protected:
25   - stim::vec3<T> P;
26   - stim::vec3<T> N;
27   - stim::vec3<T> U;
  25 + vec3<T> P;
  26 + vec3<T> N;
  27 + vec3<T> U;
28 28  
29 29 ///Initializes the plane with standard coordinates.
30 30 ///
31 31 CUDA_CALLABLE void init()
32 32 {
33   - P = stim::vec3<T>(0, 0, 0);
34   - N = stim::vec3<T>(0, 0, 1);
35   - U = stim::vec3<T>(1, 0, 0);
  33 + P = vec3<T>(0, 0, 0);
  34 + N = vec3<T>(0, 0, 1);
  35 + U = vec3<T>(1, 0, 0);
36 36 }
37 37  
38 38 public:
... ... @@ -60,7 +60,7 @@ class plane
60 60 {
61 61 init();
62 62 P = c;
63   - stim::vec3<T> n = (c - a).cross(b - a);
  63 + vec3<T> n = (c - a).cross(b - a);
64 64 try
65 65 {
66 66 if(n.len() != 0)
... ... @@ -165,9 +165,9 @@ class plane
165 165  
166 166 }
167 167  
168   - CUDA_CALLABLE stim::plane<T> operator-()
  168 + CUDA_CALLABLE plane<T> operator-()
169 169 {
170   - stim::plane<T> p = *this;
  170 + plane<T> p = *this;
171 171  
172 172 //negate the normal vector
173 173 p.N = -p.N;
... ...
stim/math/quaternion.h
  Show/Hide comments   View file @b28f312
1   -#ifndef RTS_QUATERNION_H
2   -#define RTS_QUATERNION_H
  1 +#ifndef TIRA_QUATERNION_H
  2 +#define TIRA_QUATERNION_H
3 3  
4 4 #include <stim/math/matrix_sq.h>
5 5 #include <stim/cuda/cudatools/callable.h>
6 6  
7   -namespace stim{
  7 +namespace tira {
8 8  
9 9 template<typename T>
10 10 class quaternion
... ...
stim/math/vec3.h
  Show/Hide comments   View file @b28f312
1   -#ifndef STIM_VEC3_H
2   -#define STIM_VEC3_H
  1 +#ifndef TIRA_VEC3_H
  2 +#define TIRA_VEC3_H
3 3  
4 4  
5 5 #include <stim/cuda/cudatools/callable.h>
... ... @@ -9,7 +9,7 @@
9 9 #include <sstream>
10 10  
11 11  
12   -namespace stim{
  12 +namespace tira{
13 13  
14 14  
15 15 /// A class designed to act as a 3D vector with CUDA compatibility
... ... @@ -313,20 +313,20 @@ std::string str() const{
313 313 return result;
314 314 }
315 315 }; //end class cvec3
316   -} //end namespace stim
  316 +} //end namespace tira
317 317  
318 318  
319 319  
320 320  
321 321 /// Multiply a vector by a constant when the vector is on the right hand side
322 322 template <typename T>
323   -stim::vec3<T> operator*(T lhs, stim::vec3<T> rhs){
  323 +tira::vec3<T> operator*(T lhs, tira::vec3<T> rhs){
324 324 return rhs * lhs;
325 325 }
326 326  
327 327 //stream operator
328 328 template<typename T>
329   -std::ostream& operator<<(std::ostream& os, stim::vec3<T> const& rhs){
  329 +std::ostream& operator<<(std::ostream& os, tira::vec3<T> const& rhs){
330 330 os<<rhs.str();
331 331 return os;
332 332 }
... ...
stim/math/vector.h
  Show/Hide comments   View file @b28f312
1   -#ifndef STIM_VECTOR_H
2   -#define STIM_VECTOR_H
  1 +#ifndef TIRA_VECTOR_H
  2 +#define TIRA_VECTOR_H
3 3  
4 4 #include <iostream>
5 5 #include <cmath>
... ... @@ -11,7 +11,7 @@
11 11 #include <stim/cuda/cudatools/callable.h>
12 12 #include <stim/math/vec3.h>
13 13  
14   -namespace stim
  14 +namespace tira
15 15 {
16 16  
17 17 template <class T>
... ... @@ -338,8 +338,8 @@ struct vec : public std::vector&lt;T&gt;
338 338 }
339 339  
340 340 /// Cast to a vec3
341   - operator stim::vec3<T>(){
342   - stim::vec3<T> r;
  341 + operator tira::vec3<T>(){
  342 + tira::vec3<T> r;
343 343 size_t N = std::min(size(), (size_t)3);
344 344 for(size_t i = 0; i < N; i++)
345 345 r[i] = at(i);
... ... @@ -405,10 +405,10 @@ struct vec : public std::vector&lt;T&gt;
405 405 };
406 406  
407 407  
408   -} //end namespace rts
  408 +} //end namespace tira
409 409  
410 410 template <typename T>
411   -std::ostream& operator<<(std::ostream& os, stim::vec<T> v)
  411 +std::ostream& operator<<(std::ostream& os, tira::vec<T> v)
412 412 {
413 413 os<<v.str();
414 414 return os;
... ... @@ -417,9 +417,9 @@ std::ostream&amp; operator&lt;&lt;(std::ostream&amp; os, stim::vec&lt;T&gt; v)
417 417  
418 418 /// Multiply a vector by a constant when the vector is on the right hand side
419 419 template <typename T>
420   -stim::vec<T> operator*(T lhs, stim::vec<T> rhs)
  420 +tira::vec<T> operator*(T lhs, tira::vec<T> rhs)
421 421 {
422   - stim::vec<T> r;
  422 + tira::vec<T> r;
423 423  
424 424 return rhs * lhs;
425 425 }
... ...
stim/optics/beam.h 0 → 100644
  Show/Hide comments   View file @b28f312
  1 +#ifndef TIRA_BEAM
  2 +#define TIRA_BEAM
  3 +
  4 +#include "../math/vec3.h"
  5 +#include "../math/function.h"
  6 +#include "../optics/planewave.h"
  7 +#include <vector>
  8 +#include <numbers>
  9 +#include <random>
  10 +
  11 +namespace tira{
  12 + namespace optics{
  13 +template<typename T>
  14 +class beam {
  15 +
  16 + vec3<T> m_direction;
  17 + vec3<T> m_focus;
  18 + cvec3<T> m_E;
  19 + T m_k;
  20 +
  21 +
  22 +
  23 +public:
  24 + /// <summary>
  25 + /// Beam constructor
  26 + /// </summary>
  27 + /// <param name="d">Propagation direction of the beam (automatically normalized).</param>
  28 + /// <param name="f">Focal point for the beam.</param>
  29 + /// <param name="pol">Beam polarization (automatically normalized and orthogonalized). Only works for linear polarization, though.</param>
  30 + /// <param name="E">Complex beam amplitude (applied to the polarization to create a linearly polarized beam).</param>
  31 + /// <param name="lambda">Wavelength</param>
  32 + beam(
  33 + vec3<T> d = vec3<T>(0, 0, 1),
  34 + vec3<T> f = vec3<T>(0, 0, 0),
  35 + vec3<T> pol = vec3<T>(1, 0, 0),
  36 + std::complex<T> E = std::complex<T>(1.0, 0.0),
  37 + T lambda = 1.0
  38 + ) {
  39 +
  40 + m_direction = d.direction(); //normalize and set the plane wave direction
  41 + pol = (pol - d.dot(pol)).direction(); //orthogonalize and normalize the polarization vector
  42 + m_focus = f; //set the focal point
  43 + m_E[0] = pol[0] * E; //convert the polarization and complex amplitude to a linearly-polarized E vector
  44 + m_E[1] = pol[1] * E;
  45 + m_E[2] = pol[2] * E;
  46 + m_k = 2.0 * std::numbers::pi * lambda; //calculate and store the wavenumber
  47 + }
  48 +
  49 + /// <summary>
  50 + /// Generate a focal spot for the beam using Monte-Carlo sampling.
  51 + /// </summary>
  52 + /// <param name="NA">Numerical aperture of the focusing optics.</param>
  53 + /// <param name="samples">Number of Monte-Carlo samples.</param>
  54 + /// <param name="NA_in">NA of the internal obscuration.</param>
  55 + /// <returns></returns>
  56 + std::vector< planewave<T> > mcfocus(T NA, size_t samples, T NA_in = 0) {
  57 +
  58 + //create a rotation matrix to orient the beam along the specified direction vector
  59 + tira::matrix_sq<double, 3> R;// = tira::matrix_sq<double, 3>::identity(); //initialize the rotation matrix to the identity matrix
  60 + tira::vec3<double> Z(0, 0, 1); //create a vector along the Z direction
  61 + tira::quaternion<double> q;
  62 + q.CreateRotation(Z, m_direction);
  63 + R = q.toMatrix3();
  64 +
  65 + std::vector< planewave<T> > result(samples); //generate an array of N plane waves
  66 + double alpha = std::asin(NA); //calculate the angle of the band-pass
  67 + double alpha_in = std::asin(NA_in); //calculate the angle of the central obscuration
  68 + double cos_alpha = std::cos(alpha); //calculate the cosine of the band-pass and central obscuration
  69 + double cos_alpha_in = std::cos(alpha_in);
  70 +
  71 + //random sampling is done using Archimede's method - uniform samples are taken in cylindrical space and projected onto a sphere
  72 + std::random_device rd; //create a random number generator
  73 + std::default_random_engine eng(rd());
  74 + std::uniform_real_distribution<> theta_dist(0.0, 2.0 * std::numbers::pi); //create a distribution for theta (in cylindrical coordinates)
  75 + std::uniform_real_distribution<> z_dist(cos_alpha, cos_alpha_in); //create a uniform distribution for sampling the z-axis
  76 +
  77 + //generate a guide plane wave
  78 + planewave<T> pw(m_direction[0] * m_k, m_direction[1] * m_k, m_direction[2] * m_k, m_E[0], m_E[1], m_E[2]);
  79 +
  80 + double theta, phi; //variables to store the coordinates for theta/phi on a sphere
  81 + tira::vec3<T> dir;
  82 + T w = 1.0 / samples; //integration weight
  83 + for (size_t n = 0; n < samples; n++) {
  84 + theta = theta_dist(eng); //randomly generate a theta coordinate between 0 and 2pi
  85 + phi = std::acos(z_dist(eng)); //randomly generate a z coordinate based on the NA and project to phi
  86 + dir = R * tira::vec3<double>(1.0, theta, phi).sph2cart(); //convert the value to Cartesian coordinates and store the sample vector
  87 + result[n] = pw.refract(dir); //refract the plane wave to align with the sample direction
  88 + result[n] = result[n].translate(m_focus[0], m_focus[1], m_focus[2]); //refocus the plane wave to the focal position
  89 + result[n] = result[n] * w; //apply the integration weight (1/N)
  90 + }
  91 +
  92 + return result; //return the sample array
  93 + }
  94 +};
  95 +} //end namespace optics
  96 +} //end namespace tira
  97 +
  98 +#endif
... ...
stim/optics/planewave.h
  Show/Hide comments   View file @b28f312
1   -#ifndef STIM_PLANEWAVE_H
2   -#define STIM_PLANEWAVE_H
  1 +#ifndef TIRA_PLANEWAVE_H
  2 +#define TIRA_PLANEWAVE_H
3 3  
4 4 #include <string>
5 5 #include <sstream>
... ... @@ -11,29 +11,8 @@
11 11 #include "../math/plane.h"
12 12 #include <complex>
13 13  
14   -/// Calculate whether or not a vector v intersects the front (1) or back (-1) of a plane.
15   -/// This function returns -1 if the vector lies within the plane (is orthogonal to the surface normal)
16   -/*template <typename T>
17   -int plane_face(stim::vec3<T> v, stim::vec3<T> plane_normal) {
18   - T dotprod = v.dot(plane_normal); //calculate the dot product
19 14  
20   - if (dotprod < 0) return 1;
21   - if (dotprod > 0) return -1;
22   - return 0;
23   -}
24   -
25   -/// Calculate the component of a vector v that is perpendicular to a plane defined by a normal.
26   -template <typename T>
27   -stim::vec3<T> plane_perpendicular(stim::vec3<T> v, stim::vec3<T> plane_normal) {
28   - return plane_normal * v.dot(plane_normal);
29   -}
30   -
31   -template <typename T>
32   -stim::vec3<T> plane_parallel(stim::vec3<T> v, stim::vec3<T> plane_normal) {
33   - return v - plane_perpendicular(v, plane_normal);
34   -}*/
35   -
36   -namespace stim{
  15 +namespace tira{
37 16 namespace optics{
38 17  
39 18  
... ... @@ -77,121 +56,6 @@ namespace stim{
77 56 }
78 57  
79 58  
80   -/*template<typename T>
81   -class cvec3 {
82   -public:
83   - std::complex<T> x;
84   - std::complex<T> y;
85   - std::complex<T> z;
86   -
87   - cvec3(std::complex<T> _x, std::complex<T> _y, std::complex<T> _z) {
88   - x = _x;
89   - y = _y;
90   - z = _z;
91   - }
92   -};*/
93   -
94   -/*template<typename T>
95   -class cvec3 {
96   -public:
97   - std::complex<T> m_v[3];
98   -
99   - cvec3(std::complex<T> x, std::complex<T> y, std::complex<T> z) {
100   - m_v[0] = x;
101   - m_v[1] = y;
102   - m_v[2] = z;
103   - }
104   -
105   - cvec3() : cvec3(0, 0, 0) {}
106   -
107   - void operator()(std::complex<T> x, std::complex<T> y, std::complex<T> z) {
108   - m_v[0] = x;
109   - m_v[1] = y;
110   - m_v[2] = z;
111   - }
112   -
113   - std::complex<T> operator[](int i) const { return m_v[i]; }
114   - std::complex<T>& operator[](int i) { return m_v[i]; }
115   -
116   - /// Calculate the 2-norm of the complex vector
117   - T norm2() {
118   - T xx = std::real(m_v[0] * std::conj(m_v[0]));
119   - T yy = std::real(m_v[1] * std::conj(m_v[1]));
120   - T zz = std::real(m_v[2] * std::conj(m_v[2]));
121   - return std::sqrt(xx + yy + zz);
122   - }
123   -
124   - /// Returns the normalized direction vector
125   - cvec3 direction() {
126   - cvec3 result;
127   -
128   - std::complex<T> length = norm2();
129   - result[0] = m_v[0] / length;
130   - result[1] = m_v[1] / length;
131   - result[2] = m_v[2] / length;
132   -
133   - return result;
134   - }
135   -
136   - std::string str() {
137   - std::stringstream ss;
138   - ss << m_v[0] << ", " << m_v[1] << ", " << m_v[2];
139   - return ss.str();
140   - }
141   -
142   - //copy constructor
143   - cvec3(const cvec3 &other) {
144   - m_v[0] = other.m_v[0];
145   - m_v[1] = other.m_v[1];
146   - m_v[2] = other.m_v[2];
147   - }
148   -
149   - /// Assignment operator
150   - cvec3& operator=(const cvec3 &rhs) {
151   - m_v[0] = rhs.m_v[0];
152   - m_v[1] = rhs.m_v[1];
153   - m_v[2] = rhs.m_v[2];
154   -
155   - return *this;
156   - }
157   -
158   - /// Calculate and return the cross product between this vector and another
159   - cvec3 cross(cvec3 rhs) {
160   - cvec3 result;
161   -
162   - //compute the cross product (only valid for 3D vectors)
163   - result[0] = (m_v[1] * rhs[2] - m_v[2] * rhs[1]);
164   - result[1] = (m_v[2] * rhs[0] - m_v[0] * rhs[2]);
165   - result[2] = (m_v[1] * rhs[1] - m_v[1] * rhs[0]);
166   -
167   - return result;
168   - }
169   -
170   - /// Calculate and return the dot product between this vector and another
171   - std::complex<T> dot(cvec3 rhs) {
172   - return m_v[0] * rhs[0] + m_v[1] * rhs[1] + m_v[2] * rhs[2];
173   - }
174   -
175   - /// Arithmetic multiplication: returns the vector scaled by a constant value
176   - template<typename R>
177   - cvec3 operator*(R rhs) const
178   - {
179   - cvec3 result;
180   - result[0] = m_v[0] * rhs;
181   - result[1] = m_v[1] * rhs;
182   - result[2] = m_v[2] * rhs;
183   -
184   - return result;
185   - }
186   -
187   -};
188   -template <typename T>
189   -std::ostream& operator<<(std::ostream& os, stim::optics::cvec3<T> p) {
190   - os << p.str();
191   - return os;
192   -}
193   -*/
194   -
195 59 template<typename T>
196 60 class planewave{
197 61  
... ... @@ -201,12 +65,12 @@ protected:
201 65 cvec3<T> m_E; //amplitude (for a scalar plane wave, only E0[0] is used)
202 66  
203 67 /// Bend a plane wave via refraction, given that the new propagation direction is known
204   - CUDA_CALLABLE planewave<T> bend(stim::vec3<T> v) const {
  68 + CUDA_CALLABLE planewave<T> bend(vec3<T> v) const {
205 69  
206   - stim::vec3<T> k_real(m_k.get(0).real(), m_k.get(1).real(), m_k.get(2).real()); //force the vector to be real (can only refract real directions)
  70 + vec3<T> k_real(m_k.get(0).real(), m_k.get(1).real(), m_k.get(2).real()); //force the vector to be real (can only refract real directions)
207 71  
208   - stim::vec3<T> kn_hat = v.direction(); //normalize the new k
209   - stim::vec3<T> k_hat = k_real.direction(); //normalize the current k
  72 + vec3<T> kn_hat = v.direction(); //normalize the new k
  73 + vec3<T> k_hat = k_real.direction(); //normalize the current k
210 74  
211 75 planewave<T> new_p; //create a new plane wave
212 76  
... ... @@ -227,10 +91,10 @@ protected:
227 91 }
228 92  
229 93 vec3<T> r = k_hat.cross(kn_hat); //compute the rotation vector
230   - T theta = asin(r.len()); //compute the angle of the rotation about r
231   - quaternion<T> q; //create a quaternion to capture the rotation
  94 + T theta = asin(r.len()); //compute the angle of the rotation about r
  95 + quaternion<T> q; //create a quaternion to capture the rotation
232 96 q.CreateRotation(theta, r.direction());
233   - stim::matrix_sq<T, 3> R = q.toMatrix3();
  97 + matrix_sq<T, 3> R = q.toMatrix3();
234 98 vec3< std::complex<T> > E(m_E.get(0), m_E.get(1), m_E.get(2));
235 99 vec3< std::complex<T> > E0n = R * E; //apply the rotation to E0
236 100 //new_p.m_k = kn_hat * kmag();
... ... @@ -252,8 +116,6 @@ public:
252 116  
253 117  
254 118 ///constructor: create a plane wave propagating along k
255   - //CUDA_CALLABLE planewave(vec<T> kvec = stim::vec<T>(0, 0, stim::TAU),
256   - // vec< complex<T> > E = stim::vec<T>(1, 0, 0))
257 119 CUDA_CALLABLE planewave(std::complex<T> kx, std::complex<T> ky, std::complex<T> kz,
258 120 std::complex<T> Ex, std::complex<T> Ey, std::complex<T> Ez) {
259 121  
... ... @@ -313,22 +175,10 @@ public:
313 175 return result;
314 176 }
315 177  
316   - /*int scatter(vec3<T> surface_normal, vec3<T> surface_point, T ni, std::complex<T> nt,
317   - planewave& Pi, planewave& Pr, planewave& Pt) {
318   -
319   -
320   - return 0;
321   -
322   - }*/
323   -
324 178 CUDA_CALLABLE T kmag() const {
325 179 return std::sqrt(std::real(m_k.get(0) * std::conj(m_k.get(0)) + m_k.get(1) * std::conj(m_k.get(1)) + m_k.get(2) * std::conj(m_k.get(2))));
326 180 }
327 181  
328   - CUDA_CALLABLE planewave<T> refract(stim::vec3<T> kn) const {
329   - return bend(kn);
330   - }
331   -
332 182 /// Return a plane wave with the origin translated by (x, y, z)
333 183 CUDA_CALLABLE planewave<T> translate(T x, T y, T z) const {
334 184 planewave<T> result;
... ... @@ -350,6 +200,17 @@ public:
350 200 return *this;
351 201 }
352 202  
  203 + ///return a plane wave with the applied refractive index (scales the k-vector by the input)
  204 + CUDA_CALLABLE planewave<T> ri(T n) {
  205 + planewave<T> result;
  206 + result.m_E = m_E;
  207 + result.m_k = m_k * n;
  208 + return result;
  209 + }
  210 + CUDA_CALLABLE planewave<T> refract(vec3<T> kn) const {
  211 + return bend(kn);
  212 + }
  213 +
353 214 /*CUDA_CALLABLE T lambda() const{
354 215 return stim::TAU / k.len();
355 216 }
... ... @@ -390,9 +251,7 @@ public:
390 251 return planewave<T>(k * (nt / ni), E0);
391 252 }
392 253  
393   - CUDA_CALLABLE planewave<T> refract(stim::vec<T> kn) const{
394   - return bend(kn);
395   - }
  254 +
396 255  
397 256 /// Calculate the result of a plane wave hitting an interface between two refractive indices
398 257  
... ... @@ -408,29 +267,28 @@ public:
408 267 /// Calculate the scattering result when nr = n1/n0
409 268  
410 269 /// @param P is a plane representing the position and orientation of the surface
411   - /// @param r is the ration n1/n0
  270 + /// @param r is the ratio n1/n0
412 271 /// @param n1 is the refractive index inside the surface (in the direction away from the normal)
413 272 /// @param r is the reflected component of the plane wave
414 273 /// @param t is the transmitted component of the plane wave
415 274  
416   - void scatter(vec3<T> plane_normal, vec3<T> plane_position, std::complex<T> nr, planewave<T>& r, planewave<T>& t) {
  275 + int scatter(vec3<T> plane_normal, vec3<T> plane_position, std::complex<T> nr, planewave<T>& r, planewave<T>& t) {
417 276  
418 277 if (m_k[0].imag() != 0.0 || m_k[1].imag() != 0.0 || m_k[2].imag() != 0) {
419 278 std::cout << "ERROR: cannot scatter a plane wave with an imaginary k-vector." << std::endl;
420 279 }
421 280  
422   - stim::vec3<T> ki(m_k[0].real(), m_k[1].real(), m_k[2].real()); //force the current k vector to be real
423   - stim::vec3<T> kr;
424   - stim::cvec3<T> kt, Ei, Er, Et;
  281 + vec3<T> ki(m_k[0].real(), m_k[1].real(), m_k[2].real()); //force the current k vector to be real
  282 + vec3<T> kr;
  283 + cvec3<T> kt, Ei, Er, Et;
425 284  
426 285 plane_normal = plane_normal.direction();
427   - stim::vec3<T> k_dir = ki.direction(); // calculate the direction of the incident plane wave
428   -
429   - int facing = plane_face(k_dir, plane_normal); //determine which direction the plane wave is coming in
  286 + vec3<T> k_dir = ki.direction(); //calculate the direction of the incident plane wave
430 287  
431   - if (facing == -1) { //if the wave hits the back of the plane, invert the plane and nr
  288 + //int facing = plane_face(k_dir, plane_normal); //determine which direction the plane wave is coming in
  289 + if (k_dir.dot(plane_normal) > 0) { //if the wave hits the back of the plane, invert the plane and nr
432 290 std::cout << "ERROR: k-vector intersects the wrong side of the boundary." << std::endl;
433   - return -1; //the plane wave is impacting the wrong side of the surface
  291 + return -1; //the plane wave is impacting the wrong side of the surface
434 292 }
435 293  
436 294 //use Snell's Law to calculate the transmitted angle
... ... @@ -443,17 +301,15 @@ public:
443 301 std::complex<T> rp = (1.0 - nr) / (1.0 + nr); //compute the Fresnel coefficients
444 302 std::complex<T> tp = 2.0 / (1.0 + nr);
445 303  
446   - //ki = kv * ni; //calculate the incident k-vector (scaled by the incident refractive index)
447 304 kr = -ki; //the reflection vector is the inverse of the incident vector
448 305 kt[0] = ki[0] * nr;
449 306 kt[1] = ki[1] * nr;
450 307 kt[2] = ki[2] * nr;
451   - //Ei = Ev; //compute the E vectors for all three plane waves based on the Fresnel coefficients
452   - Er = m_E * rp;
  308 +
  309 + Er = m_E * rp; //compute the E vectors based on the Fresnel coefficients
453 310 Et = m_E * tp;
454 311  
455 312 //calculate the phase offset based on the plane positions
456   - //T phase_i = plane_position.dot(kv - ki); //compute the phase offset for each plane wave
457 313 T phase_r = plane_position.dot(ki - kr);
458 314 std::complex<T> phase_t =
459 315 plane_position[0] * (ki[0] - kt[0]) +
... ... @@ -465,25 +321,25 @@ public:
465 321 T sin_theta_r = sin_theta_i;
466 322 T theta_r = theta_i;
467 323  
468   - std::complex<T> sin_theta_t = nr * sin(theta_i); //compute the sine of theta_t using Snell's law
  324 + std::complex<T> sin_theta_t = (1.0/nr) * sin(theta_i); //compute the sine of theta_t using Snell's law
469 325 std::complex<T> cos_theta_t = std::sqrt(1.0 - sin_theta_t * sin_theta_t);
470 326 std::complex<T> theta_t = asin(sin_theta_t); //compute the cosine of theta_t
471 327  
472 328 //Define the basis vectors for the calculation (plane of incidence)
473   - stim::vec3<T> z_hat = -plane_normal;
474   - stim::vec3<T> y_hat = plane_parallel(k_dir, plane_normal);
475   - stim::vec3<T> x_hat = y_hat.cross(z_hat);
  329 + vec3<T> z_hat = -plane_normal;
  330 + vec3<T> plane_perpendicular = plane_normal * k_dir.dot(plane_normal);
  331 + vec3<T> y_hat = (k_dir - plane_perpendicular).direction();
  332 + vec3<T> x_hat = y_hat.cross(z_hat);
476 333  
477 334 //calculate the k-vector magnitudes
478   - //T kv_mag = kv.norm2();
479 335 T ki_mag = ki.norm2();
480 336 T kr_mag = ki_mag;
481 337 std::complex<T> kt_mag = ki_mag * nr;
482 338  
483 339 //calculate the k vector directions
484   - stim::vec3<T> ki_dir = y_hat * sin_theta_i + z_hat * cos_theta_i;
485   - stim::vec3<T> kr_dir = y_hat * sin_theta_r - z_hat * cos_theta_r;
486   - stim::cvec3<T> kt_dir;
  340 + vec3<T> ki_dir = y_hat * sin_theta_i + z_hat * cos_theta_i;
  341 + vec3<T> kr_dir = y_hat * sin_theta_r - z_hat * cos_theta_r;
  342 + cvec3<T> kt_dir;
487 343 kt_dir[0] = y_hat[0] * sin_theta_t + z_hat[0] * cos_theta_t;
488 344 kt_dir[1] = y_hat[1] * sin_theta_t + z_hat[1] * cos_theta_t;
489 345 kt_dir[2] = y_hat[2] * sin_theta_t + z_hat[2] * cos_theta_t;
... ... @@ -500,34 +356,22 @@ public:
500 356 std::complex<T> tp = ((2.0 * sin_theta_t * cos_theta_i) / (std::sin(theta_t + theta_i) * std::cos(theta_t - theta_i)));
501 357  
502 358 //calculate the p component directions for each E vector
503   - stim::vec3<T> Eip_dir = y_hat * cos_theta_i - z_hat * sin_theta_i;
504   - stim::vec3<T> Erp_dir = y_hat * cos_theta_r + z_hat * sin_theta_r;
505   - stim::cvec3<T> Etp_dir;
  359 + vec3<T> Eip_dir = y_hat * cos_theta_i - z_hat * sin_theta_i;
  360 + vec3<T> Erp_dir = y_hat * cos_theta_r + z_hat * sin_theta_r;
  361 + cvec3<T> Etp_dir;
506 362 Etp_dir[0] = y_hat[0] * cos_theta_t - z_hat[0] * sin_theta_t;
507 363 Etp_dir[1] = y_hat[1] * cos_theta_t - z_hat[1] * sin_theta_t;
508 364 Etp_dir[2] = y_hat[2] * cos_theta_t - z_hat[2] * sin_theta_t;
509 365  
510 366 //calculate the s and t components of each E vector
511   - //std::complex<T> E_mag = Ev.norm2();
512 367 std::complex<T> Ei_s = m_E.dot(x_hat);
513   - //std::complex<T> Ei_p = std::sqrt(E_mag * E_mag - Ei_s * Ei_s);
514 368 std::complex<T> Ei_p = m_E.dot(Eip_dir);
515 369 std::complex<T> Er_s = rs * Ei_s;
516 370 std::complex<T> Er_p = rp * Ei_p;
517 371 std::complex<T> Et_s = ts * Ei_s;
518 372 std::complex<T> Et_p = tp * Ei_p;
519 373  
520   -
521   -
522 374 //calculate the E vector for each plane wave
523   - //Ei[0] = Eip_dir[0] * Ei_p + x_hat[0] * Ei_s;
524   - //Ei[1] = Eip_dir[1] * Ei_p + x_hat[1] * Ei_s;
525   - //Ei[2] = Eip_dir[2] * Ei_p + x_hat[2] * Ei_s;
526   -
527   - //std::cout << Ev << std::endl;
528   - //std::cout << Ei << std::endl;
529   - //std::cout << std::endl;
530   -
531 375 Er[0] = Erp_dir[0] * Er_p + x_hat[0] * Er_s;
532 376 Er[1] = Erp_dir[1] * Er_p + x_hat[1] * Er_s;
533 377 Er[2] = Erp_dir[2] * Er_p + x_hat[2] * Er_s;
... ... @@ -539,134 +383,21 @@ public:
539 383  
540 384  
541 385 //calculate the phase offset based on the plane positions
542   - //T phase_i = plane_position.dot(kv - ki); //compute the phase offset for each plane wave
543 386 T phase_r = plane_position.dot(ki - kr);
544 387 std::complex<T> phase_t =
545 388 plane_position[0] * (ki[0] - kt[0]) +
546 389 plane_position[1] * (ki[1] - kt[1]) +
547 390 plane_position[2] * (ki[2] - kt[2]);
548 391  
549   -
550   - //Ei = Ei * std::exp(std::complex<T>(0, 1));
  392 + //apply the phase offset
551 393 Er = Er * std::exp(std::complex<T>(0, 1) * phase_r);
552 394 Et = Et * std::exp(std::complex<T>(0, 1) * phase_t);
553 395  
554   - //Pi = stim::optics::planewave<T>(ki[0], ki[1], ki[2], Ei[0], Ei[1], Ei[2]);
555   - r = stim::optics::planewave<T>(kr[0], kr[1], kr[2], Er[0], Er[1], Er[2]);
556   - t = stim::optics::planewave<T>(kt[0], kt[1], kt[2], Et[0], Et[1], Et[2]);
  396 + //generate the reflected and transmitted waves
  397 + r = planewave<T>(kr[0], kr[1], kr[2], Er[0], Er[1], Er[2]);
  398 + t = planewave<T>(kt[0], kt[1], kt[2], Et[0], Et[1], Et[2]);
557 399  
558 400 return 0;
559   - /*TODO: Generate a real vector from the current K vector - this will not address complex k-vectors
560   - vec3< std::complex<T> > k(3);
561   - k[0] = m_k[0];
562   - k[1] = m_k[1];
563   - k[2] = m_k[2];
564   -
565   - vec3< std::complex<T> > E(3);
566   - E[0] = m_E[0];// std::complex<T>(m_E[0].real(), m_E[0].imag());
567   - E[1] = m_E[1];// std::complex<T>(m_E[1].real(), m_E[1].imag());
568   - E[2] = m_E[2];// std::complex<T>(m_E[2].real(), m_E[2].imag());
569   -
570   - std::complex<T> cOne(1, 0);
571   - std::complex<T> cTwo(2, 0);
572   -
573   - T kmag = m_k.norm2();
574   -
575   - int facing = plane_face(k, plane_normal); //determine which direction the plane wave is coming in
576   -
577   - if(facing == -1){ //if the wave hits the back of the plane, invert the plane and nr
578   - plane_normal = -plane_normal; //flip the plane
579   - nr = cOne / nr; //invert the refractive index (now nr = n0/n1)
580   - }
581   -
582   - //use Snell's Law to calculate the transmitted angle
583   - //vec3<T> plane_normal = -P.n();
584   - T cos_theta_i = k.norm().dot(-plane_normal); //compute the cosine of theta_i
585   - T theta_i = acos(cos_theta_i); //compute theta_i
586   - std::complex<T> sin_theta_t = (cOne / nr) * sin(theta_i); //compute the sine of theta_t using Snell's law
587   - std::complex<T> theta_t = asin(sin_theta_t); //compute the cosine of theta_t
588   -
589   - //bool tir = false; //flag for total internal reflection
590   - //if(theta_t != theta_t){
591   - // tir = true;
592   - // theta_t = stim::PI / (T)2;
593   - //}
594   -
595   - //handle the degenerate case where theta_i is 0 (the plane wave hits head-on)
596   - if(theta_i == 0){
597   - std::complex<T> rp = (cOne - nr) / (cOne + nr); //compute the Fresnel coefficients
598   - std::complex<T> tp = (cOne * cTwo) / (cOne + nr);
599   - vec3<T> kr = -k;
600   - vec3<T> kt = k * nr; //set the k vectors for theta_i = 0
601   - vec3< std::complex<T> > Er = E * rp; //compute the E vectors
602   - vec3< std::complex<T> > Et = E * tp;
603   - T phase_t = plane_position.dot(k - kt); //compute the phase offset
604   - T phase_r = plane_position.dot(k - kr);
605   -
606   - //create the plane waves
607   - //r = planewave<T>(kr, Er, phase_r);
608   - //t = planewave<T>(kt, Et, phase_t);
609   -
610   - r = planewave(kr[0], kr[1], kr[2], Er[0], Er[1], Er[2]);
611   - t = planewave(kt[0], kt[1], kt[2], Et[0], Et[1], Et[2]);
612   -
613   - return;
614   - }
615   -
616   -
617   - //compute the Fresnel coefficients
618   - T rp, rs, tp, ts;
619   - rp = tan(theta_t - theta_i) / tan(theta_t + theta_i);
620   - rs = sin(theta_t - theta_i) / sin(theta_t + theta_i);
621   -
622   - //if (tir) {
623   - // tp = ts = 0;
624   - //}
625   - //else
626   - {
627   - tp = ( 2 * sin(theta_t) * cos(theta_i) ) / ( sin(theta_t + theta_i) * cos(theta_t - theta_i) );
628   - ts = ( 2 * sin(theta_t) * cos(theta_i) ) / sin(theta_t + theta_i);
629   - }
630   -
631   - //compute the coordinate space for the plane of incidence
632   - vec3<T> z_hat = -plane_normal;
633   - vec3<T> y_hat = plane_parallel(k, plane_normal).norm();
634   - vec3<T> x_hat = y_hat.cross(z_hat).norm();
635   -
636   - //compute the k vectors for r and t
637   - vec3<T> kr, kt;
638   - kr = ( y_hat * sin(theta_i) - z_hat * cos(theta_i) ) * kmag;
639   - kt = ( y_hat * sin(theta_t) + z_hat * cos(theta_t) ) * kmag * nr;
640   -
641   - //compute the magnitude of the p- and s-polarized components of the incident E vector
642   - std::complex<T> Ei_s = E.dot(x_hat);
643   - //int sign = sgn(E0.dot(y_hat));
644   - vec3< std::complex<T> > cx_hat = x_hat;
645   - std::complex<T> Ei_p = (E - cx_hat * Ei_s).len();// *(T)sign;
646   - //compute the magnitude of the p- and s-polarized components of the reflected E vector
647   - std::complex<T> Er_s = Ei_s * rs;
648   - std::complex<T> Er_p = Ei_p * rp;
649   - //compute the magnitude of the p- and s-polarized components of the transmitted E vector
650   - std::complex<T> Et_s = Ei_s * ts;
651   - std::complex<T> Et_p = Ei_p * tp;
652   -
653   - //compute the reflected E vector
654   - vec3< std::complex<T> > Er = vec3< std::complex<T> >(y_hat * cos(theta_i) + z_hat * sin(theta_i)) * Er_p + cx_hat * Er_s;
655   - //compute the transmitted E vector
656   - vec3< std::complex<T> > Et = vec3< std::complex<T> >(y_hat * cos(theta_t) - z_hat * sin(theta_t)) * Et_p + cx_hat * Et_s;
657   -
658   - T phase_t = plane_position.dot(k - kt);
659   - T phase_r = plane_position.dot(k - kr);
660   -
661   - //create the plane waves
662   - //r.m_k = kr;
663   - //r.m_E = Er * exp( complex<T>(0, phase_r) );
664   -
665   - //t.m_k = kt;
666   - //t.m_E = Et * exp( complex<T>(0, phase_t) );
667   -
668   - r = planewave(kr[0], kr[1], kr[2], Er[0], Er[1], Er[2]);
669   - t = planewave(kt[0], kt[1], kt[2], Et[0], Et[1], Et[2]);*/
670 401 }
671 402  
672 403 std::string str()
... ... @@ -678,10 +409,10 @@ public:
678 409 }
679 410 }; //end planewave class
680 411 } //end namespace optics
681   -} //end namespace stim
  412 +} //end namespace tira
682 413  
683 414 template <typename T>
684   -std::ostream& operator<<(std::ostream& os, stim::optics::planewave<T> p)
  415 +std::ostream& operator<<(std::ostream& os, tira::optics::planewave<T> p)
685 416 {
686 417 os<<p.str();
687 418 return os;
... ...