codebase / stimlib

Commit ecfd14dfcbaf0a9cdcfb883ccd359bc2e6de5ab4

Authored by David Mayerich
1 parent 559d0fcb

implemented D field in halfspace

Showing 11 changed files with 485 additions and 280 deletions   Show diff stats
Inline Side-by-side
math/complex.h
  Show/Hide comments   View file @ecfd14d
... ... @@ -67,10 +67,17 @@ struct complex
67 67 return result;
68 68 }
69 69  
  70 + //returns the complex signum (-1, 0, 1)
  71 + CUDA_CALLABLE int sgn(){
  72 + if(r > 0) return 1;
  73 + else if(r < 0) return -1;
  74 + else return (0 < i - i < 0);
  75 + }
  76 +
70 77 //ARITHMETIC OPERATORS--------------------
71 78  
72 79 //binary + operator (returns the result of adding two complex values)
73   - CUDA_CALLABLE complex<T> operator+ (const complex<T> rhs)
  80 + CUDA_CALLABLE complex<T> operator+ (const complex<T> rhs) const
74 81 {
75 82 complex<T> result;
76 83 result.r = r + rhs.r;
... ... @@ -78,7 +85,7 @@ struct complex
78 85 return result;
79 86 }
80 87  
81   - CUDA_CALLABLE complex<T> operator+ (const T rhs)
  88 + CUDA_CALLABLE complex<T> operator+ (const T rhs) const
82 89 {
83 90 complex<T> result;
84 91 result.r = r + rhs;
... ... @@ -87,7 +94,7 @@ struct complex
87 94 }
88 95  
89 96 //binary - operator (returns the result of adding two complex values)
90   - CUDA_CALLABLE complex<T> operator- (const complex<T> rhs)
  97 + CUDA_CALLABLE complex<T> operator- (const complex<T> rhs) const
91 98 {
92 99 complex<T> result;
93 100 result.r = r - rhs.r;
... ... @@ -105,7 +112,7 @@ struct complex
105 112 }
106 113  
107 114 //binary MULTIPLICATION operators (returns the result of multiplying complex values)
108   - CUDA_CALLABLE complex<T> operator* (const complex<T> rhs)
  115 + CUDA_CALLABLE complex<T> operator* (const complex<T> rhs) const
109 116 {
110 117 complex<T> result;
111 118 result.r = r * rhs.r - i * rhs.i;
... ... @@ -118,7 +125,7 @@ struct complex
118 125 }
119 126  
120 127 //binary DIVISION operators (returns the result of dividing complex values)
121   - CUDA_CALLABLE complex<T> operator/ (const complex<T> rhs)
  128 + CUDA_CALLABLE complex<T> operator/ (const complex<T> rhs) const
122 129 {
123 130 complex<T> result;
124 131 T denom = rhs.r * rhs.r + rhs.i * rhs.i;
... ... @@ -257,7 +264,7 @@ struct complex
257 264 return result;
258 265 }
259 266  
260   - std::string toStr()
  267 + std::string str()
261 268 {
262 269 std::stringstream ss;
263 270 ss<<"("<<r<<","<<i<<")";
... ... @@ -280,6 +287,13 @@ struct complex
280 287 return false;
281 288 }
282 289  
  290 + CUDA_CALLABLE bool operator!=(T rhs)
  291 + {
  292 + if(r != rhs || i != 0)
  293 + return true;
  294 + return false;
  295 + }
  296 +
283 297 //CASTING operators
284 298 template < typename otherT >
285 299 operator complex<otherT>()
... ... @@ -451,7 +465,7 @@ CUDA_CALLABLE rts::complex&lt;A&gt; cos(const rts::complex&lt;A&gt; x)
451 465 template<class A>
452 466 std::ostream& operator<<(std::ostream& os, rts::complex<A> x)
453 467 {
454   - os<<x.toStr();
  468 + os<<x.str();
455 469 return os;
456 470 }
457 471  
... ...
math/matrix.h
  Show/Hide comments   View file @ecfd14d
... ... @@ -40,9 +40,10 @@ struct matrix
40 40 return *this;
41 41 }
42 42  
43   - CUDA_CALLABLE vec<T, N> operator*(vec<T, N> rhs)
  43 + template<typename Y>
  44 + CUDA_CALLABLE vec<Y, N> operator*(vec<Y, N> rhs)
44 45 {
45   - vec<T, N> result;
  46 + vec<Y, N> result;
46 47  
47 48 for(int r=0; r<N; r++)
48 49 for(int c=0; c<N; c++)
... ...
math/plane.h
  Show/Hide comments   View file @ecfd14d
... ... @@ -58,6 +58,13 @@ public:
58 58 N = N.norm();
59 59 }
60 60  
  61 + template< typename U >
  62 + CUDA_CALLABLE operator plane<U, D>(){
  63 +
  64 + plane<U, D> result(N, P);
  65 + return result;
  66 + }
  67 +
61 68 CUDA_CALLABLE vec<T, D> norm(){
62 69 return N;
63 70 }
... ... @@ -105,6 +112,11 @@ public:
105 112 return v - perpendicular(v);
106 113 }
107 114  
  115 + CUDA_CALLABLE void decompose(vec<T, D> v, vec<T, D>& para, vec<T, D>& perp){
  116 + perp = N * v.dot(N);
  117 + para = v - perp;
  118 + }
  119 +
108 120 //get both the parallel and perpendicular components of a vector v w.r.t. the plane
109 121 CUDA_CALLABLE void project(vec<T, D> v, vec<T, D> &v_par, vec<T, D> &v_perp){
110 122  
... ... @@ -129,6 +141,16 @@ public:
129 141  
130 142 }
131 143  
  144 + CUDA_CALLABLE rts::plane<T, D> operator-()
  145 + {
  146 + rts::plane<T, D> p = *this;
  147 +
  148 + //negate the normal vector
  149 + p.N = -p.N;
  150 +
  151 + return p;
  152 + }
  153 +
132 154 //output a string
133 155 std::string str(){
134 156 std::stringstream ss;
... ... @@ -138,7 +160,7 @@ public:
138 160 }
139 161  
140 162 ///////Friendship
141   - friend CUDA_CALLABLE rts::plane<T, D> operator- <> (rts::plane<T, D> v);
  163 + //friend CUDA_CALLABLE rts::plane<T, D> operator- <> (rts::plane<T, D> v);
142 164  
143 165  
144 166  
... ... @@ -149,16 +171,8 @@ public:
149 171 //arithmetic operators
150 172  
151 173 //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;
  174 +//template <typename T, int D>
159 175  
160   - return p;
161   -}
162 176  
163 177  
164 178  
... ...
math/realfield.cuh
  Show/Hide comments   View file @ecfd14d
... ... @@ -35,6 +35,24 @@ __global__ void gpu_gaussian(T* dest, unsigned int r0, unsigned int r1, T mean,
35 35  
36 36 namespace rts{
37 37  
  38 +//multiply R = X * Y
  39 +template<typename T>
  40 +__global__ void gpu_field_multiply(T* R, T* X, T* Y, unsigned int r0, unsigned int r1){
  41 +
  42 + int iu = blockIdx.x * blockDim.x + threadIdx.x;
  43 + int iv = blockIdx.y * blockDim.y + threadIdx.y;
  44 +
  45 + //make sure that the thread indices are in-bounds
  46 + if(iu >= r0 || iv >= r1) return;
  47 +
  48 + //compute the index into the field
  49 + int i = iv*r0 + iu;
  50 +
  51 + //calculate and store the result
  52 + R[i] = X[i] * Y[i];
  53 +
  54 +}
  55 +
38 56 template<typename P, unsigned int N = 1, bool positive = false>
39 57 class realfield{
40 58  
... ... @@ -107,7 +125,7 @@ public:
107 125 //std::cout<<"realfield DESTRUCTOR"<<std::endl;
108 126 }
109 127  
110   - P* ptr(unsigned int n)
  128 + P* ptr(unsigned int n = 0)
111 129 {
112 130 if(n < N)
113 131 return X[n];
... ... @@ -212,6 +230,27 @@ public:
212 230 return *this;
213 231 }
214 232  
  233 + //multiply two fields (element-wise multiplication)
  234 + realfield<P, N, positive> operator* (const realfield & rhs){
  235 +
  236 + int maxThreads = rts::maxThreadsPerBlock(); //compute the optimal block size
  237 + int SQRT_BLOCK = (int)std::sqrt((float)maxThreads);
  238 +
  239 + //create one thread for each detector pixel
  240 + dim3 dimBlock(SQRT_BLOCK, SQRT_BLOCK);
  241 + dim3 dimGrid((R[0] + SQRT_BLOCK -1)/SQRT_BLOCK, (R[1] + SQRT_BLOCK - 1)/SQRT_BLOCK);
  242 +
  243 + //create a scalar field to store the result
  244 + realfield<P, N, positive> result(R[0], R[1]);
  245 +
  246 + for(int n=0; n<N; n++)
  247 + rts::gpu_field_multiply <<<dimGrid, dimBlock>>> (result.X[n], X[n], rhs.X[n], R[0], R[1]);
  248 +
  249 + return result;
  250 +
  251 +
  252 + }
  253 +
215 254 ///copy constructor
216 255 realfield(const realfield &rhs)
217 256 {
... ...
math/vector.h
  Show/Hide comments   View file @ecfd14d
... ... @@ -56,6 +56,19 @@ struct vec
56 56 v[i] = other.v[i];
57 57 }
58 58  
  59 +
  60 + template< typename U >
  61 + CUDA_CALLABLE operator vec<U, N>(){
  62 + vec<U, N> result;
  63 + for(int i=0; i<N; i++)
  64 + result.v[i] = v[i];
  65 +
  66 + return result;
  67 + }
  68 +
  69 + //template<class U>
  70 + //friend vec<U, N>::operator vec<T, N>();
  71 +
59 72 CUDA_CALLABLE T len() const
60 73 {
61 74 //compute and return the vector length
... ... @@ -64,7 +77,7 @@ struct vec
64 77 {
65 78 sum_sq += v[i] * v[i];
66 79 }
67   - return std::sqrt(sum_sq);
  80 + return sqrt(sum_sq);
68 81  
69 82 }
70 83  
... ... @@ -190,6 +203,13 @@ struct vec
190 203 v[i] = rhs;
191 204 return *this;
192 205 }
  206 +
  207 + template<typename Y>
  208 + CUDA_CALLABLE vec<T, N> & operator=(vec<Y, N> rhs){
  209 + for(int i=0; i<N; i++)
  210 + v[i] = rhs.v[i];
  211 + return *this;
  212 + }
193 213 //unary minus
194 214 CUDA_CALLABLE vec<T, N> operator-() const{
195 215 vec<T, N> r;
... ... @@ -209,7 +229,7 @@ struct vec
209 229 return false;
210 230 }
211 231  
212   - std::string toStr() const
  232 + std::string str() const
213 233 {
214 234 std::stringstream ss;
215 235  
... ... @@ -239,7 +259,7 @@ struct vec
239 259 template <typename T, int N>
240 260 std::ostream& operator<<(std::ostream& os, rts::vec<T, N> v)
241 261 {
242   - os<<v.toStr();
  262 + os<<v.str();
243 263 return os;
244 264 }
245 265  
... ...
optics/beam.h
  Show/Hide comments   View file @ecfd14d
... ... @@ -184,7 +184,8 @@ public:
184 184 {
185 185 std::stringstream ss;
186 186 ss<<"Beam:"<<std::endl;
187   - ss<<" Central Plane Wave: "<<beam::E0<<" e^i ( "<<beam::k<<" . r )"<<std::endl;
  187 + //ss<<" Central Plane Wave: "<<beam::E0<<" e^i ( "<<beam::k<<" . r )"<<std::endl;
  188 + ss<<" Central Plane Wave: "<<beam::k<<std::endl;
188 189 if(_na[0] == 0)
189 190 ss<<" NA: "<<_na[1];
190 191 else
... ...
optics/efield.cuh
  Show/Hide comments   View file @ecfd14d
... ... @@ -7,15 +7,8 @@
7 7 #include "../optics/planewave.h"
8 8 #include "../cuda/devices.h"
9 9 #include "../optics/beam.h"
  10 +#include "../math/rect.h"
10 11  
11   -namespace rts{
12   -template<typename T> class halfspace;
13   -
14   -template<typename T> class efield;
15   -}
16   -
17   -template<typename T>
18   -void operator<<(rts::efield<T> &ef, rts::halfspace<T> hs);
19 12  
20 13 namespace rts{
21 14  
... ... @@ -107,8 +100,8 @@ __global__ void gpu_efield_polarization(complex&lt;T&gt;* X, complex&lt;T&gt;* Y, complex&lt;T&gt;
107 100 //compute the field polarization
108 101 Px[i] = X[i].abs();
109 102 Py[i] = Y[i].abs();
110   -
111 103 Pz[i] = Z[i].abs();
  104 +
112 105 }
113 106  
114 107 /* This function computes the sum of two complex fields and stores the result in *dest
... ... @@ -132,7 +125,7 @@ __global__ void gpu_efield_sum(complex&lt;T&gt;* dest, complex&lt;T&gt;* src, unsigned int r
132 125 /* This class implements a discrete representation of an electromagnetic field
133 126 in 2D. The majority of this representation is done on the GPU.
134 127 */
135   -template<typename P>
  128 +template<typename T>
136 129 class efield
137 130 {
138 131 protected:
... ... @@ -140,17 +133,17 @@ protected:
140 133 bool vector;
141 134  
142 135 //gpu pointer to the field data
143   - rts::complex<P>* X;
144   - rts::complex<P>* Y;
145   - rts::complex<P>* Z;
  136 + rts::complex<T>* X;
  137 + rts::complex<T>* Y;
  138 + rts::complex<T>* Z;
146 139  
147 140 //resolution of the discrete field
148 141 unsigned int R[2];
149 142  
150 143 //shape of the 2D field
151   - rect<P> pos;
  144 + rect<T> pos;
152 145  
153   - void from_planewave(planewave<P> p)
  146 + void from_planewave(planewave<T> p)
154 147 {
155 148 unsigned int SQRT_BLOCK = 16;
156 149 //create one thread for each detector pixel
... ... @@ -158,7 +151,7 @@ protected:
158 151 dim3 dimGrid((R[0] + SQRT_BLOCK -1)/SQRT_BLOCK, (R[1] + SQRT_BLOCK - 1)/SQRT_BLOCK);
159 152  
160 153  
161   - gpu_planewave2efield<P> <<<dimGrid, dimBlock>>> (X, Y, Z, R[0], R[1], p, pos);
  154 + gpu_planewave2efield<T> <<<dimGrid, dimBlock>>> (X, Y, Z, R[0], R[1], p, pos);
162 155 }
163 156  
164 157 void init(unsigned int res0, unsigned int res1, bool _vector)
... ... @@ -170,16 +163,16 @@ protected:
170 163 R[1] = res1;
171 164  
172 165 //allocate memory on the gpu
173   - cudaMalloc(&X, sizeof(rts::complex<P>) * R[0] * R[1]);
174   - cudaMemset(X, 0, sizeof(rts::complex<P>) * R[0] * R[1]);
  166 + cudaMalloc(&X, sizeof(rts::complex<T>) * R[0] * R[1]);
  167 + cudaMemset(X, 0, sizeof(rts::complex<T>) * R[0] * R[1]);
175 168  
176 169 if(vector)
177 170 {
178   - cudaMalloc(&Y, sizeof(rts::complex<P>) * R[0] * R[1]);
179   - cudaMemset(Y, 0, sizeof(rts::complex<P>) * R[0] * R[1]);
  171 + cudaMalloc(&Y, sizeof(rts::complex<T>) * R[0] * R[1]);
  172 + cudaMemset(Y, 0, sizeof(rts::complex<T>) * R[0] * R[1]);
180 173  
181   - cudaMalloc(&Z, sizeof(rts::complex<P>) * R[0] * R[1]);
182   - cudaMemset(Z, 0, sizeof(rts::complex<P>) * R[0] * R[1]);
  174 + cudaMalloc(&Z, sizeof(rts::complex<T>) * R[0] * R[1]);
  175 + cudaMemset(Z, 0, sizeof(rts::complex<T>) * R[0] * R[1]);
183 176 }
184 177 }
185 178  
... ... @@ -190,14 +183,14 @@ protected:
190 183 if(Z != NULL) cudaFree(Z);
191 184 }
192 185  
193   - void shallowcpy(const rts::efield<P> & src)
  186 + void shallowcpy(const rts::efield<T> & src)
194 187 {
195 188 vector = src.vector;
196 189 R[0] = src.R[0];
197 190 R[1] = src.R[1];
198 191 }
199 192  
200   - void deepcpy(const rts::efield<P> & src)
  193 + void deepcpy(const rts::efield<T> & src)
201 194 {
202 195 //perform a shallow copy
203 196 shallowcpy(src);
... ... @@ -205,18 +198,18 @@ protected:
205 198 //allocate memory on the gpu
206 199 if(src.X != NULL)
207 200 {
208   - cudaMalloc(&X, sizeof(rts::complex<P>) * R[0] * R[1]);
209   - cudaMemcpy(X, src.X, sizeof(rts::complex<P>) * R[0] * R[1], cudaMemcpyDeviceToDevice);
  201 + cudaMalloc(&X, sizeof(rts::complex<T>) * R[0] * R[1]);
  202 + cudaMemcpy(X, src.X, sizeof(rts::complex<T>) * R[0] * R[1], cudaMemcpyDeviceToDevice);
210 203 }
211 204 if(src.Y != NULL)
212 205 {
213   - cudaMalloc(&Y, sizeof(rts::complex<P>) * R[0] * R[1]);
214   - cudaMemcpy(Y, src.Y, sizeof(rts::complex<P>) * R[0] * R[1], cudaMemcpyDeviceToDevice);
  206 + cudaMalloc(&Y, sizeof(rts::complex<T>) * R[0] * R[1]);
  207 + cudaMemcpy(Y, src.Y, sizeof(rts::complex<T>) * R[0] * R[1], cudaMemcpyDeviceToDevice);
215 208 }
216 209 if(src.Z != NULL)
217 210 {
218   - cudaMalloc(&Z, sizeof(rts::complex<P>) * R[0] * R[1]);
219   - cudaMemcpy(Z, src.Z, sizeof(rts::complex<P>) * R[0] * R[1], cudaMemcpyDeviceToDevice);
  211 + cudaMalloc(&Z, sizeof(rts::complex<T>) * R[0] * R[1]);
  212 + cudaMemcpy(Z, src.Z, sizeof(rts::complex<T>) * R[0] * R[1], cudaMemcpyDeviceToDevice);
220 213 }
221 214 }
222 215  
... ... @@ -224,7 +217,7 @@ public:
224 217 efield(unsigned int res0, unsigned int res1, bool _vector = true)
225 218 {
226 219 init(res0, res1, _vector);
227   - //pos = rts::rect<P>(rts::vec<P>(-10, 0, -10), rts::vec<P>(-10, 0, 10), rts::vec<P>(10, 0, 10));
  220 + //pos = rts::rect<T>(rts::vec<T>(-10, 0, -10), rts::vec<T>(-10, 0, 10), rts::vec<T>(10, 0, 10));
228 221 }
229 222  
230 223 //destructor
... ... @@ -236,20 +229,40 @@ public:
236 229 ///Clear the field - set all points to zero
237 230 void clear()
238 231 {
239   - cudaMemset(X, 0, sizeof(rts::complex<P>) * R[0] * R[1]);
  232 + cudaMemset(X, 0, sizeof(rts::complex<T>) * R[0] * R[1]);
240 233  
241 234 if(vector)
242 235 {
243   - cudaMemset(Y, 0, sizeof(rts::complex<P>) * R[0] * R[1]);
244   - cudaMemset(Z, 0, sizeof(rts::complex<P>) * R[0] * R[1]);
  236 + cudaMemset(Y, 0, sizeof(rts::complex<T>) * R[0] * R[1]);
  237 + cudaMemset(Z, 0, sizeof(rts::complex<T>) * R[0] * R[1]);
245 238 }
246 239 }
247 240  
248   - void position(rect<P> _p)
  241 + void position(rect<T> _p)
249 242 {
250 243 pos = _p;
251 244 }
252 245  
  246 + //access functions
  247 + complex<T>* x(){
  248 + return X;
  249 + }
  250 + complex<T>* y(){
  251 + return Y;
  252 + }
  253 + complex<T>* z(){
  254 + return Z;
  255 + }
  256 + unsigned int Ru(){
  257 + return R[0];
  258 + }
  259 + unsigned int Rv(){
  260 + return R[1];
  261 + }
  262 + rect<T> p(){
  263 + return pos;
  264 + }
  265 +
253 266 std::string str()
254 267 {
255 268 stringstream ss;
... ... @@ -262,7 +275,7 @@ public:
262 275 }
263 276  
264 277 //assignment operator: assignment from another electric field
265   - efield<P> & operator= (const efield<P> & rhs)
  278 + efield<T> & operator= (const efield<T> & rhs)
266 279 {
267 280 destroy(); //destroy any previous information about this field
268 281 deepcpy(rhs); //create a deep copy
... ... @@ -270,7 +283,7 @@ public:
270 283 }
271 284  
272 285 //assignment operator: build an electric field from a plane wave
273   - efield<P> & operator= (const planewave<P> & rhs)
  286 + efield<T> & operator= (const planewave<T> & rhs)
274 287 {
275 288  
276 289 clear(); //clear any previous field data
... ... @@ -279,7 +292,7 @@ public:
279 292 }
280 293  
281 294 //assignment operator: add an existing electric field
282   - efield<P> & operator+= (const efield<P> & rhs)
  295 + efield<T> & operator+= (const efield<T> & rhs)
283 296 {
284 297 //if this field and the source field represent the same regions in space
285 298 if(R[0] == rhs.R[0] && R[1] == rhs.R[1] && pos == rhs.pos)
... ... @@ -309,10 +322,10 @@ public:
309 322 }
310 323  
311 324 //assignment operator: build an electric field from a beam
312   - efield<P> & operator= (const rts::beam<P> & rhs)
  325 + efield<T> & operator= (const rts::beam<T> & rhs)
313 326 {
314 327 //get a vector of monte-carlo samples
315   - std::vector< rts::planewave<P> > p_list = rhs.mc();
  328 + std::vector< rts::planewave<T> > p_list = rhs.mc();
316 329  
317 330 clear(); //clear any previous field data
318 331 for(unsigned int i = 0; i < p_list.size(); i++)
... ... @@ -322,13 +335,13 @@ public:
322 335  
323 336  
324 337 //return a scalar field representing field magnitude
325   - realfield<P, 1, true> mag()
  338 + realfield<T, 1, true> mag()
326 339 {
327 340 int maxThreads = rts::maxThreadsPerBlock(); //compute the optimal block size
328 341 int SQRT_BLOCK = (int)std::sqrt((float)maxThreads);
329 342  
330 343 //create a scalar field to store the result
331   - realfield<P, 1, true> M(R[0], R[1]);
  344 + realfield<T, 1, true> M(R[0], R[1]);
332 345  
333 346 //create one thread for each detector pixel
334 347 dim3 dimBlock(SQRT_BLOCK, SQRT_BLOCK);
... ... @@ -341,13 +354,13 @@ public:
341 354 }
342 355  
343 356 //return a sacalar field representing the field magnitude at an infinitely small point in time
344   - realfield<P, 1, true> real_mag()
  357 + realfield<T, 1, true> real_mag()
345 358 {
346 359 int maxThreads = rts::maxThreadsPerBlock(); //compute the optimal block size
347 360 int SQRT_BLOCK = (int)std::sqrt((float)maxThreads);
348 361  
349 362 //create a scalar field to store the result
350   - realfield<P, 1, true> M(R[0], R[1]);
  363 + realfield<T, 1, true> M(R[0], R[1]);
351 364  
352 365 //create one thread for each detector pixel
353 366 dim3 dimBlock(SQRT_BLOCK, SQRT_BLOCK);
... ... @@ -360,7 +373,7 @@ public:
360 373 }
361 374  
362 375 //return a vector field representing field polarization
363   - realfield<P, 3, true> polarization()
  376 + realfield<T, 3, true> polarization()
364 377 {
365 378 if(!vector)
366 379 {
... ... @@ -375,7 +388,7 @@ public:
375 388 dim3 dimGrid((R[0] + SQRT_BLOCK -1)/SQRT_BLOCK, (R[1] + SQRT_BLOCK - 1)/SQRT_BLOCK);
376 389  
377 390  
378   - realfield<P, 3, true> Pol(R[0], R[1]); //create a vector field to store the result
  391 + realfield<T, 3, true> Pol(R[0], R[1]); //create a vector field to store the result
379 392  
380 393 //compute the polarization and store it in the vector field
381 394 gpu_efield_polarization<float> <<<dimGrid, dimBlock>>> (X, Y, Z, R[0], R[1], Pol.ptr(0), Pol.ptr(1), Pol.ptr(2));
... ... @@ -384,7 +397,7 @@ public:
384 397 }
385 398  
386 399 ////////FRIENDSHIP
387   - friend void operator<< <>(rts::efield<P> &ef, rts::halfspace<P> hs);
  400 + //friend void operator<< <T>(rts::efield<T> &ef, rts::halfspace<T> hs);
388 401  
389 402  
390 403 };
... ...
optics/esphere.cuh
  Show/Hide comments   View file @ecfd14d
... ... @@ -135,7 +135,7 @@ public:
135 135  
136 136 //determine important parameters for the scattering domain
137 137 unsigned int sR = ceil(sqrt( (P)(pow((P)esphere::R[0],2) + pow((P)esphere::R[1],2))) );
138   - unsigned int thetaR = 256;
  138 + //unsigned int thetaR = 256;
139 139  
140 140 /////////////////////continue scattering code here/////////////////////////
141 141  
... ...
optics/halfspace.cuh
  Show/Hide comments   View file @ecfd14d
1   -#ifndef RTS_HALFSPACE_CUH
2   -#define RTS_HALFSPACE_CUH
  1 +#ifndef RTS_HALFSPACE_H
  2 +#define RTS_HALFSPACE_H
  3 +
  4 +#include "../math/plane.h"
3 5  
4   -#include "../optics/halfspace.h"
5   -#include "../optics/efield.cuh"
6 6  
7 7 namespace rts{
8 8  
  9 +//GPU kernel to compute the electric field
9 10 template<typename T>
10 11 __global__ void gpu_halfspace_pw2ef(complex<T>* X, complex<T>* Y, complex<T>* Z, unsigned int r0, unsigned int r1,
11 12 plane<T> P, planewave<T> w, rect<T> q, bool at_surface = false)
... ... @@ -36,56 +37,303 @@ __global__ void gpu_halfspace_pw2ef(complex&lt;T&gt;* X, complex&lt;T&gt;* Y, complex&lt;T&gt;* Z,
36 37 //vec<T> r(p[0], p[1], p[2]);
37 38  
38 39 complex<T> x( 0.0f, w.kvec().dot(p) );
39   - complex<T> phase( 0.0f, w.phase());
  40 + //complex<T> phase( 0.0f, w.phase());
40 41  
41 42 if(Y == NULL) //if this is a scalar simulation
42 43 X[i] += w.E().len() * exp(x); //use the vector magnitude as the plane wave amplitude
43 44 else
44 45 {
45 46 //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);
  47 + X[i] += w.E()[0] * exp(x);// * exp(phase);
  48 + Y[i] += w.E()[1] * exp(x);// * exp(phase);
  49 + Z[i] += w.E()[2] * exp(x);// * exp(phase);
  50 + }
  51 +}
  52 +
  53 +//GPU kernel to compute the electric displacement field
  54 +template<typename T>
  55 +__global__ void gpu_halfspace_pw2df(complex<T>* X, complex<T>* Y, complex<T>* Z, unsigned int r0, unsigned int r1,
  56 + plane<T> P, planewave<T> w, rect<T> q, T n, bool at_surface = false)
  57 +{
  58 + int iu = blockIdx.x * blockDim.x + threadIdx.x;
  59 + int iv = blockIdx.y * blockDim.y + threadIdx.y;
  60 +
  61 + //make sure that the thread indices are in-bounds
  62 + if(iu >= r0 || iv >= r1) return;
  63 +
  64 + //compute the index into the field
  65 + int i = iv*r0 + iu;
  66 +
  67 + //get the current position
  68 + vec<T> p = q( (T)iu/(T)r0, (T)iv/(T)r1 );
  69 +
  70 + if(at_surface){
  71 + if(P.side(p) > 0)
  72 + return;
  73 + }
  74 + else{
  75 + if(P.side(p) >= 0)
  76 + return;
  77 + }
  78 +
  79 + //if the current position is on the wrong side of the plane
  80 +
  81 + //vec<T> r(p[0], p[1], p[2]);
  82 +
  83 + complex<T> x( 0.0f, w.kvec().dot(p) );
  84 + //complex<T> phase( 0.0f, w.phase());
  85 +
  86 + //vec< complex<T> > testE(1, 0, 0);
  87 +
  88 +
  89 +
  90 + if(Y == NULL) //if this is a scalar simulation
  91 + X[i] += w.E().len() * exp(x); //use the vector magnitude as the plane wave amplitude
  92 + else
  93 + {
  94 + plane< complex<T> > cplane = plane< complex<T>, 3>(P);
  95 + vec< complex<T> > E_para;// = cplane.parallel(w.E());
  96 + vec< complex<T> > E_perp;// = cplane.perpendicular(w.E()) * (n*n);
  97 + cplane.decompose(w.E(), E_para, E_perp);
  98 + T epsilon = n*n;
  99 +
  100 + X[i] += (E_para[0] + E_perp[0] * epsilon) * exp(x);
  101 + Y[i] += (E_para[1] + E_perp[1] * epsilon) * exp(x);
  102 + Z[i] += (E_para[2] + E_perp[2] * epsilon) * exp(x);
49 103 }
50 104 }
  105 +
  106 +//computes a scalar field containing the refractive index of the half-space at each point
  107 +template<typename T>
  108 +__global__ void gpu_halfspace_n(T* n, unsigned int r0, unsigned int r1, rect<T> q, plane<T> P, T n0, T n1){
  109 +
  110 + int iu = blockIdx.x * blockDim.x + threadIdx.x;
  111 + int iv = blockIdx.y * blockDim.y + threadIdx.y;
  112 +
  113 + //make sure that the thread indices are in-bounds
  114 + if(iu >= r0 || iv >= r1) return;
  115 +
  116 + //compute the index into the field
  117 + int i = iv*r0 + iu;
  118 +
  119 + //get the current position
  120 + vec<T> p = q( (T)iu/(T)r0, (T)iv/(T)r1 );
  121 +
  122 + //set the appropriate refractive index
  123 + if(P.side(p) < 0) n[i] = n0;
  124 + else n[i] = n1;
51 125 }
52 126  
53 127 template<class T>
54   -void operator<<(rts::efield<T> &ef, rts::halfspace<T> hs){
  128 +class halfspace
  129 +{
  130 +private:
  131 + rts::plane<T> S; //surface plane splitting the half space
  132 + rts::complex<T> n0; //refractive index at the front of the plane
  133 + rts::complex<T> n1; //refractive index at the back of the plane
  134 +
  135 + //lists of waves in front (pw0) and behind (pw1) the plane
  136 + std::vector< rts::planewave<T> > w0;
  137 + std::vector< rts::planewave<T> > w1;
55 138  
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;
  139 + //rts::planewave<T> pi; //incident plane wave
  140 + //rts::planewave<T> pr; //plane wave reflected from the surface
  141 + //rts::planewave<T> pt; //plane wave transmitted through the surface
60 142  
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);
  143 + void init(){
  144 + n0 = 1.0;
  145 + n1 = 1.5;
  146 + }
65 147  
  148 + //compute which side of the interface is hit by the incoming plane wave (0 = front, 1 = back)
  149 + bool facing(planewave<T> p){
  150 + if(p.kvec().dot(S.norm()) > 0)
  151 + return 1;
  152 + else
  153 + return 0;
  154 + }
66 155  
  156 + T calc_theta_i(vec<T> v){
67 157  
68   - //render each plane wave
  158 + vec<T> v_hat = v.norm();
69 159  
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);
  160 + //compute the cosine of the angle between k_hat and the plane normal
  161 + T cos_theta_i = v_hat.dot(S.norm());
  162 +
  163 + return acos(abs(cos_theta_i));
74 164 }
75 165  
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);
  166 + T calc_theta_t(T ni_nt, T theta_i){
  167 +
  168 + T sin_theta_t = ni_nt * sin(theta_i);
  169 + return asin(sin_theta_t);
80 170 }
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 171  
85 172  
86   - //return ef;
87   -}
  173 +public:
  174 +
  175 + //constructors
  176 + halfspace(){
  177 + init();
  178 + }
  179 +
  180 + halfspace(T na, T nb){
  181 + init();
  182 + n0 = na;
  183 + n1 = nb;
  184 + }
  185 +
  186 + halfspace(T na, T nb, plane<T> p){
  187 + n0 = na;
  188 + n1 = nb;
  189 + S = p;
  190 + }
  191 +
  192 + //compute the transmitted and reflective waves given the incident (vacuum) plane wave p
  193 + void incident(rts::planewave<T> p){
  194 +
  195 + planewave<T> r, t;
  196 + p.scatter(S, n1.real()/n0.real(), r, t);
  197 +
  198 + //std::cout<<"i+r: "<<p.r()[0] + r.r()[0]<<std::endl;
  199 + //std::cout<<"t: "<<t.r()[0]<<std::endl;
  200 +
  201 + if(facing(p)){
  202 + w1.push_back(p);
  203 +
  204 + if(r.E().len() != 0)
  205 + w1.push_back(r);
  206 + if(t.E().len() != 0)
  207 + w0.push_back(t);
  208 + }
  209 + else{
  210 + w0.push_back(p);
  211 +
  212 + if(r.E().len() != 0)
  213 + w0.push_back(r);
  214 + if(t.E().len() != 0)
  215 + w1.push_back(t);
  216 + }
  217 + }
  218 +
  219 + void incident(rts::beam<T> b, unsigned int N = 10000){
  220 +
  221 + //generate a plane wave decomposition for the beam
  222 + std::vector< planewave<T> > pw_list = b.mc(N);
  223 +
  224 + //calculate the reflected and refracted waves for each incident wave
  225 + for(unsigned int w = 0; w < pw_list.size(); w++){
  226 + incident(pw_list[w]);
  227 + }
  228 + }
  229 +
  230 + //return the electric field at the specified resolution and position
  231 + rts::efield<T> E(unsigned int r0, unsigned int r1, rect<T> R){
  232 +
  233 + int maxThreads = rts::maxThreadsPerBlock(); //compute the optimal block size
  234 + int SQRT_BLOCK = (int)std::sqrt((float)maxThreads);
  235 +
  236 + //create one thread for each detector pixel
  237 + dim3 dimBlock(SQRT_BLOCK, SQRT_BLOCK);
  238 + dim3 dimGrid((r0 + SQRT_BLOCK -1)/SQRT_BLOCK, (r1 + SQRT_BLOCK - 1)/SQRT_BLOCK);
  239 +
  240 + //create an electric field
  241 + rts::efield<T> ef(r0, r1);
  242 + ef.position(R);
88 243  
  244 + //render each plane wave
  245 +
  246 + //plane waves at the surface front
  247 + for(unsigned int w = 0; w < w0.size(); w++){
  248 + //std::cout<<"w0 "<<w<<": "<<hs.w0[w].str()<<std::endl;
  249 + rts::gpu_halfspace_pw2ef<T> <<<dimGrid, dimBlock>>> (ef.x(), ef.y(), ef.z(), r0, r1, S, w0[w], ef.p());
  250 + }
  251 +
  252 + //plane waves at the surface back
  253 + for(unsigned int w = 0; w < w1.size(); w++){
  254 + //std::cout<<"w1 "<<w<<": "<<hs.w1[w].str()<<std::endl;
  255 + rts::gpu_halfspace_pw2ef<T> <<<dimGrid, dimBlock>>> (ef.x(), ef.y(), ef.z(), r0, r1, -S, w1[w], ef.p(), true);
  256 + }
  257 +
  258 + return ef;
  259 + }
  260 +
  261 + //return the electric displacement at the specified resolution and position
  262 + rts::efield<T> D(unsigned int r0, unsigned int r1, rect<T> R){
  263 +
  264 + int maxThreads = rts::maxThreadsPerBlock(); //compute the optimal block size
  265 + int SQRT_BLOCK = (int)std::sqrt((float)maxThreads);
  266 +
  267 + //create one thread for each detector pixel
  268 + dim3 dimBlock(SQRT_BLOCK, SQRT_BLOCK);
  269 + dim3 dimGrid((r0 + SQRT_BLOCK -1)/SQRT_BLOCK, (r1 + SQRT_BLOCK - 1)/SQRT_BLOCK);
  270 +
  271 + //create a complex vector field
  272 + rts::efield<T> df(r0, r1);
  273 + df.position(R);
  274 +
  275 + //render each plane wave
  276 +
  277 + //plane waves at the surface front
  278 + for(unsigned int w = 0; w < w0.size(); w++){
  279 + //std::cout<<"w0 "<<w<<": "<<hs.w0[w].str()<<std::endl;
  280 + rts::gpu_halfspace_pw2df<T> <<<dimGrid, dimBlock>>> (df.x(), df.y(), df.z(), r0, r1, S, w0[w], df.p(), n0.real());
  281 + }
  282 +
  283 + //plane waves at the surface back
  284 + for(unsigned int w = 0; w < w1.size(); w++){
  285 + //std::cout<<"w1 "<<w<<": "<<hs.w1[w].str()<<std::endl;
  286 + rts::gpu_halfspace_pw2df<T> <<<dimGrid, dimBlock>>> (df.x(), df.y(), df.z(), r0, r1, -S, w1[w], df.p(), n1.real(), true);
  287 + }
  288 +
  289 + return df;
  290 + }
  291 +
  292 + realfield<T, 1, true> nfield(unsigned int Ru, unsigned int Rv, rect<T> p){
  293 +
  294 + int maxThreads = rts::maxThreadsPerBlock(); //compute the optimal block size
  295 + int SQRT_BLOCK = (int)std::sqrt((float)maxThreads);
  296 +
  297 + //create one thread for each detector pixel
  298 + dim3 dimBlock(SQRT_BLOCK, SQRT_BLOCK);
  299 + dim3 dimGrid((Ru + SQRT_BLOCK -1)/SQRT_BLOCK, (Rv + SQRT_BLOCK - 1)/SQRT_BLOCK);
  300 +
  301 + realfield<T, 1, true> n(Ru, Rv); //create a scalar field to store the result of n
  302 +
  303 + rts::gpu_halfspace_n<T> <<<dimGrid, dimBlock>>> (n.ptr(), Ru, Rv, p, S, n0.real(), n1.real());
  304 +
  305 + return n;
  306 + }
  307 +
  308 + std::string str(){
  309 + std::stringstream ss;
  310 + ss<<"Half Space-------------"<<std::endl;
  311 + ss<<"n0: "<<n0<<std::endl;
  312 + ss<<"n1: "<<n1<<std::endl;
  313 + ss<<S.str();
  314 +
  315 +
  316 + if(w0.size() != 0 || w1.size() != 0){
  317 + ss<<std::endl<<"Plane Waves:"<<std::endl;
  318 + for(unsigned int w = 0; w < w0.size(); w++){
  319 + ss<<"w0 = "<<w<<": "<<w0[w]<<std::endl;
  320 + }
  321 + for(unsigned int w = 0; w < w1.size(); w++){
  322 + ss<<"w1 = "<<w<<": "<<w1[w]<<std::endl;
  323 + }
  324 + }
  325 +
  326 + return ss.str();
  327 + }
  328 + ////////FRIENDSHIP
  329 + //friend void operator<< <> (rts::efield<T> &ef, rts::halfspace<T> hs);
  330 +
  331 +};
  332 +
  333 +
  334 +
  335 +
  336 +}
89 337  
90 338  
91 339 #endif
92 340 \ No newline at end of file
... ...
optics/halfspace.h deleted
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
150 0 \ No newline at end of file
optics/planewave.h
  Show/Hide comments   View file @ecfd14d
... ... @@ -22,10 +22,10 @@ class planewave{
22 22 protected:
23 23  
24 24 vec<T> k; //k = tau / lambda
25   - vec<T> E0; //amplitude
26   - T phi;
  25 + vec< complex<T> > E0; //amplitude
  26 + //T phi;
27 27  
28   - planewave<T> bend(rts::vec<T> kn) const{
  28 + CUDA_CALLABLE planewave<T> bend(rts::vec<T> kn) const{
29 29  
30 30 vec<T> kn_hat = kn.norm(); //normalize the new k
31 31 vec<T> k_hat = k.norm(); //normalize the current k
... ... @@ -76,7 +76,7 @@ protected:
76 76 q.CreateRotation(theta, r.norm());
77 77  
78 78 //apply the rotation to E0
79   - vec<T> E0n = q.toMatrix3() * E0;
  79 + vec< complex<T> > E0n = q.toMatrix3() * E0;
80 80  
81 81 new_p.k = kn_hat * kmag();
82 82 new_p.E0 = E0n;
... ... @@ -95,20 +95,22 @@ public:
95 95 }*/
96 96 ///constructor: create a plane wave propagating along k, polarized along _E0, at frequency _omega
97 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)
  98 + vec< complex<T> > E = rts::vec<T>(1, 0, 0), T phase = 0)
100 99 {
101   - phi = phase;
  100 + //phi = phase;
  101 +
102 102 k = kvec;
103   - vec<T> k_hat = k.norm();
  103 + vec< complex<T> > k_hat = k.norm();
104 104  
105 105 if(E.len() == 0) //if the plane wave has an amplitude of 0
106 106 E0 = vec<T>(0); //just return it
107 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
  108 + vec< complex<T> > s = (k_hat.cross(E)).norm(); //compute an orthogonal side vector
  109 + vec< complex<T> > E_hat = (s.cross(k)).norm(); //compute a normalized E0 direction vector
110 110 E0 = E_hat * E_hat.dot(E); //compute the projection of _E0 onto E0_hat
111 111 }
  112 +
  113 + E0 = E0 * exp( complex<T>(0, phase) );
112 114 }
113 115  
114 116 ///multiplication operator: scale E0
... ... @@ -129,7 +131,7 @@ public:
129 131 return k.len();
130 132 }
131 133  
132   - CUDA_CALLABLE vec<T> E(){
  134 + CUDA_CALLABLE vec< complex<T> > E(){
133 135 return E0;
134 136 }
135 137  
... ... @@ -137,13 +139,13 @@ public:
137 139 return k;
138 140 }
139 141  
140   - CUDA_CALLABLE T phase(){
  142 + /*CUDA_CALLABLE T phase(){
141 143 return phi;
142 144 }
143 145  
144 146 CUDA_CALLABLE void phase(T p){
145 147 phi = p;
146   - }
  148 + }*/
147 149  
148 150 CUDA_CALLABLE vec< complex<T> > pos(vec<T> p = vec<T>(0, 0, 0)){
149 151 vec< complex<T> > result;
... ... @@ -199,8 +201,8 @@ public:
199 201 T tp = 2 / (1 + nr);
200 202 vec<T> kr = -k;
201 203 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 + vec< complex<T> > Er = E0 * rp; //compute the E vectors
  205 + vec< complex<T> > Et = E0 * tp;
204 206 T phase_t = P.p().dot(k - kt); //compute the phase offset
205 207 T phase_r = P.p().dot(k - kr);
206 208 //std::cout<<"Degeneracy: Head-On"<<std::endl;
... ... @@ -242,16 +244,18 @@ public:
242 244 kt = ( y_hat * sin(theta_t) + z_hat * cos(theta_t) ) * kmag() * nr;
243 245  
244 246 //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();
  247 + complex<T> Ei_s = E0.dot(x_hat);
  248 + //int sgn = (0 < E0.dot(y_hat)) - (E0.dot(y_hat) < 0);
  249 + int sgn = E0.dot(y_hat).sgn();
  250 + vec< complex<T> > cx_hat = x_hat;
  251 + complex<T> Ei_p = ( E0 - cx_hat * Ei_s ).len() * sgn;
248 252 //T Ei_p = ( E0 - x_hat * Ei_s ).len();
249 253 //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;
  254 + complex<T> Er_s = Ei_s * rs;
  255 + complex<T> Er_p = Ei_p * rp;
252 256 //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;
  257 + complex<T> Et_s = Ei_s * ts;
  258 + complex<T> Et_p = Ei_p * tp;
255 259  
256 260 //std::cout<<"E0: "<<E0<<std::endl;
257 261 //std::cout<<"E0 dot y_hat: "<<E0.dot(y_hat)<<std::endl;
... ... @@ -262,9 +266,9 @@ public:
262 266  
263 267  
264 268 //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;
  269 + vec< complex<T> > Er = vec< complex<T> >(y_hat * cos(theta_i) + z_hat * sin(theta_i)) * Er_p + cx_hat * Er_s;
266 270 //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;
  271 + vec< complex<T> > Et = vec< complex<T> >(y_hat * cos(theta_t) - z_hat * sin(theta_t)) * Et_p + cx_hat * Et_s;
268 272  
269 273 T phase_t = P.p().dot(k - kt);
270 274 T phase_r = P.p().dot(k - kr);
... ... @@ -275,15 +279,15 @@ public:
275 279  
276 280 //create the plane waves
277 281 r.k = kr;
278   - r.E0 = Er;
279   - r.phi = phase_r;
  282 + r.E0 = Er * exp( complex<T>(0, phase_r) );
  283 + //r.phi = phase_r;
280 284  
281 285 //t = bend(kt);
282 286 //t.k = t.k * nr;
283 287  
284 288 t.k = kt;
285   - t.E0 = Et;
286   - t.phi = phase_t;
  289 + t.E0 = Et * exp( complex<T>(0, phase_t) );
  290 + //t.phi = phase_t;
287 291 //std::cout<<"i: "<<str()<<std::endl;
288 292 //std::cout<<"r: "<<r.str()<<std::endl;
289 293 //std::cout<<"t: "<<t.str()<<std::endl;
... ...