codebase / stimlib

Commit 559d0fcbc33712d4de94f257b6b2718ab9f5c729

Authored by David Mayerich
1 parent d609550e

wave interactions with a planar surface

Showing 12 changed files with 851 additions and 213 deletions   Show diff stats
Inline Side-by-side
math/bessel.h
  Show/Hide comments   View file @559d0fc
... ... @@ -759,28 +759,28 @@ int bessjyv(P v,P x,P &vm,P *jv,P *yv,
759 759 }
760 760 vm = n + v0;
761 761 return 0;
762   -}
763   -
764   -template<typename P>
  762 +}
  763 +
  764 +template<typename P>
765 765 int bessjyv_sph(int v, P z, P &vm, P* cjv,
766   - P* cyv, P* cjvp, P* cyvp)
767   -{
768   - //first, compute the bessel functions of fractional order
769   - bessjyv(v + 0.5, z, vm, cjv, cyv, cjvp, cyvp);
770   -
771   - //iterate through each and scale
772   - for(int n = 0; n<=v; n++)
773   - {
774   -
775   - cjv[n] = cjv[n] * sqrt(PI/(z * 2.0));
776   - cyv[n] = cyv[n] * sqrt(PI/(z * 2.0));
777   -
778   - cjvp[n] = -1.0 / (z * 2.0) * cjv[n] + cjvp[n] * sqrt(PI / (z * 2.0));
779   - cyvp[n] = -1.0 / (z * 2.0) * cyv[n] + cyvp[n] * sqrt(PI / (z * 2.0));
780   - }
781   -
782   - return 0;
783   -
  766 + P* cyv, P* cjvp, P* cyvp)
  767 +{
  768 + //first, compute the bessel functions of fractional order
  769 + bessjyv(v + 0.5, z, vm, cjv, cyv, cjvp, cyvp);
  770 +
  771 + //iterate through each and scale
  772 + for(int n = 0; n<=v; n++)
  773 + {
  774 +
  775 + cjv[n] = cjv[n] * sqrt(rtsPI/(z * 2.0));
  776 + cyv[n] = cyv[n] * sqrt(rtsPI/(z * 2.0));
  777 +
  778 + cjvp[n] = -1.0 / (z * 2.0) * cjv[n] + cjvp[n] * sqrt(rtsPI / (z * 2.0));
  779 + cyvp[n] = -1.0 / (z * 2.0) * cyv[n] + cyvp[n] * sqrt(rtsPI / (z * 2.0));
  780 + }
  781 +
  782 + return 0;
  783 +
784 784 }
785 785  
786 786 template<typename P>
... ... @@ -1485,28 +1485,28 @@ int cbessjyva(P v,complex&lt;P&gt; z,P &amp;vm,complex&lt;P&gt;*cjv,
1485 1485 }
1486 1486 vm = n+v0;
1487 1487 return 0;
1488   -}
1489   -
1490   -template<typename P>
  1488 +}
  1489 +
  1490 +template<typename P>
1491 1491 int cbessjyva_sph(int v,complex<P> z,P &vm,complex<P>*cjv,
1492   - complex<P>*cyv,complex<P>*cjvp,complex<P>*cyvp)
1493   -{
1494   - //first, compute the bessel functions of fractional order
1495   - cbessjyva<P>(v + 0.5, z, vm, cjv, cyv, cjvp, cyvp);
1496   -
1497   - //iterate through each and scale
1498   - for(int n = 0; n<=v; n++)
1499   - {
1500   -
1501   - cjv[n] = cjv[n] * sqrt(PI/(z * 2.0));
1502   - cyv[n] = cyv[n] * sqrt(PI/(z * 2.0));
1503   -
1504   - cjvp[n] = -1.0 / (z * 2.0) * cjv[n] + cjvp[n] * sqrt(PI / (z * 2.0));
1505   - cyvp[n] = -1.0 / (z * 2.0) * cyv[n] + cyvp[n] * sqrt(PI / (z * 2.0));
1506   - }
1507   -
1508   - return 0;
1509   -
  1492 + complex<P>*cyv,complex<P>*cjvp,complex<P>*cyvp)
  1493 +{
  1494 + //first, compute the bessel functions of fractional order
  1495 + cbessjyva<P>(v + 0.5, z, vm, cjv, cyv, cjvp, cyvp);
  1496 +
  1497 + //iterate through each and scale
  1498 + for(int n = 0; n<=v; n++)
  1499 + {
  1500 +
  1501 + cjv[n] = cjv[n] * sqrt(rtsPI/(z * 2.0));
  1502 + cyv[n] = cyv[n] * sqrt(rtsPI/(z * 2.0));
  1503 +
  1504 + cjvp[n] = -1.0 / (z * 2.0) * cjv[n] + cjvp[n] * sqrt(rtsPI / (z * 2.0));
  1505 + cyvp[n] = -1.0 / (z * 2.0) * cyv[n] + cyvp[n] * sqrt(rtsPI / (z * 2.0));
  1506 + }
  1507 +
  1508 + return 0;
  1509 +
1510 1510 }
1511 1511  
1512 1512 } //end namespace rts
... ...
math/complex.h
  Show/Hide comments   View file @559d0fc
... ... @@ -23,7 +23,14 @@ struct complex
23 23 CUDA_CALLABLE complex()
24 24 {
25 25 r = 0;
26   - i = 0;
  26 + i = 0;
  27 + }
  28 +
  29 + //constructor when given real and imaginary values
  30 + CUDA_CALLABLE complex(T r, T i = 0)
  31 + {
  32 + this->r = r;
  33 + this->i = i;
27 34 }
28 35  
29 36 //access methods
... ... @@ -48,12 +55,7 @@ struct complex
48 55 return i_val;
49 56 }
50 57  
51   - //constructor when given real and imaginary values
52   - CUDA_CALLABLE complex(T r, T i)
53   - {
54   - this->r = r;
55   - this->i = i;
56   - }
  58 +
57 59  
58 60 //return the current value multiplied by i
59 61 CUDA_CALLABLE complex<T> imul()
... ...
math/constants.h
  Show/Hide comments   View file @559d0fc
1 1 #ifndef RTS_CONSTANTS_H
2 2 #define RTS_CONSTANTS_H
3 3  
4   -namespace rts{
5   -
6   - double PI = 3.14159;
7   - double TAU = 2 * PI;
8   -}
  4 +#define rtsPI 3.14159
  5 +#define rtsTAU 2 * rtsPI
9 6  
10 7 #endif
11 8 \ No newline at end of file
... ...
math/plane.h 0 → 100644
  Show/Hide comments   View file @559d0fc
  1 +#ifndef RTS_PLANE_H
  2 +#define RTS_PLANE_H
  3 +
  4 +#include <iostream>
  5 +#include "../math/vector.h"
  6 +#include "rts/cuda/callable.h"
  7 +
  8 +
  9 +namespace rts{
  10 +template <typename T, int D> class plane;
  11 +}
  12 +
  13 +template <typename T, int D>
  14 +CUDA_CALLABLE rts::plane<T, D> operator-(rts::plane<T, D> v);
  15 +
  16 +namespace rts{
  17 +
  18 +template <class T, int D = 3>
  19 +class plane{
  20 +
  21 + //a plane is defined by a point and a normal
  22 +
  23 +private:
  24 +
  25 + vec<T, D> P; //point on the plane
  26 + vec<T, D> N; //plane normal
  27 +
  28 + CUDA_CALLABLE void init(){
  29 + P = vec<T, D>(0, 0, 0);
  30 + N = vec<T, D>(0, 0, 1);
  31 + }
  32 +
  33 +
  34 +public:
  35 +
  36 + //default constructor
  37 + CUDA_CALLABLE plane(){
  38 + init();
  39 + }
  40 +
  41 + CUDA_CALLABLE plane(vec<T, D> n, vec<T, D> p = vec<T, D>(0, 0, 0)){
  42 + P = p;
  43 + N = n.norm();
  44 + }
  45 +
  46 + CUDA_CALLABLE plane(T z_pos){
  47 + init();
  48 + P[2] = z_pos;
  49 + }
  50 +
  51 + //create a plane from three points (a triangle)
  52 + CUDA_CALLABLE plane(vec<T, D> a, vec<T, D> b, vec<T, D> c){
  53 + P = c;
  54 + N = (c - a).cross(b - a);
  55 + if(N.len() == 0) //handle the degenerate case when two vectors are the same, N = 0
  56 + N = 0;
  57 + else
  58 + N = N.norm();
  59 + }
  60 +
  61 + CUDA_CALLABLE vec<T, D> norm(){
  62 + return N;
  63 + }
  64 +
  65 + CUDA_CALLABLE vec<T, D> p(){
  66 + return P;
  67 + }
  68 +
  69 + //flip the plane front-to-back
  70 + CUDA_CALLABLE plane<T, D> flip(){
  71 + plane<T, D> result = *this;
  72 + result.N = -result.N;
  73 + return result;
  74 + }
  75 +
  76 + //determines how a vector v intersects the plane (1 = intersects front, 0 = within plane, -1 = intersects back)
  77 + CUDA_CALLABLE int face(vec<T, D> v){
  78 +
  79 + T dprod = v.dot(N); //get the dot product between v and N
  80 +
  81 + //conditional returns the appropriate value
  82 + if(dprod < 0)
  83 + return 1;
  84 + else if(dprod > 0)
  85 + return -1;
  86 + else
  87 + return 0;
  88 + }
  89 +
  90 + //determine on which side of the plane a point lies (1 = front, 0 = on the plane, -1 = back)
  91 + CUDA_CALLABLE int side(vec<T, D> p){
  92 +
  93 + vec<T, D> v = p - P; //get the vector from P to the query point p
  94 +
  95 + return face(v);
  96 + }
  97 +
  98 + //compute the component of v that is perpendicular to the plane
  99 + CUDA_CALLABLE vec<T, D> perpendicular(vec<T, D> v){
  100 + return N * v.dot(N);
  101 + }
  102 +
  103 + //compute the projection of v in the plane
  104 + CUDA_CALLABLE vec<T, D> parallel(vec<T, D> v){
  105 + return v - perpendicular(v);
  106 + }
  107 +
  108 + //get both the parallel and perpendicular components of a vector v w.r.t. the plane
  109 + CUDA_CALLABLE void project(vec<T, D> v, vec<T, D> &v_par, vec<T, D> &v_perp){
  110 +
  111 + v_perp = v.dot(N);
  112 + v_par = v - v_perp;
  113 + }
  114 +
  115 + //compute the reflection of v off of the plane
  116 + CUDA_CALLABLE vec<T, D> reflect(vec<T, D> v){
  117 +
  118 + //compute the reflection using N_prime as the plane normal
  119 + vec<T, D> par = parallel(v);
  120 + vec<T, D> r = (-v) + par * 2;
  121 +
  122 + /*std::cout<<"----------------REFLECT-----------------------------"<<std::endl;
  123 + std::cout<<str()<<std::endl;
  124 + std::cout<<"v: "<<v<<std::endl;
  125 + std::cout<<"r: "<<r<<std::endl;
  126 + std::cout<<"Perpendicular: "<<perpendicular(v)<<std::endl;
  127 + std::cout<<"Parallel: "<<par<<std::endl;*/
  128 + return r;
  129 +
  130 + }
  131 +
  132 + //output a string
  133 + std::string str(){
  134 + std::stringstream ss;
  135 + ss<<"P: "<<P<<std::endl;
  136 + ss<<"N: "<<N;
  137 + return ss.str();
  138 + }
  139 +
  140 + ///////Friendship
  141 + friend CUDA_CALLABLE rts::plane<T, D> operator- <> (rts::plane<T, D> v);
  142 +
  143 +
  144 +
  145 +};
  146 +
  147 +}
  148 +
  149 +//arithmetic operators
  150 +
  151 +//negative operator flips the plane (front to back)
  152 +template <typename T, int D>
  153 +CUDA_CALLABLE rts::plane<T, D> operator-(rts::plane<T, D> p_rhs)
  154 +{
  155 + rts::plane<T, D> p = p_rhs;
  156 +
  157 + //negate the normal vector
  158 + p.N = -p.N;
  159 +
  160 + return p;
  161 +}
  162 +
  163 +
  164 +
  165 +#endif
0 166 \ No newline at end of file
... ...
math/realfield.cuh
  Show/Hide comments   View file @559d0fc
... ... @@ -84,7 +84,7 @@ public:
84 84 {
85 85 R[0] = R[1] = 0;
86 86 init();
87   - std::cout<<"realfield CONSTRUCTOR"<<std::endl;
  87 + //std::cout<<"realfield CONSTRUCTOR"<<std::endl;
88 88 }
89 89 realfield(unsigned int x, unsigned int y)
90 90 {
... ... @@ -98,13 +98,13 @@ public:
98 98 }
99 99 //shape = rect<P>(vec<P>(-1, -1, 0), vec<P>(-1, 1, 0), vec<P>(1, 1, 0)); //default geometry
100 100 clear(); //zero the field
101   - std::cout<<"realfield CONSTRUCTOR"<<std::endl;
  101 + //std::cout<<"realfield CONSTRUCTOR"<<std::endl;
102 102 }
103 103  
104 104 ~realfield()
105 105 {
106 106 destroy();
107   - std::cout<<"realfield DESTRUCTOR"<<std::endl;
  107 + //std::cout<<"realfield DESTRUCTOR"<<std::endl;
108 108 }
109 109  
110 110 P* ptr(unsigned int n)
... ... @@ -207,7 +207,7 @@ public:
207 207 //copy the slice
208 208 HANDLE_ERROR(cudaMemcpy(X[n], rhs.X[n], sizeof(P) * R[0] * R[1], cudaMemcpyDeviceToDevice));
209 209 }
210   - std::cout<<"Assignment operator."<<std::endl;
  210 + //std::cout<<"Assignment operator."<<std::endl;
211 211  
212 212 return *this;
213 213 }
... ... @@ -234,7 +234,7 @@ public:
234 234 }
235 235 }
236 236  
237   - std::cout<<"realfield COPY CONSTRUCTOR"<<std::endl;
  237 + //std::cout<<"realfield COPY CONSTRUCTOR"<<std::endl;
238 238 }
239 239  
240 240 /*void gaussian(P mean, P std, unsigned int n=0) //creates a 3D gaussian using component n
... ...
math/vector.h
  Show/Hide comments   View file @559d0fc
... ... @@ -26,9 +26,10 @@ struct vec
26 26 }
27 27  
28 28 //efficiency constructor, makes construction easier for 1D-4D vectors
29   - CUDA_CALLABLE vec(T x)
  29 + CUDA_CALLABLE vec(T rhs)
30 30 {
31   - v[0] = x;
  31 + for(int i=0; i<N; i++)
  32 + v[i] = rhs;
32 33 }
33 34 CUDA_CALLABLE vec(T x, T y)
34 35 {
... ... @@ -49,9 +50,10 @@ struct vec
49 50 v[3] = w;
50 51 }
51 52  
52   - CUDA_CALLABLE vec(const T(&data)[N])
53   - {
54   - memcpy(v, data, sizeof(T) * N);
  53 + //copy constructor
  54 + CUDA_CALLABLE vec( const vec<T, N>& other){
  55 + for(int i=0; i<N; i++)
  56 + v[i] = other.v[i];
55 57 }
56 58  
57 59 CUDA_CALLABLE T len() const
... ... @@ -135,12 +137,21 @@ struct vec
135 137 //arithmetic
136 138 CUDA_CALLABLE vec<T, N> operator+(vec<T, N> rhs) const
137 139 {
138   - vec<T, N> result;
  140 + vec<T, N> result;
139 141  
140   - for(int i=0; i<N; i++)
141   - result.v[i] = v[i] + rhs.v[i];
  142 + for(int i=0; i<N; i++)
  143 + result.v[i] = v[i] + rhs.v[i];
142 144  
143   - return result;
  145 + return result;
  146 + }
  147 + CUDA_CALLABLE vec<T, N> operator+(T rhs) const
  148 + {
  149 + vec<T, N> result;
  150 +
  151 + for(int i=0; i<N; i++)
  152 + result.v[i] = v[i] + rhs;
  153 +
  154 + return result;
144 155 }
145 156 CUDA_CALLABLE vec<T, N> operator-(vec<T, N> rhs) const
146 157 {
... ... @@ -174,15 +185,21 @@ struct vec
174 185 v[i] = v[i] * rhs;
175 186 return *this;
176 187 }
  188 + CUDA_CALLABLE vec<T, N> & operator=(T rhs){
  189 + for(int i=0; i<N; i++)
  190 + v[i] = rhs;
  191 + return *this;
  192 + }
  193 + //unary minus
  194 + CUDA_CALLABLE vec<T, N> operator-() const{
  195 + vec<T, N> r;
177 196  
178   - //conversion from a point
179   - /*CUDA_CALLABLE vector<T, N> & operator=(point<T, N> rhs)
180   - {
181   - for(int n=0; n<N; n++)
182   - v[n] = rhs.p[n];
  197 + //negate the vector
  198 + for(int i=0; i<N; i++)
  199 + r.v[i] = -v[i];
183 200  
184   - return *this;
185   - }*/
  201 + return r;
  202 + }
186 203  
187 204 CUDA_CALLABLE bool operator==(vec<T, N> rhs) const
188 205 {
... ... @@ -226,18 +243,7 @@ std::ostream&amp; operator&lt;&lt;(std::ostream&amp; os, rts::vec&lt;T, N&gt; v)
226 243 return os;
227 244 }
228 245  
229   -//arithmetic operators
230   -template <typename T, int N>
231   -CUDA_CALLABLE rts::vec<T, N> operator-(rts::vec<T, N> v)
232   -{
233   - rts::vec<T, N> r;
234   -
235   - //negate the vector
236   - for(int i=0; i<N; i++)
237   - r.v[i] = -v.v[i];
238 246  
239   - return r;
240   -}
241 247  
242 248 template <typename T, int N>
243 249 CUDA_CALLABLE rts::vec<T, N> operator*(T lhs, rts::vec<T, N> rhs)
... ...
optics/beam.h
  Show/Hide comments   View file @559d0fc
... ... @@ -16,10 +16,10 @@ public:
16 16  
17 17 private:
18 18  
19   - P na[2]; //numerical aperature of the focusing optics
  19 + P _na[2]; //numerical aperature of the focusing optics
20 20 vec<P> f; //focal point
21 21 function<P, P> apod; //apodization function
22   - unsigned int apod_res; //resolution of complex apodization filters
  22 + unsigned int apod_res; //resolution of apodization filter functions
23 23  
24 24 void apod_uniform()
25 25 {
... ... @@ -69,28 +69,44 @@ public:
69 69  
70 70 ///constructor: build a default beam (NA=1.0)
71 71 beam(
72   - vec<P> _k = rts::vec<P>(0, 0, TAU),
  72 + vec<P> k = rts::vec<P>(0, 0, rtsTAU),
73 73 vec<P> _E0 = rts::vec<P>(1, 0, 0),
74 74 beam_type _apod = Uniform)
75   - : planewave<P>(_k, _E0)
  75 + : planewave<P>(k, _E0)
76 76 {
77   - na[0] = (P)0.0;
78   - na[1] = (P)1.0;
  77 + _na[0] = (P)0.0;
  78 + _na[1] = (P)1.0;
79 79 f = vec<P>( (P)0, (P)0, (P)0 );
80 80 apod_res = 256; //set the default resolution for apodization filters
81 81 set_apod(_apod); //set the apodization function type
82 82 }
83 83  
  84 + beam<P> refract(rts::vec<P> kn) const{
  85 +
  86 + beam<P> new_beam;
  87 + new_beam._na[0] = _na[0];
  88 + new_beam._na[1] = _na[1];
  89 +
  90 +
  91 + rts::planewave<P> pw = planewave<P>::bend(kn);
  92 + //std::cout<<pw.str()<<std::endl;
  93 +
  94 + new_beam.k = pw.kvec();
  95 + new_beam.E0 = pw.E();
  96 +
  97 + return new_beam;
  98 + }
  99 +
84 100 ///Numerical Aperature functions
85   - void NA(P _na)
  101 + void NA(P na)
86 102 {
87   - na[0] = (P)0;
88   - na[1] = _na;
  103 + _na[0] = (P)0;
  104 + _na[1] = na;
89 105 }
90   - void NA(P _na0, P _na1)
  106 + void NA(P na0, P na1)
91 107 {
92   - na[0] = _na0;
93   - na[1] = _na1;
  108 + _na[0] = na0;
  109 + _na[1] = na1;
94 110 }
95 111  
96 112 /*string str() :
... ... @@ -116,7 +132,12 @@ public:
116 132 ///compute the rotation operator to transform (0, 0, 1) to k
117 133 P cos_angle = k_hat.dot(rts::vec<P>(0, 0, 1));
118 134 rts::matrix<P, 3> rotation;
119   - if(cos_angle != 1.0)
  135 +
  136 + //if the cosine of the angle is -1, the rotation is just a flip across the z axis
  137 + if(cos_angle == -1){
  138 + rotation(2, 2) = -1;
  139 + }
  140 + else if(cos_angle != 1.0)
120 141 {
121 142 rts::vec<P> r_axis = rts::vec<P>(0, 0, 1).cross(k_hat).norm(); //compute the axis of rotation
122 143 P angle = acos(cos_angle); //compute the angle of rotation
... ... @@ -127,8 +148,8 @@ public:
127 148  
128 149 //find the phi values associated with the cassegrain ring
129 150 P PHI[2];
130   - PHI[0] = (P)asin(na[0]);
131   - PHI[1] = (P)asin(na[1]);
  151 + PHI[0] = (P)asin(_na[0]);
  152 + PHI[1] = (P)asin(_na[1]);
132 153  
133 154 //calculate the z-axis cylinder coordinates associated with these angles
134 155 P Z[2];
... ... @@ -152,12 +173,28 @@ public:
152 173 vec<P> k_prime = rotation * cart; //create a sample vector
153 174  
154 175 //store a wave refracted along the given direction
155   - samples.push_back(beam::refract(k_prime) * apod(phi/PHI[1]));
  176 + //std::cout<<"k prime: "<<rotation<<std::endl;
  177 + samples.push_back(planewave<P>::refract(k_prime) * apod(phi/PHI[1]));
156 178 }
157 179  
158 180 return samples;
159 181 }
160 182  
  183 + std::string str()
  184 + {
  185 + std::stringstream ss;
  186 + ss<<"Beam:"<<std::endl;
  187 + ss<<" Central Plane Wave: "<<beam::E0<<" e^i ( "<<beam::k<<" . r )"<<std::endl;
  188 + if(_na[0] == 0)
  189 + ss<<" NA: "<<_na[1];
  190 + else
  191 + ss<<" NA: "<<_na[0]<<" -- "<<_na[1];
  192 +
  193 + return ss.str();
  194 + }
  195 +
  196 +
  197 +
161 198 };
162 199  
163 200 }
... ...
optics/efield.cuh
  Show/Hide comments   View file @559d0fc
... ... @@ -8,7 +8,14 @@
8 8 #include "../cuda/devices.h"
9 9 #include "../optics/beam.h"
10 10  
  11 +namespace rts{
  12 +template<typename T> class halfspace;
11 13  
  14 +template<typename T> class efield;
  15 +}
  16 +
  17 +template<typename T>
  18 +void operator<<(rts::efield<T> &ef, rts::halfspace<T> hs);
12 19  
13 20 namespace rts{
14 21  
... ... @@ -63,6 +70,27 @@ __global__ void gpu_efield_magnitude(complex&lt;T&gt;* X, complex&lt;T&gt;* Y, complex&lt;T&gt;* Z
63 70 }
64 71  
65 72 template<typename T>
  73 +__global__ void gpu_efield_real_magnitude(complex<T>* X, complex<T>* Y, complex<T>* Z, unsigned int r0, unsigned int r1, T* M)
  74 +{
  75 + int iu = blockIdx.x * blockDim.x + threadIdx.x;
  76 + int iv = blockIdx.y * blockDim.y + threadIdx.y;
  77 +
  78 + //make sure that the thread indices are in-bounds
  79 + if(iu >= r0 || iv >= r1) return;
  80 +
  81 + //compute the index into the field
  82 + int i = iv*r0 + iu;
  83 +
  84 + if(Y == NULL) //if this is a scalar simulation
  85 + M[i] = X[i].abs(); //compute the magnitude of the X component
  86 + else
  87 + {
  88 + M[i] = rts::vec<T>(X[i].real(), Y[i].real(), Z[i].real()).len();
  89 + //M[i] = Z[i].abs();
  90 + }
  91 +}
  92 +
  93 +template<typename T>
66 94 __global__ void gpu_efield_polarization(complex<T>* X, complex<T>* Y, complex<T>* Z,
67 95 unsigned int r0, unsigned int r1,
68 96 T* Px, T* Py, T* Pz)
... ... @@ -79,6 +107,7 @@ __global__ void gpu_efield_polarization(complex&lt;T&gt;* X, complex&lt;T&gt;* Y, complex&lt;T&gt;
79 107 //compute the field polarization
80 108 Px[i] = X[i].abs();
81 109 Py[i] = Y[i].abs();
  110 +
82 111 Pz[i] = Z[i].abs();
83 112 }
84 113  
... ... @@ -146,7 +175,6 @@ protected:
146 175  
147 176 if(vector)
148 177 {
149   - std::cout<<"vector:";
150 178 cudaMalloc(&Y, sizeof(rts::complex<P>) * R[0] * R[1]);
151 179 cudaMemset(Y, 0, sizeof(rts::complex<P>) * R[0] * R[1]);
152 180  
... ... @@ -280,6 +308,7 @@ public:
280 308 return *this; //return this object
281 309 }
282 310  
  311 + //assignment operator: build an electric field from a beam
283 312 efield<P> & operator= (const rts::beam<P> & rhs)
284 313 {
285 314 //get a vector of monte-carlo samples
... ... @@ -311,6 +340,25 @@ public:
311 340 return M;
312 341 }
313 342  
  343 + //return a sacalar field representing the field magnitude at an infinitely small point in time
  344 + realfield<P, 1, true> real_mag()
  345 + {
  346 + int maxThreads = rts::maxThreadsPerBlock(); //compute the optimal block size
  347 + int SQRT_BLOCK = (int)std::sqrt((float)maxThreads);
  348 +
  349 + //create a scalar field to store the result
  350 + realfield<P, 1, true> M(R[0], R[1]);
  351 +
  352 + //create one thread for each detector pixel
  353 + dim3 dimBlock(SQRT_BLOCK, SQRT_BLOCK);
  354 + dim3 dimGrid((R[0] + SQRT_BLOCK -1)/SQRT_BLOCK, (R[1] + SQRT_BLOCK - 1)/SQRT_BLOCK);
  355 +
  356 + //compute the magnitude and store it in a scalar field
  357 + gpu_efield_real_magnitude<float> <<<dimGrid, dimBlock>>> (X, Y, Z, R[0], R[1], M.ptr(0));
  358 +
  359 + return M;
  360 + }
  361 +
314 362 //return a vector field representing field polarization
315 363 realfield<P, 3, true> polarization()
316 364 {
... ... @@ -335,6 +383,8 @@ public:
335 383 return Pol; //return the vector field
336 384 }
337 385  
  386 + ////////FRIENDSHIP
  387 + friend void operator<< <>(rts::efield<P> &ef, rts::halfspace<P> hs);
338 388  
339 389  
340 390 };
... ...
optics/esphere.cuh
  Show/Hide comments   View file @559d0fc
... ... @@ -134,7 +134,7 @@ public:
134 134 calcAB(rhs.kmag(), Nl, A, B); //compute the scattering coefficients
135 135  
136 136 //determine important parameters for the scattering domain
137   - unsigned int sR = ceil(sqrt( (P)(pow(esphere::R[0],2) + pow(esphere::R[1],2))) );
  137 + unsigned int sR = ceil(sqrt( (P)(pow((P)esphere::R[0],2) + pow((P)esphere::R[1],2))) );
138 138 unsigned int thetaR = 256;
139 139  
140 140 /////////////////////continue scattering code here/////////////////////////
... ...
optics/halfspace.cuh 0 → 100644
  Show/Hide comments   View file @559d0fc
  1 +#ifndef RTS_HALFSPACE_CUH
  2 +#define RTS_HALFSPACE_CUH
  3 +
  4 +#include "../optics/halfspace.h"
  5 +#include "../optics/efield.cuh"
  6 +
  7 +namespace rts{
  8 +
  9 +template<typename T>
  10 +__global__ void gpu_halfspace_pw2ef(complex<T>* X, complex<T>* Y, complex<T>* Z, unsigned int r0, unsigned int r1,
  11 + plane<T> P, planewave<T> w, rect<T> q, bool at_surface = false)
  12 +{
  13 + int iu = blockIdx.x * blockDim.x + threadIdx.x;
  14 + int iv = blockIdx.y * blockDim.y + threadIdx.y;
  15 +
  16 + //make sure that the thread indices are in-bounds
  17 + if(iu >= r0 || iv >= r1) return;
  18 +
  19 + //compute the index into the field
  20 + int i = iv*r0 + iu;
  21 +
  22 + //get the current position
  23 + vec<T> p = q( (T)iu/(T)r0, (T)iv/(T)r1 );
  24 +
  25 + if(at_surface){
  26 + if(P.side(p) > 0)
  27 + return;
  28 + }
  29 + else{
  30 + if(P.side(p) >= 0)
  31 + return;
  32 + }
  33 +
  34 + //if the current position is on the wrong side of the plane
  35 +
  36 + //vec<T> r(p[0], p[1], p[2]);
  37 +
  38 + complex<T> x( 0.0f, w.kvec().dot(p) );
  39 + complex<T> phase( 0.0f, w.phase());
  40 +
  41 + if(Y == NULL) //if this is a scalar simulation
  42 + X[i] += w.E().len() * exp(x); //use the vector magnitude as the plane wave amplitude
  43 + else
  44 + {
  45 + //X[i] = Y[i] = Z[i] = 1;
  46 + X[i] += w.E()[0] * exp(x) * exp(phase);
  47 + Y[i] += w.E()[1] * exp(x) * exp(phase);
  48 + Z[i] += w.E()[2] * exp(x) * exp(phase);
  49 + }
  50 +}
  51 +}
  52 +
  53 +template<class T>
  54 +void operator<<(rts::efield<T> &ef, rts::halfspace<T> hs){
  55 +
  56 + //std::cout<<"---------RENDER HALFSPACE--------------"<<std::endl;
  57 + //output a parameter of the efield
  58 + //std::cout<<ef.pos<<std::endl;
  59 + //std::cout<<hs.S.str()<<std::endl;
  60 +
  61 + unsigned int SQRT_BLOCK = 16;
  62 + //create one thread for each detector pixel
  63 + dim3 dimBlock(SQRT_BLOCK, SQRT_BLOCK);
  64 + dim3 dimGrid((ef.R[0] + SQRT_BLOCK -1)/SQRT_BLOCK, (ef.R[1] + SQRT_BLOCK - 1)/SQRT_BLOCK);
  65 +
  66 +
  67 +
  68 + //render each plane wave
  69 +
  70 + //plane waves at the surface front
  71 + for(unsigned int w = 0; w < hs.w0.size(); w++){
  72 + //std::cout<<"w0 "<<w<<": "<<hs.w0[w].str()<<std::endl;
  73 + rts::gpu_halfspace_pw2ef<T> <<<dimGrid, dimBlock>>> (ef.X, ef.Y, ef.Z, ef.R[0], ef.R[1], hs.S, hs.w0[w], ef.pos);
  74 + }
  75 +
  76 + //plane waves at the surface back
  77 + for(unsigned int w = 0; w < hs.w1.size(); w++){
  78 + //std::cout<<"w1 "<<w<<": "<<hs.w1[w].str()<<std::endl;
  79 + rts::gpu_halfspace_pw2ef<T> <<<dimGrid, dimBlock>>> (ef.X, ef.Y, ef.Z, ef.R[0], ef.R[1], -hs.S, hs.w1[w], ef.pos, true);
  80 + }
  81 + //rts::gpu_halfspace_pw2ef<T> <<<dimGrid, dimBlock>>> (ef.X, ef.Y, ef.Z, ef.R[0], ef.R[1], hs.S, hs.pi, ef.pos);
  82 + //rts::gpu_halfspace_pw2ef<T> <<<dimGrid, dimBlock>>> (ef.X, ef.Y, ef.Z, ef.R[0], ef.R[1], hs.S, hs.pr, ef.pos);
  83 + //rts::gpu_halfspace_pw2ef<T> <<<dimGrid, dimBlock>>> (ef.X, ef.Y, ef.Z, ef.R[0], ef.R[1], -hs.S, hs.pt, ef.pos, true);
  84 +
  85 +
  86 + //return ef;
  87 +}
  88 +
  89 +
  90 +
  91 +#endif
0 92 \ No newline at end of file
... ...
optics/halfspace.h 0 → 100644
  Show/Hide comments   View file @559d0fc
  1 +#ifndef RTS_HALFSPACE_H
  2 +#define RTS_HALFSPACE_H
  3 +
  4 +namespace rts{
  5 +template<typename T> class halfspace;
  6 +
  7 +template<typename T> class efield;
  8 +}
  9 +
  10 +template<typename T>
  11 +void operator<<(rts::efield<T> &ef, rts::halfspace<T> hs);
  12 +
  13 +namespace rts{
  14 +
  15 +template<class T>
  16 +class halfspace
  17 +{
  18 +private:
  19 + rts::plane<T> S; //surface plane splitting the half space
  20 + rts::complex<T> n0; //refractive index at the front of the plane
  21 + rts::complex<T> n1; //refractive index at the back of the plane
  22 +
  23 + //lists of waves in front (pw0) and behind (pw1) the plane
  24 + std::vector< rts::planewave<T> > w0;
  25 + std::vector< rts::planewave<T> > w1;
  26 +
  27 + //rts::planewave<T> pi; //incident plane wave
  28 + //rts::planewave<T> pr; //plane wave reflected from the surface
  29 + //rts::planewave<T> pt; //plane wave transmitted through the surface
  30 +
  31 + void init(){
  32 + n0 = 1.0;
  33 + n1 = 1.5;
  34 + }
  35 +
  36 + //compute which side of the interface is hit by the incoming plane wave (0 = front, 1 = back)
  37 + bool facing(planewave<T> p){
  38 + if(p.kvec().dot(S.norm()) > 0)
  39 + return 1;
  40 + else
  41 + return 0;
  42 + }
  43 +
  44 + T calc_theta_i(vec<T> v){
  45 +
  46 + vec<T> v_hat = v.norm();
  47 +
  48 + //compute the cosine of the angle between k_hat and the plane normal
  49 + T cos_theta_i = v_hat.dot(S.norm());
  50 +
  51 + return acos(abs(cos_theta_i));
  52 + }
  53 +
  54 + T calc_theta_t(T ni_nt, T theta_i){
  55 +
  56 + T sin_theta_t = ni_nt * sin(theta_i);
  57 + return asin(sin_theta_t);
  58 + }
  59 +
  60 +
  61 +public:
  62 +
  63 + //constructors
  64 + halfspace(){
  65 + init();
  66 + }
  67 +
  68 + halfspace(T na, T nb){
  69 + init();
  70 + n0 = na;
  71 + n1 = nb;
  72 + }
  73 +
  74 + halfspace(T na, T nb, plane<T> p){
  75 + n0 = na;
  76 + n1 = nb;
  77 + S = p;
  78 + }
  79 +
  80 + //compute the transmitted and reflective waves given the incident (vacuum) plane wave p
  81 + void incident(rts::planewave<T> p){
  82 +
  83 + planewave<T> r, t;
  84 + p.scatter(S, n1.real()/n0.real(), r, t);
  85 +
  86 + //std::cout<<"i+r: "<<p.r()[0] + r.r()[0]<<std::endl;
  87 + //std::cout<<"t: "<<t.r()[0]<<std::endl;
  88 +
  89 + if(facing(p)){
  90 + w1.push_back(p);
  91 +
  92 + if(r.E().len() != 0)
  93 + w1.push_back(r);
  94 + if(t.E().len() != 0)
  95 + w0.push_back(t);
  96 + }
  97 + else{
  98 + w0.push_back(p);
  99 +
  100 + if(r.E().len() != 0)
  101 + w0.push_back(r);
  102 + if(t.E().len() != 0)
  103 + w1.push_back(t);
  104 + }
  105 + }
  106 +
  107 + void incident(rts::beam<T> b, unsigned int N = 10000){
  108 +
  109 + //generate a plane wave decomposition for the beam
  110 + std::vector< planewave<T> > pw_list = b.mc(N);
  111 +
  112 + //calculate the reflected and refracted waves for each incident wave
  113 + for(unsigned int w = 0; w < pw_list.size(); w++){
  114 + incident(pw_list[w]);
  115 + }
  116 + }
  117 +
  118 + std::string str(){
  119 + std::stringstream ss;
  120 + ss<<"Half Space-------------"<<std::endl;
  121 + ss<<"n0: "<<n0<<std::endl;
  122 + ss<<"n1: "<<n1<<std::endl;
  123 + ss<<S.str();
  124 +
  125 +
  126 + if(w0.size() != 0 || w1.size() != 0){
  127 + ss<<std::endl<<"Plane Waves:"<<std::endl;
  128 + for(unsigned int w = 0; w < w0.size(); w++){
  129 + ss<<"w0 = "<<w<<": "<<w0[w]<<std::endl;
  130 + }
  131 + for(unsigned int w = 0; w < w1.size(); w++){
  132 + ss<<"w1 = "<<w<<": "<<w1[w]<<std::endl;
  133 + }
  134 + }
  135 +
  136 + return ss.str();
  137 + }
  138 + ////////FRIENDSHIP
  139 + friend void operator<< <> (rts::efield<T> &ef, rts::halfspace<T> hs);
  140 +
  141 +};
  142 +
  143 +
  144 +
  145 +
  146 +}
  147 +
  148 +
  149 +#endif
0 150 \ No newline at end of file
... ...
optics/planewave.h
  Show/Hide comments   View file @559d0fc
... ... @@ -7,6 +7,8 @@
7 7 #include "../math/vector.h"
8 8 #include "../math/quaternion.h"
9 9 #include "../math/constants.h"
  10 +#include "../math/plane.h"
  11 +#include "rts/cuda/callable.h"
10 12  
11 13 /*Basic conversions used here (assuming a vacuum)
12 14 lambda =
... ... @@ -14,160 +16,299 @@
14 16  
15 17 namespace rts{
16 18  
17   -template<typename P>
  19 +template<typename T>
18 20 class planewave{
19 21  
  22 +protected:
  23 +
  24 + vec<T> k; //k = tau / lambda
  25 + vec<T> E0; //amplitude
  26 + T phi;
  27 +
  28 + planewave<T> bend(rts::vec<T> kn) const{
  29 +
  30 + vec<T> kn_hat = kn.norm(); //normalize the new k
  31 + vec<T> k_hat = k.norm(); //normalize the current k
  32 +
  33 + //std::cout<<"PLANE WAVE BENDING------------------"<<std::endl;
  34 + //std::cout<<"kn_hat: "<<kn_hat<<" k_hat: "<<k_hat<<std::endl;
  35 +
  36 + planewave<T> new_p; //create a new plane wave
  37 +
  38 + //if kn is equal to k or -k, handle the degenerate case
  39 + T k_dot_kn = k_hat.dot(kn_hat);
  40 +
  41 + //if k . n < 0, then the bend is a reflection
  42 + //flip k_hat
  43 + if(k_dot_kn < 0) k_hat = -k_hat;
  44 +
  45 + //std::cout<<"k dot kn: "<<k_dot_kn<<std::endl;
  46 +
  47 + //std::cout<<"k_dot_kn: "<<k_dot_kn<<std::endl;
  48 + if(k_dot_kn == -1){
  49 + new_p.k = -k;
  50 + new_p.E0 = E0;
  51 + return new_p;
  52 + }
  53 + else if(k_dot_kn == 1){
  54 + new_p.k = k;
  55 + new_p.E0 = E0;
  56 + return new_p;
  57 + }
  58 +
  59 + vec<T> r = k_hat.cross(kn_hat); //compute the rotation vector
  60 +
  61 + //std::cout<<"r: "<<r<<std::endl;
  62 +
  63 + T theta = asin(r.len()); //compute the angle of the rotation about r
  64 +
  65 +
  66 +
  67 + //deal with a zero vector (both k and kn point in the same direction)
  68 + //if(theta == (T)0)
  69 + //{
  70 + // new_p = *this;
  71 + // return new_p;
  72 + //}
  73 +
  74 + //create a quaternion to capture the rotation
  75 + quaternion<T> q;
  76 + q.CreateRotation(theta, r.norm());
  77 +
  78 + //apply the rotation to E0
  79 + vec<T> E0n = q.toMatrix3() * E0;
  80 +
  81 + new_p.k = kn_hat * kmag();
  82 + new_p.E0 = E0n;
  83 +
  84 + return new_p;
  85 + }
  86 +
20 87 public:
21   - vec<P> k; //k = tau / lambda
22   - vec<P> E0; //amplitude
23 88  
24 89  
25 90 ///constructor: create a plane wave propagating along z, polarized along x
26   - /*planewave(P lambda = (P)1)
  91 + /*planewave(T lambda = (T)1)
27 92 {
28   - k = rts::vec<P>(0, 0, 1) * (TAU/lambda);
29   - E0 = rts::vec<P>(1, 0, 0);
  93 + k = rts::vec<T>(0, 0, 1) * (TAU/lambda);
  94 + E0 = rts::vec<T>(1, 0, 0);
30 95 }*/
31   - ///constructor: create a plane wave propagating along _k, polarized along _E0, at frequency _omega
32   - planewave(vec<P> _k = rts::vec<P>(0, 0, TAU), vec<P> _E0 = rts::vec<P>(1, 0, 0))
  96 + ///constructor: create a plane wave propagating along k, polarized along _E0, at frequency _omega
  97 + CUDA_CALLABLE planewave(vec<T> kvec = rts::vec<T>(0, 0, rtsTAU),
  98 + vec<T> E = rts::vec<T>(1, 0, 0),
  99 + T phase = 0)
33 100 {
34   - k = _k;
35   - vec<P> k_hat = _k.norm();
36   -
37   - vec<P> s = (k.cross(_E0)).norm(); //compute an orthogonal side vector
38   - vec<P> E0_hat = (s.cross(k)).norm(); //compute a normalized E0 direction vector
39   - E0 = E0_hat * E0_hat.dot(_E0); //compute the projection of _E0 onto E0_hat
  101 + phi = phase;
  102 + k = kvec;
  103 + vec<T> k_hat = k.norm();
  104 +
  105 + if(E.len() == 0) //if the plane wave has an amplitude of 0
  106 + E0 = vec<T>(0); //just return it
  107 + else{
  108 + vec<T> s = (k.cross(E)).norm(); //compute an orthogonal side vector
  109 + vec<T> E_hat = (s.cross(k)).norm(); //compute a normalized E0 direction vector
  110 + E0 = E_hat * E_hat.dot(E); //compute the projection of _E0 onto E0_hat
  111 + }
40 112 }
41 113  
42 114 ///multiplication operator: scale E0
43   - planewave<P> & operator* (const P & rhs)
  115 + CUDA_CALLABLE planewave<T> & operator* (const T & rhs)
44 116 {
45 117  
46 118 E0 = E0 * rhs;
47 119 return *this;
48 120 }
49 121  
50   - P lambda() const
  122 + CUDA_CALLABLE T lambda() const
51 123 {
52   - return TAU / k.len();
  124 + return rtsTAU / k.len();
53 125 }
54 126  
55   - P kmag() const
  127 + CUDA_CALLABLE T kmag() const
56 128 {
57 129 return k.len();
58 130 }
59 131  
  132 + CUDA_CALLABLE vec<T> E(){
  133 + return E0;
  134 + }
  135 +
  136 + CUDA_CALLABLE vec<T> kvec(){
  137 + return k;
  138 + }
  139 +
  140 + CUDA_CALLABLE T phase(){
  141 + return phi;
  142 + }
  143 +
  144 + CUDA_CALLABLE void phase(T p){
  145 + phi = p;
  146 + }
  147 +
  148 + CUDA_CALLABLE vec< complex<T> > pos(vec<T> p = vec<T>(0, 0, 0)){
  149 + vec< complex<T> > result;
  150 +
  151 + T kdp = k.dot(p);
  152 + complex<T> x = complex<T>(0, kdp);
  153 + complex<T> expx = exp(x);
  154 +
  155 + result[0] = E0[0] * expx;
  156 + result[1] = E0[1] * expx;
  157 + result[2] = E0[2] * expx;
  158 +
  159 + return result;
  160 + }
  161 +
  162 + //scales k based on a transition from material ni to material nt
  163 + CUDA_CALLABLE planewave<T> n(T ni, T nt){
  164 + return planewave<T>(k * (nt / ni), E0);
  165 + }
60 166  
61   - planewave<P> refract(rts::vec<P> kn) const
  167 + CUDA_CALLABLE planewave<T> refract(rts::vec<T> kn) const
62 168 {
63   - vec<P> kn_hat = kn.norm(); //normalize the new k
64   - vec<P> k_hat = k.norm(); //normalize the current k
  169 + return bend(kn);
  170 + }
65 171  
66   - vec<P> r = k_hat.cross(kn_hat); //compute the rotation vector
  172 + void scatter(rts::plane<T> P, T nr, planewave<T> &r, planewave<T> &t){
67 173  
68   - P theta = asin(r.len()); //compute the angle of the rotation about r
  174 + int facing = P.face(k); //determine which direction the plane wave is coming in
69 175  
70   - planewave<P> new_p; //create a new plane wave
  176 + //if(facing == 0) //if the wave is tangent to the plane, return an identical wave
  177 + // return *this;
  178 + //else
  179 + if(facing == -1){ //if the wave hits the back of the plane, invert the plane and nr
  180 + P = P.flip(); //flip the plane
  181 + nr = 1/nr; //invert the refractive index (now nr = n0/n1)
  182 + }
71 183  
72   - //deal with a zero vector (both k and kn point in the same direction)
73   - if(theta == (P)0)
74   - {
75   - new_p = *this;
76   - return new_p;
  184 + //use Snell's Law to calculate the transmitted angle
  185 + T cos_theta_i = k.norm().dot(-P.norm()); //compute the cosine of theta_i
  186 + T theta_i = acos(cos_theta_i); //compute theta_i
  187 + T sin_theta_t = (1/nr) * sin(theta_i); //compute the sine of theta_t using Snell's law
  188 + T theta_t = asin(sin_theta_t); //compute the cosine of theta_t
  189 +
  190 + bool tir = false; //flag for total internal reflection
  191 + if(theta_t != theta_t){
  192 + tir = true;
  193 + theta_t = rtsPI / (T)2;
77 194 }
78 195  
79   - //create a quaternion to capture the rotation
80   - quaternion<P> q;
81   - q.CreateRotation(theta, r.norm());
  196 + //handle the degenerate case where theta_i is 0 (the plane wave hits head-on)
  197 + if(theta_i == 0){
  198 + T rp = (1 - nr) / (1 + nr); //compute the Fresnel coefficients
  199 + T tp = 2 / (1 + nr);
  200 + vec<T> kr = -k;
  201 + vec<T> kt = k * nr; //set the k vectors for theta_i = 0
  202 + vec<T> Er = E0 * rp; //compute the E vectors
  203 + vec<T> Et = E0 * tp;
  204 + T phase_t = P.p().dot(k - kt); //compute the phase offset
  205 + T phase_r = P.p().dot(k - kr);
  206 + //std::cout<<"Degeneracy: Head-On"<<std::endl;
  207 + //std::cout<<"rs: "<<rp<<" rp: "<<rp<<" ts: "<<tp<<" tp: "<<tp<<std::endl;
  208 + //std::cout<<"phase r: "<<phase_r<<" phase t: "<<phase_t<<std::endl;
  209 +
  210 + //create the plane waves
  211 + r = planewave<T>(kr, Er, phase_r);
  212 + t = planewave<T>(kt, Et, phase_t);
  213 +
  214 + //std::cout<<"i + r: "<<pos()[0] + r.pos()[0]<<pos()[1] + r.pos()[1]<<pos()[2] + r.pos()[2]<<std::endl;
  215 + //std::cout<<"t: "<<t.pos()[0]<<t.pos()[1]<<t.pos()[2]<<std::endl;
  216 + //std::cout<<"--------------------------------"<<std::endl;
  217 + return;
  218 + }
82 219  
83   - //apply the rotation to E0
84   - vec<P> E0n = q.toMatrix3() * E0;
85 220  
86   - new_p.k = kn_hat * k.len();
87   - new_p.E0 = E0n;
  221 + //compute the Fresnel coefficients
  222 + T rp, rs, tp, ts;
  223 + rp = tan(theta_t - theta_i) / tan(theta_t + theta_i);
  224 + rs = sin(theta_t - theta_i) / sin(theta_t + theta_i);
  225 +
  226 + if(tir){
  227 + tp = ts = 0;
  228 + }
  229 + else{
  230 + tp = ( 2 * sin(theta_t) * cos(theta_i) ) / ( sin(theta_t + theta_i) * cos(theta_t - theta_i) );
  231 + ts = ( 2 * sin(theta_t) * cos(theta_i) ) / sin(theta_t + theta_i);
  232 + }
88 233  
89   - return new_p;
  234 + //compute the coordinate space for the plane of incidence
  235 + vec<T> z_hat = -P.norm();
  236 + vec<T> y_hat = P.parallel(k).norm();
  237 + vec<T> x_hat = y_hat.cross(z_hat).norm();
  238 +
  239 + //compute the k vectors for r and t
  240 + vec<T> kr, kt;
  241 + kr = ( y_hat * sin(theta_i) - z_hat * cos(theta_i) ) * kmag();
  242 + kt = ( y_hat * sin(theta_t) + z_hat * cos(theta_t) ) * kmag() * nr;
  243 +
  244 + //compute the magnitude of the p- and s-polarized components of the incident E vector
  245 + T Ei_s = E0.dot(x_hat);
  246 + int sgn = (0 < E0.dot(y_hat)) - (E0.dot(y_hat) < 0);
  247 + T Ei_p = sgn * ( E0 - x_hat * Ei_s ).len();
  248 + //T Ei_p = ( E0 - x_hat * Ei_s ).len();
  249 + //compute the magnitude of the p- and s-polarized components of the reflected E vector
  250 + T Er_s = Ei_s * rs;
  251 + T Er_p = Ei_p * rp;
  252 + //compute the magnitude of the p- and s-polarized components of the transmitted E vector
  253 + T Et_s = Ei_s * ts;
  254 + T Et_p = Ei_p * tp;
  255 +
  256 + //std::cout<<"E0: "<<E0<<std::endl;
  257 + //std::cout<<"E0 dot y_hat: "<<E0.dot(y_hat)<<std::endl;
  258 + //std::cout<<"theta i: "<<theta_i<<" theta t: "<<theta_t<<std::endl;
  259 + //std::cout<<"x_hat: "<<x_hat<<" y_hat: "<<y_hat<<" z_hat: "<<z_hat<<std::endl;
  260 + //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;
  261 + //std::cout<<"rs: "<<rs<<" rp: "<<rp<<" ts: "<<ts<<" tp: "<<tp<<std::endl;
  262 +
  263 +
  264 + //compute the reflected E vector
  265 + vec<T> Er = (y_hat * cos(theta_i) + z_hat * sin(theta_i)) * Er_p + x_hat * Er_s;
  266 + //compute the transmitted E vector
  267 + vec<T> Et = (y_hat * cos(theta_t) - z_hat * sin(theta_t)) * Et_p + x_hat * Et_s;
  268 +
  269 + T phase_t = P.p().dot(k - kt);
  270 + T phase_r = P.p().dot(k - kr);
  271 +
  272 + //std::cout<<"phase r: "<<phase_r<<" phase t: "<<phase_t<<std::endl;
  273 +
  274 + //std::cout<<"phase: "<<phase<<std::endl;
  275 +
  276 + //create the plane waves
  277 + r.k = kr;
  278 + r.E0 = Er;
  279 + r.phi = phase_r;
  280 +
  281 + //t = bend(kt);
  282 + //t.k = t.k * nr;
  283 +
  284 + t.k = kt;
  285 + t.E0 = Et;
  286 + t.phi = phase_t;
  287 + //std::cout<<"i: "<<str()<<std::endl;
  288 + //std::cout<<"r: "<<r.str()<<std::endl;
  289 + //std::cout<<"t: "<<t.str()<<std::endl;
  290 +
  291 + //std::cout<<"i + r: "<<pos()[0] + r.pos()[0]<<pos()[1] + r.pos()[1]<<pos()[2] + r.pos()[2]<<std::endl;
  292 + //std::cout<<"t: "<<t.pos()[0]<<t.pos()[1]<<t.pos()[2]<<std::endl;
  293 + //std::cout<<"--------------------------------"<<std::endl;
90 294  
91   - /*vec<P> kn_hat = kn.norm(); //normalize new_k
92   - vec<P> k_hat = k.norm();
93   -
94   - //compute the side vector (around which we will be rotating)
95   - vec<P> E0_hat = E0.norm();
96   - rts::vec<P> s = k_hat.cross(E0_hat);
97   -
98   - //compute the projection of k' onto the k-E plane
99   - rts::vec<P> s_prime = s * (kn_hat.dot(s));
100   -
101   - //compute the angle between
102   - vec<P> kp_hat = (kn_hat - s_prime).norm();
103   - P theta = acos(k_hat.dot( kp_hat ));
104   - if(kn_hat.dot(E0_hat) < 0)
105   - theta = -theta;
106   -
107   - //rotate E0 around s by theta
108   - quaternion<P> q;
109   - q.CreateRotation(theta, s);
110   - rts::vec<P> E0_prime = q.toMatrix3() * E0;
111   -
112   - //create the refracted plane wave
113   - planewave<P> new_p;
114   - new_p.E0 = E0_prime;
115   - new_p.k = kn_hat * k.len();
116   -
117   - return new_p;*/
118   -
119   - /*vec<P> kn_hat = kn.norm(); //normalize kn
120   - vec<P> k_hat = k.norm();
121   - vec<P> E0_hat = E0.norm();
122   -
123   - vec<P> B = k_hat.cross(E0_hat);
124   - planewave<P> newp;
125   - newp.k = kn_hat * k.len();
126   - newp.E0 = B.cross(kn_hat).norm();
127   - std::cout<<newp.k.norm().dot(newp.E0.norm())<<std::endl;
128   - return newp;*/
129   -/*
130   - //compute the side vector (around which we will be rotating)
131   - rts::vec<P> s_hat = k_hat.cross(E0_hat);
132   - //std::cout<<s.len()<<std::endl;
133   -
134   - //project kn_hat into the k-E0 plane
135   - rts::vec<P> sp = s_hat * (kn_hat.dot(s_hat)); //project k_new onto s
136   - rts::vec<P> kp = (kn_hat - sp); //correct k_new so it lies on the E0-k plane
137   - rts::vec<P> kp_hat = kp.norm();
138   -
139   - //compute the angle and direction between k_prime and k
140   - P theta = acos(k_hat.dot(kp_hat));
141   - if(kp_hat.dot(E0_hat) < 0)
142   - theta = -theta;
143   -
144   - //rotate E0 around s by theta
145   - quaternion<P> q;
146   - q.CreateRotation(theta, s_hat);
147   - rts::vec<P> E0n = q.toMatrix3() * E0;
148   - rts::vec<P> E0n_hat = E0n.norm();
149   -
150   - //std::cout<<s_hat.dot(kp_hat)<<" "<<s_hat.dot(E0n_hat)<<" "<<s_hat.dot(E0_hat)<<" "<<s_hat.dot(k_hat)<<" "<<
151   - // E0_hat.dot(k_hat)<<" "<<k_hat.dot(kp_hat)<<" "<<E0_hat.dot(E0n_hat)<<" "<<E0n_hat.dot(kp_hat)<<std::endl;
152   -
153   - //create the refracted plane wave
154   - //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;
155   -
156   - //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;
157   -
158   - planewave<P> new_p(kn_hat * k.len(), E0n); //create the plane wave
159   - std::cout<<"|E0n|: "<<new_p.E0.len()<<" E0n*kn: "<<(new_p.E0.norm()).dot(new_p.k.norm())<<std::endl;
160   -
161   - return new_p;*/
162 295 }
163 296  
164 297 std::string str()
165 298 {
166 299 std::stringstream ss;
167   - ss<<E0<<" e^i ( "<<k<<" . r )";
  300 + ss<<"Plane Wave:"<<std::endl;
  301 + ss<<" "<<E0<<" e^i ( "<<k<<" . r )";
168 302 return ss.str();
169 303 }
170 304 };
171 305 }
172 306  
  307 +template <typename T>
  308 +std::ostream& operator<<(std::ostream& os, rts::planewave<T> p)
  309 +{
  310 + os<<p.str();
  311 + return os;
  312 +}
  313 +
173 314 #endif
174 315 \ No newline at end of file
... ...