From 559d0fcbc33712d4de94f257b6b2718ab9f5c729 Mon Sep 17 00:00:00 2001
From: David Mayerich <mayerich@illinois.edu>
Date: Fri, 25 Jul 2014 16:32:36 -0500
Subject: [PATCH] wave interactions with a planar surface

---
 math/bessel.h        |  84 ++++++++++++++++++++++++++++++++++++++++++------------------------------------------
 math/complex.h       |  16 +++++++++-------
 math/constants.h     |   7 ++-----
 math/plane.h         | 165 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 math/realfield.cuh   |  10 +++++-----
 math/vector.h        |  60 +++++++++++++++++++++++++++++++++---------------------------
 optics/beam.h        |  69 +++++++++++++++++++++++++++++++++++++++++++++++++++++----------------
 optics/efield.cuh    |  52 +++++++++++++++++++++++++++++++++++++++++++++++++++-
 optics/esphere.cuh   |   2 +-
 optics/halfspace.cuh |  91 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 optics/halfspace.h   | 149 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 optics/planewave.h   | 359 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-------------------------------------------------------------------------------------------------------------
 12 files changed, 851 insertions(+), 213 deletions(-)
 create mode 100644 math/plane.h
 create mode 100644 optics/halfspace.cuh
 create mode 100644 optics/halfspace.h

diff --git a/math/bessel.h b/math/bessel.h
index 88bf4bf..9d12b8b 100644
--- a/math/bessel.h
+++ b/math/bessel.h
@@ -759,28 +759,28 @@ int bessjyv(P v,P x,P &vm,P *jv,P *yv,
     }
     vm = n + v0;
     return 0;
-}
-
-template<typename P>
+}
+
+template<typename P>
 int bessjyv_sph(int v, P z, P &vm, P* cjv,
-    P* cyv, P* cjvp, P* cyvp)
-{
-    //first, compute the bessel functions of fractional order
-    bessjyv(v + 0.5, z, vm, cjv, cyv, cjvp, cyvp);
-
-    //iterate through each and scale
-    for(int n = 0; n<=v; n++)
-    {
-
-        cjv[n] = cjv[n] * sqrt(PI/(z * 2.0));
-        cyv[n] = cyv[n] * sqrt(PI/(z * 2.0));
-
-        cjvp[n] = -1.0 / (z * 2.0) * cjv[n] + cjvp[n] * sqrt(PI / (z * 2.0));
-        cyvp[n] = -1.0 / (z * 2.0) * cyv[n] + cyvp[n] * sqrt(PI / (z * 2.0));
-    }
-
-	return 0;
-
+    P* cyv, P* cjvp, P* cyvp)
+{
+    //first, compute the bessel functions of fractional order
+    bessjyv(v + 0.5, z, vm, cjv, cyv, cjvp, cyvp);
+
+    //iterate through each and scale
+    for(int n = 0; n<=v; n++)
+    {
+
+        cjv[n] = cjv[n] * sqrt(rtsPI/(z * 2.0));
+        cyv[n] = cyv[n] * sqrt(rtsPI/(z * 2.0));
+
+        cjvp[n] = -1.0 / (z * 2.0) * cjv[n] + cjvp[n] * sqrt(rtsPI / (z * 2.0));
+        cyvp[n] = -1.0 / (z * 2.0) * cyv[n] + cyvp[n] * sqrt(rtsPI / (z * 2.0));
+    }
+
+	return 0;
+
 }
 
 template<typename P>
@@ -1485,28 +1485,28 @@ int cbessjyva(P v,complex<P> z,P &vm,complex<P>*cjv,
     }
     vm = n+v0;
     return 0;
-}
-
-template<typename P>
+}
+
+template<typename P>
 int cbessjyva_sph(int v,complex<P> z,P &vm,complex<P>*cjv,
-    complex<P>*cyv,complex<P>*cjvp,complex<P>*cyvp)
-{
-    //first, compute the bessel functions of fractional order
-    cbessjyva<P>(v + 0.5, z, vm, cjv, cyv, cjvp, cyvp);
-
-    //iterate through each and scale
-    for(int n = 0; n<=v; n++)
-    {
-
-        cjv[n] = cjv[n] * sqrt(PI/(z * 2.0));
-        cyv[n] = cyv[n] * sqrt(PI/(z * 2.0));
-
-        cjvp[n] = -1.0 / (z * 2.0) * cjv[n] + cjvp[n] * sqrt(PI / (z * 2.0));
-        cyvp[n] = -1.0 / (z * 2.0) * cyv[n] + cyvp[n] * sqrt(PI / (z * 2.0));
-    }
-
-	return 0;
-
+    complex<P>*cyv,complex<P>*cjvp,complex<P>*cyvp)
+{
+    //first, compute the bessel functions of fractional order
+    cbessjyva<P>(v + 0.5, z, vm, cjv, cyv, cjvp, cyvp);
+
+    //iterate through each and scale
+    for(int n = 0; n<=v; n++)
+    {
+
+        cjv[n] = cjv[n] * sqrt(rtsPI/(z * 2.0));
+        cyv[n] = cyv[n] * sqrt(rtsPI/(z * 2.0));
+
+        cjvp[n] = -1.0 / (z * 2.0) * cjv[n] + cjvp[n] * sqrt(rtsPI / (z * 2.0));
+        cyvp[n] = -1.0 / (z * 2.0) * cyv[n] + cyvp[n] * sqrt(rtsPI / (z * 2.0));
+    }
+
+	return 0;
+
 }
 
 }	//end namespace rts
diff --git a/math/complex.h b/math/complex.h
index 12f7975..e5414be 100644
--- a/math/complex.h
+++ b/math/complex.h
@@ -23,7 +23,14 @@ struct complex
     CUDA_CALLABLE complex()
     {
         r = 0;
-	i = 0;
+	   i = 0;
+    }
+
+    //constructor when given real and imaginary values
+    CUDA_CALLABLE complex(T r, T i = 0)
+    {
+        this->r = r;
+        this->i = i;
     }
 
 	//access methods
@@ -48,12 +55,7 @@ struct complex
 		return i_val;
 	}
 
-    //constructor when given real and imaginary values
-    CUDA_CALLABLE complex(T r, T i)
-    {
-        this->r = r;
-        this->i = i;
-    }
+    
 
     //return the current value multiplied by i
     CUDA_CALLABLE complex<T> imul()
diff --git a/math/constants.h b/math/constants.h
index 598c152..3dae3a5 100644
--- a/math/constants.h
+++ b/math/constants.h
@@ -1,10 +1,7 @@
 #ifndef RTS_CONSTANTS_H
 #define RTS_CONSTANTS_H
 
-namespace rts{
-
-	double PI = 3.14159;
-	double TAU = 2 * PI;
-}
+#define rtsPI 	3.14159
+#define rtsTAU	2 * rtsPI
 
 #endif
\ No newline at end of file
diff --git a/math/plane.h b/math/plane.h
new file mode 100644
index 0000000..f7753ac
--- /dev/null
+++ b/math/plane.h
@@ -0,0 +1,165 @@
+#ifndef RTS_PLANE_H
+#define RTS_PLANE_H
+
+#include <iostream>
+#include "../math/vector.h"
+#include "rts/cuda/callable.h"
+
+
+namespace rts{
+template <typename T, int D> class plane;
+}
+
+template <typename T, int D>
+CUDA_CALLABLE rts::plane<T, D> operator-(rts::plane<T, D> v);
+
+namespace rts{
+
+template <class T, int D = 3>
+class plane{
+
+	//a plane is defined by a point and a normal
+
+private:
+
+	vec<T, D> P;	//point on the plane
+	vec<T, D> N;	//plane normal
+
+	CUDA_CALLABLE void init(){
+		P = vec<T, D>(0, 0, 0);
+		N = vec<T, D>(0, 0, 1);
+	}
+
+
+public:
+
+	//default constructor
+	CUDA_CALLABLE plane(){
+		init();
+	}
+	
+	CUDA_CALLABLE plane(vec<T, D> n, vec<T, D> p = vec<T, D>(0, 0, 0)){
+		P = p;
+		N = n.norm();
+	}
+
+	CUDA_CALLABLE plane(T z_pos){
+		init();
+		P[2] = z_pos;
+	}
+
+	//create a plane from three points (a triangle)
+	CUDA_CALLABLE plane(vec<T, D> a, vec<T, D> b, vec<T, D> c){
+		P = c;
+		N = (c - a).cross(b - a);
+		if(N.len() == 0)	//handle the degenerate case when two vectors are the same, N = 0
+			N = 0;
+		else
+			N = N.norm();
+	}
+
+	CUDA_CALLABLE vec<T, D> norm(){
+		return N;
+	}
+
+	CUDA_CALLABLE vec<T, D> p(){
+		return P;
+	}
+
+	//flip the plane front-to-back
+	CUDA_CALLABLE plane<T, D> flip(){
+		plane<T, D> result = *this;
+		result.N = -result.N;
+		return result;
+	}
+
+	//determines how a vector v intersects the plane (1 = intersects front, 0 = within plane, -1 = intersects back)
+	CUDA_CALLABLE int face(vec<T, D> v){
+		
+		T dprod = v.dot(N);		//get the dot product between v and N
+
+		//conditional returns the appropriate value
+		if(dprod < 0)
+			return 1;
+		else if(dprod > 0)
+			return -1;
+		else
+			return 0;
+	}
+
+	//determine on which side of the plane a point lies (1 = front, 0 = on the plane, -1 = back)
+	CUDA_CALLABLE int side(vec<T, D> p){
+
+		vec<T, D> v = p - P;	//get the vector from P to the query point p
+
+		return face(v);
+	}
+
+	//compute the component of v that is perpendicular to the plane
+	CUDA_CALLABLE vec<T, D> perpendicular(vec<T, D> v){
+		return N * v.dot(N);
+	}
+
+	//compute the projection of v in the plane
+	CUDA_CALLABLE vec<T, D> parallel(vec<T, D> v){
+		return v - perpendicular(v);
+	}
+
+	//get both the parallel and perpendicular components of a vector v w.r.t. the plane
+	CUDA_CALLABLE void project(vec<T, D> v, vec<T, D> &v_par, vec<T, D> &v_perp){
+
+		v_perp = v.dot(N);
+		v_par = v - v_perp;
+	}
+
+	//compute the reflection of v off of the plane
+	CUDA_CALLABLE vec<T, D> reflect(vec<T, D> v){
+
+		//compute the reflection using N_prime as the plane normal
+		vec<T, D> par = parallel(v);
+		vec<T, D> r = (-v) + par * 2;
+
+		/*std::cout<<"----------------REFLECT-----------------------------"<<std::endl;
+		std::cout<<str()<<std::endl;
+		std::cout<<"v: "<<v<<std::endl;
+		std::cout<<"r: "<<r<<std::endl;
+		std::cout<<"Perpendicular: "<<perpendicular(v)<<std::endl;
+		std::cout<<"Parallel: "<<par<<std::endl;*/
+		return r;
+
+	}
+
+	//output a string
+	std::string str(){
+		std::stringstream ss;
+		ss<<"P: "<<P<<std::endl;
+		ss<<"N: "<<N;
+		return ss.str();
+	}
+
+	///////Friendship
+	friend CUDA_CALLABLE rts::plane<T, D> operator- <> (rts::plane<T, D> v);
+
+
+
+};
+
+}
+
+//arithmetic operators
+
+//negative operator flips the plane (front to back)
+template <typename T, int D>
+CUDA_CALLABLE rts::plane<T, D> operator-(rts::plane<T, D> p_rhs)
+{
+	rts::plane<T, D> p = p_rhs;
+
+	//negate the normal vector
+	p.N = -p.N;
+
+	return p;
+}
+
+
+
+#endif
\ No newline at end of file
diff --git a/math/realfield.cuh b/math/realfield.cuh
index 5a1d906..8e0c459 100644
--- a/math/realfield.cuh
+++ b/math/realfield.cuh
@@ -84,7 +84,7 @@ public:
 	{
 		R[0] = R[1] = 0;
 		init();
-		std::cout<<"realfield CONSTRUCTOR"<<std::endl;
+		//std::cout<<"realfield CONSTRUCTOR"<<std::endl;
 	}
 	realfield(unsigned int x, unsigned int y)
 	{
@@ -98,13 +98,13 @@ public:
 		}
 		//shape = rect<P>(vec<P>(-1, -1, 0), vec<P>(-1, 1, 0), vec<P>(1, 1, 0));	//default geometry
 		clear();		//zero the field
-		std::cout<<"realfield CONSTRUCTOR"<<std::endl;
+		//std::cout<<"realfield CONSTRUCTOR"<<std::endl;
     }
 
 	~realfield()
     {
 		destroy();
-		std::cout<<"realfield DESTRUCTOR"<<std::endl;
+		//std::cout<<"realfield DESTRUCTOR"<<std::endl;
     }
 
 	P* ptr(unsigned int n)
@@ -207,7 +207,7 @@ public:
 			//copy the slice
 			HANDLE_ERROR(cudaMemcpy(X[n], rhs.X[n], sizeof(P) * R[0] * R[1], cudaMemcpyDeviceToDevice));
 		}
-        std::cout<<"Assignment operator."<<std::endl;
+        //std::cout<<"Assignment operator."<<std::endl;
 
         return *this;
     }
@@ -234,7 +234,7 @@ public:
 			}
 		}
 
-		std::cout<<"realfield COPY CONSTRUCTOR"<<std::endl;
+		//std::cout<<"realfield COPY CONSTRUCTOR"<<std::endl;
 	}
 
 	/*void gaussian(P mean, P std, unsigned int n=0)	//creates a 3D gaussian using component n
diff --git a/math/vector.h b/math/vector.h
index 6af0789..fdf567b 100644
--- a/math/vector.h
+++ b/math/vector.h
@@ -26,9 +26,10 @@ struct vec
 	}
 
 	//efficiency constructor, makes construction easier for 1D-4D vectors
-	CUDA_CALLABLE vec(T x)
+	CUDA_CALLABLE vec(T rhs)
 	{
-		v[0] = x;
+		for(int i=0; i<N; i++)
+			v[i] = rhs;
 	}
 	CUDA_CALLABLE vec(T x, T y)
 	{
@@ -49,9 +50,10 @@ struct vec
 		v[3] = w;
 	}
 
-	CUDA_CALLABLE vec(const T(&data)[N])
-	{
-		memcpy(v, data, sizeof(T) * N);
+	//copy constructor
+	CUDA_CALLABLE vec( const vec<T, N>& other){
+		for(int i=0; i<N; i++)
+			v[i] = other.v[i];
 	}
 
 	CUDA_CALLABLE T len() const
@@ -135,12 +137,21 @@ struct vec
 	//arithmetic
 	CUDA_CALLABLE vec<T, N> operator+(vec<T, N> rhs) const
 	{
-        vec<T, N> result;
+		vec<T, N> result;
 
-        for(int i=0; i<N; i++)
-            result.v[i] = v[i] + rhs.v[i];
+		for(int i=0; i<N; i++)
+		    result.v[i] = v[i] + rhs.v[i];
 
-        return result;
+		return result;
+	}
+	CUDA_CALLABLE vec<T, N> operator+(T rhs) const
+	{
+		vec<T, N> result;
+
+		for(int i=0; i<N; i++)
+		    result.v[i] = v[i] + rhs;
+
+		return result;
 	}
 	CUDA_CALLABLE vec<T, N> operator-(vec<T, N> rhs) const
 	{
@@ -174,15 +185,21 @@ struct vec
 			v[i] = v[i] * rhs;
 		return *this;
 	}
+	CUDA_CALLABLE vec<T, N> & operator=(T rhs){
+		for(int i=0; i<N; i++)
+			v[i] = rhs;
+		return *this;
+	}
+	//unary minus
+	CUDA_CALLABLE vec<T, N> operator-() const{
+		vec<T, N> r;
 
-	//conversion from a point
-	/*CUDA_CALLABLE vector<T, N> & operator=(point<T, N> rhs)
-	{
-		for(int n=0; n<N; n++)
-			v[n] = rhs.p[n];
+		//negate the vector
+		for(int i=0; i<N; i++)
+		    r.v[i] = -v[i];
 
-		return *this;
-	}*/
+		return r;
+	}
 
 	CUDA_CALLABLE bool operator==(vec<T, N> rhs) const
 	{
@@ -226,18 +243,7 @@ std::ostream& operator<<(std::ostream& os, rts::vec<T, N> v)
     return os;
 }
 
-//arithmetic operators
-template <typename T, int N>
-CUDA_CALLABLE rts::vec<T, N> operator-(rts::vec<T, N> v)
-{
-    rts::vec<T, N> r;
-
-    //negate the vector
-    for(int i=0; i<N; i++)
-        r.v[i] = -v.v[i];
 
-    return r;
-}
 
 template <typename T, int N>
 CUDA_CALLABLE rts::vec<T, N> operator*(T lhs, rts::vec<T, N> rhs)
diff --git a/optics/beam.h b/optics/beam.h
index 9e86f4c..eb8ac15 100644
--- a/optics/beam.h
+++ b/optics/beam.h
@@ -16,10 +16,10 @@ public:
 
 private:
 	
-	P na[2];		//numerical aperature of the focusing optics	
+	P _na[2];		//numerical aperature of the focusing optics	
 	vec<P> f;		//focal point	
 	function<P, P> apod;	//apodization function
-	unsigned int apod_res;	//resolution of complex apodization filters
+	unsigned int apod_res;	//resolution of apodization filter functions
 
 	void apod_uniform()
 	{
@@ -69,28 +69,44 @@ public:
 
 	///constructor: build a default beam (NA=1.0)
 	beam(
-		vec<P> _k = rts::vec<P>(0, 0, TAU), 
+		vec<P> k = rts::vec<P>(0, 0, rtsTAU), 
 		vec<P> _E0 = rts::vec<P>(1, 0, 0), 
 		beam_type _apod = Uniform)
-		: planewave<P>(_k, _E0)
+		: planewave<P>(k, _E0)
 	{
-		na[0] = (P)0.0;
-		na[1] = (P)1.0;
+		_na[0] = (P)0.0;
+		_na[1] = (P)1.0;
 		f = vec<P>( (P)0, (P)0, (P)0 );
 		apod_res = 256;						//set the default resolution for apodization filters
 		set_apod(_apod);						//set the apodization function type
 	}
 
+	beam<P> refract(rts::vec<P> kn) const{
+
+		beam<P> new_beam;
+		new_beam._na[0] = _na[0];
+		new_beam._na[1] = _na[1];
+
+
+		rts::planewave<P> pw = planewave<P>::bend(kn);
+		//std::cout<<pw.str()<<std::endl;
+
+		new_beam.k = pw.kvec();
+		new_beam.E0 = pw.E();
+
+		return new_beam;
+	}
+
 	///Numerical Aperature functions
-	void NA(P _na)
+	void NA(P na)
 	{
-		na[0] = (P)0;
-		na[1] = _na;
+		_na[0] = (P)0;
+		_na[1] = na;
 	}
-	void NA(P _na0, P _na1)
+	void NA(P na0, P na1)
 	{
-		na[0] = _na0;
-		na[1] = _na1;
+		_na[0] = na0;
+		_na[1] = na1;
 	}
 
 	/*string str() : 
@@ -116,7 +132,12 @@ public:
 		///compute the rotation operator to transform (0, 0, 1) to k
 		P cos_angle = k_hat.dot(rts::vec<P>(0, 0, 1));
 		rts::matrix<P, 3> rotation;
-		if(cos_angle != 1.0)
+
+		//if the cosine of the angle is -1, the rotation is just a flip across the z axis
+		if(cos_angle == -1){
+			rotation(2, 2) = -1;
+		}
+		else if(cos_angle != 1.0)
 		{
 			rts::vec<P> r_axis = rts::vec<P>(0, 0, 1).cross(k_hat).norm();	//compute the axis of rotation
 			P angle = acos(cos_angle);							//compute the angle of rotation
@@ -127,8 +148,8 @@ public:
 
 		//find the phi values associated with the cassegrain ring
 		P PHI[2];
-		PHI[0] = (P)asin(na[0]);
-		PHI[1] = (P)asin(na[1]);
+		PHI[0] = (P)asin(_na[0]);
+		PHI[1] = (P)asin(_na[1]);
 
 		//calculate the z-axis cylinder coordinates associated with these angles
 		P Z[2];
@@ -152,12 +173,28 @@ public:
 			vec<P> k_prime = rotation * cart;				//create a sample vector
 
 			//store a wave refracted along the given direction
-			samples.push_back(beam::refract(k_prime) * apod(phi/PHI[1]));
+			//std::cout<<"k prime: "<<rotation<<std::endl;
+			samples.push_back(planewave<P>::refract(k_prime) * apod(phi/PHI[1]));
 		}
 
 		return samples;
 	}
 
+	std::string str()
+	{
+		std::stringstream ss;
+		ss<<"Beam:"<<std::endl;
+		ss<<"	Central Plane Wave: "<<beam::E0<<" e^i ( "<<beam::k<<" . r )"<<std::endl;
+		if(_na[0] == 0)
+			ss<<"	NA: "<<_na[1];
+		else
+			ss<<"	NA: "<<_na[0]<<" -- "<<_na[1];
+
+		return ss.str();
+	}
+
+
+
 };
 
 }
diff --git a/optics/efield.cuh b/optics/efield.cuh
index 65d6e1c..d710035 100644
--- a/optics/efield.cuh
+++ b/optics/efield.cuh
@@ -8,7 +8,14 @@
 #include "../cuda/devices.h"
 #include "../optics/beam.h"
 
+namespace rts{
+template<typename T> class halfspace;
 
+template<typename T> class efield;
+}
+
+template<typename T>
+void operator<<(rts::efield<T> &ef, rts::halfspace<T> hs);
 
 namespace rts{
 
@@ -63,6 +70,27 @@ __global__ void gpu_efield_magnitude(complex<T>* X, complex<T>* Y, complex<T>* Z
 }
 
 template<typename T>
+__global__ void gpu_efield_real_magnitude(complex<T>* X, complex<T>* Y, complex<T>* Z, unsigned int r0, unsigned int r1, T* M)
+{
+    int iu = blockIdx.x * blockDim.x + threadIdx.x;
+    int iv = blockIdx.y * blockDim.y + threadIdx.y;
+
+    //make sure that the thread indices are in-bounds
+    if(iu >= r0 || iv >= r1) return;
+
+    //compute the index into the field
+    int i = iv*r0 + iu;
+
+    if(Y == NULL)                      //if this is a scalar simulation
+       M[i] = X[i].abs();              //compute the magnitude of the X component
+    else
+    {
+        M[i] = rts::vec<T>(X[i].real(), Y[i].real(), Z[i].real()).len();
+        //M[i] = Z[i].abs();
+    }
+}
+
+template<typename T>
 __global__ void gpu_efield_polarization(complex<T>* X, complex<T>* Y, complex<T>* Z, 
                                         unsigned int r0, unsigned int r1,
                                         T* Px, T* Py, T* Pz)
@@ -79,6 +107,7 @@ __global__ void gpu_efield_polarization(complex<T>* X, complex<T>* Y, complex<T>
     //compute the field polarization
     Px[i] = X[i].abs();
     Py[i] = Y[i].abs();
+
     Pz[i] = Z[i].abs();
 }
 
@@ -146,7 +175,6 @@ protected:
 
         if(vector)
         {
-            std::cout<<"vector:";
             cudaMalloc(&Y, sizeof(rts::complex<P>) * R[0] * R[1]);
 			cudaMemset(Y, 0, sizeof(rts::complex<P>) * R[0] * R[1]);
 
@@ -280,6 +308,7 @@ public:
 		return *this;		//return this object
 	}
 
+    //assignment operator: build an electric field from a beam
     efield<P> & operator= (const rts::beam<P> & rhs)
     {
         //get a vector of monte-carlo samples
@@ -311,6 +340,25 @@ public:
 		return M;
     }
 
+    //return a sacalar field representing the field magnitude at an infinitely small point in time
+    realfield<P, 1, true> real_mag()
+    {
+        int maxThreads = rts::maxThreadsPerBlock(); //compute the optimal block size
+        int SQRT_BLOCK = (int)std::sqrt((float)maxThreads);
+
+        //create a scalar field to store the result
+        realfield<P, 1, true> M(R[0], R[1]);
+
+        //create one thread for each detector pixel
+        dim3 dimBlock(SQRT_BLOCK, SQRT_BLOCK);
+        dim3 dimGrid((R[0] + SQRT_BLOCK -1)/SQRT_BLOCK, (R[1] + SQRT_BLOCK - 1)/SQRT_BLOCK);
+
+        //compute the magnitude and store it in a scalar field
+        gpu_efield_real_magnitude<float> <<<dimGrid, dimBlock>>> (X, Y, Z, R[0], R[1], M.ptr(0));
+
+        return M;
+    }
+
     //return a vector field representing field polarization
     realfield<P, 3, true> polarization()
     {
@@ -335,6 +383,8 @@ public:
         return Pol;                         //return the vector field
     }
 
+    ////////FRIENDSHIP
+    friend void operator<< <>(rts::efield<P> &ef, rts::halfspace<P> hs);
 
 
 };
diff --git a/optics/esphere.cuh b/optics/esphere.cuh
index 066102d..5ac31c1 100644
--- a/optics/esphere.cuh
+++ b/optics/esphere.cuh
@@ -134,7 +134,7 @@ public:
 		calcAB(rhs.kmag(), Nl, A, B);	//compute the scattering coefficients
 
 		//determine important parameters for the scattering domain
-		unsigned int sR = ceil(sqrt( (P)(pow(esphere::R[0],2) + pow(esphere::R[1],2))) );
+		unsigned int sR = ceil(sqrt( (P)(pow((P)esphere::R[0],2) + pow((P)esphere::R[1],2))) );
 		unsigned int thetaR = 256;
 
 		/////////////////////continue scattering code here/////////////////////////
diff --git a/optics/halfspace.cuh b/optics/halfspace.cuh
new file mode 100644
index 0000000..1b8e02e
--- /dev/null
+++ b/optics/halfspace.cuh
@@ -0,0 +1,91 @@
+#ifndef	RTS_HALFSPACE_CUH
+#define	RTS_HALFSPACE_CUH
+
+#include "../optics/halfspace.h"
+#include "../optics/efield.cuh"
+
+namespace rts{
+
+template<typename T>
+__global__ void gpu_halfspace_pw2ef(complex<T>* X, complex<T>* Y, complex<T>* Z, unsigned int r0, unsigned int r1, 
+									 plane<T> P, planewave<T> w, rect<T> q, bool at_surface = false)
+{
+    int iu = blockIdx.x * blockDim.x + threadIdx.x;
+    int iv = blockIdx.y * blockDim.y + threadIdx.y;
+
+    //make sure that the thread indices are in-bounds
+    if(iu >= r0 || iv >= r1) return;
+
+    //compute the index into the field
+    int i = iv*r0 + iu;
+
+	//get the current position
+	vec<T> p = q( (T)iu/(T)r0, (T)iv/(T)r1 );
+
+	if(at_surface){
+		if(P.side(p) > 0)
+			return;
+	}
+	else{
+		if(P.side(p) >= 0)
+			return;
+	}
+
+	//if the current position is on the wrong side of the plane
+
+	//vec<T> r(p[0], p[1], p[2]);
+
+	complex<T> x( 0.0f, w.kvec().dot(p) );
+	complex<T> phase( 0.0f, w.phase());
+
+    if(Y == NULL)                       //if this is a scalar simulation
+        X[i] += w.E().len() * exp(x);    //use the vector magnitude as the plane wave amplitude
+    else
+    {
+    	//X[i] = Y[i] = Z[i] = 1;
+        X[i] += w.E()[0] * exp(x) * exp(phase);
+        Y[i] += w.E()[1] * exp(x) * exp(phase);
+        Z[i] += w.E()[2] * exp(x) * exp(phase);
+    }
+}
+}
+
+template<class T>
+void operator<<(rts::efield<T> &ef, rts::halfspace<T> hs){
+
+	//std::cout<<"---------RENDER HALFSPACE--------------"<<std::endl;
+	//output a parameter of the efield
+	//std::cout<<ef.pos<<std::endl;
+	//std::cout<<hs.S.str()<<std::endl;
+
+	unsigned int SQRT_BLOCK = 16;
+	//create one thread for each detector pixel
+	dim3 dimBlock(SQRT_BLOCK, SQRT_BLOCK);
+	dim3 dimGrid((ef.R[0] + SQRT_BLOCK -1)/SQRT_BLOCK, (ef.R[1] + SQRT_BLOCK - 1)/SQRT_BLOCK);
+
+
+
+	//render each plane wave
+
+	//plane waves at the surface front
+	for(unsigned int w = 0; w < hs.w0.size(); w++){
+		//std::cout<<"w0 "<<w<<": "<<hs.w0[w].str()<<std::endl;
+		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);
+	}
+
+	//plane waves at the surface back
+	for(unsigned int w = 0; w < hs.w1.size(); w++){
+		//std::cout<<"w1 "<<w<<": "<<hs.w1[w].str()<<std::endl;
+		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);
+	}
+	//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);
+	//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);
+	//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);
+
+
+	//return ef;
+}
+
+
+
+#endif
\ No newline at end of file
diff --git a/optics/halfspace.h b/optics/halfspace.h
new file mode 100644
index 0000000..bb3c1e3
--- /dev/null
+++ b/optics/halfspace.h
@@ -0,0 +1,149 @@
+#ifndef	RTS_HALFSPACE_H
+#define	RTS_HALFSPACE_H
+
+namespace rts{
+template<typename T> class halfspace;
+
+template<typename T> class efield;
+}
+
+template<typename T>
+void operator<<(rts::efield<T> &ef, rts::halfspace<T> hs);
+
+namespace rts{
+
+template<class T>
+class halfspace
+{
+private:
+	rts::plane<T> S;		//surface plane splitting the half space
+	rts::complex<T> n0;		//refractive index at the front of the plane
+	rts::complex<T> n1;		//refractive index at the back of the plane
+
+	//lists of waves in front (pw0) and behind (pw1) the plane
+	std::vector< rts::planewave<T> > w0;
+	std::vector< rts::planewave<T> > w1;
+
+	//rts::planewave<T> pi;	//incident plane wave
+	//rts::planewave<T> pr;	//plane wave reflected from the surface
+	//rts::planewave<T> pt;	//plane wave transmitted through the surface
+
+	void init(){
+		n0 = 1.0;
+		n1 = 1.5;
+	}
+
+	//compute which side of the interface is hit by the incoming plane wave (0 = front, 1 = back)
+	bool facing(planewave<T> p){
+		if(p.kvec().dot(S.norm()) > 0)
+			return 1;
+		else
+			return 0;
+	}
+
+	T calc_theta_i(vec<T> v){
+
+		vec<T> v_hat = v.norm();
+
+		//compute the cosine of the angle between k_hat and the plane normal
+		T cos_theta_i = v_hat.dot(S.norm());
+
+		return acos(abs(cos_theta_i));
+	}
+
+	T calc_theta_t(T ni_nt, T theta_i){
+
+		T sin_theta_t = ni_nt * sin(theta_i);
+		return asin(sin_theta_t);
+	}
+
+
+public:
+
+	//constructors
+	halfspace(){
+		init();
+	}
+
+	halfspace(T na, T nb){
+		init();
+		n0 = na;
+		n1 = nb;
+	}
+
+	halfspace(T na, T nb, plane<T> p){
+		n0 = na;
+		n1 = nb;
+		S = p;
+	}
+
+	//compute the transmitted and reflective waves given the incident (vacuum) plane wave p
+	void incident(rts::planewave<T> p){
+
+		planewave<T> r, t;
+		p.scatter(S, n1.real()/n0.real(), r, t);
+
+		//std::cout<<"i+r: "<<p.r()[0] + r.r()[0]<<std::endl;
+		//std::cout<<"t:   "<<t.r()[0]<<std::endl;
+
+		if(facing(p)){
+			w1.push_back(p);
+
+			if(r.E().len() != 0)
+				w1.push_back(r);
+			if(t.E().len() != 0)
+				w0.push_back(t);
+		}
+		else{
+			w0.push_back(p);
+
+			if(r.E().len() != 0)
+				w0.push_back(r);
+			if(t.E().len() != 0)
+				w1.push_back(t);
+		}
+	}
+
+	void incident(rts::beam<T> b, unsigned int N = 10000){
+
+		//generate a plane wave decomposition for the beam
+		std::vector< planewave<T> > pw_list = b.mc(N);
+
+		//calculate the reflected and refracted waves for each incident wave
+		for(unsigned int w = 0; w < pw_list.size(); w++){
+			incident(pw_list[w]);
+		}
+	}
+
+	std::string str(){
+		std::stringstream ss;
+		ss<<"Half Space-------------"<<std::endl;
+		ss<<"n0: "<<n0<<std::endl;
+		ss<<"n1: "<<n1<<std::endl;
+		ss<<S.str();
+		
+
+		if(w0.size() != 0 || w1.size() != 0){
+			ss<<std::endl<<"Plane Waves:"<<std::endl;
+			for(unsigned int w = 0; w < w0.size(); w++){
+				ss<<"w0 = "<<w<<": "<<w0[w]<<std::endl;
+			}
+			for(unsigned int w = 0; w < w1.size(); w++){
+				ss<<"w1 = "<<w<<": "<<w1[w]<<std::endl;
+			}
+		}
+
+		return ss.str();
+	}
+	////////FRIENDSHIP
+    friend void operator<< <> (rts::efield<T> &ef, rts::halfspace<T> hs);
+
+};
+
+
+
+
+}
+
+
+#endif
\ No newline at end of file
diff --git a/optics/planewave.h b/optics/planewave.h
index e852164..390a2e9 100644
--- a/optics/planewave.h
+++ b/optics/planewave.h
@@ -7,6 +7,8 @@
 #include "../math/vector.h"
 #include "../math/quaternion.h"
 #include "../math/constants.h"
+#include "../math/plane.h"
+#include "rts/cuda/callable.h"
 
 /*Basic conversions used here (assuming a vacuum)
 	lambda =
@@ -14,160 +16,299 @@
 
 namespace rts{
 
-template<typename P>
+template<typename T>
 class planewave{
 
+protected:
+
+	vec<T> k;	//k = tau / lambda
+	vec<T> E0;		//amplitude
+	T phi;
+
+	planewave<T> bend(rts::vec<T> kn) const{
+
+		vec<T> kn_hat = kn.norm();				//normalize the new k
+		vec<T> k_hat = k.norm();				//normalize the current k
+
+		//std::cout<<"PLANE WAVE BENDING------------------"<<std::endl;
+		//std::cout<<"kn_hat: "<<kn_hat<<"     k_hat: "<<k_hat<<std::endl;
+
+		planewave<T> new_p;						//create a new plane wave
+
+		//if kn is equal to k or -k, handle the degenerate case
+		T k_dot_kn = k_hat.dot(kn_hat);
+
+		//if k . n < 0, then the bend is a reflection
+			//flip k_hat
+		if(k_dot_kn < 0) k_hat = -k_hat;
+
+		//std::cout<<"k dot kn: "<<k_dot_kn<<std::endl;
+
+		//std::cout<<"k_dot_kn: "<<k_dot_kn<<std::endl;
+		if(k_dot_kn == -1){
+			new_p.k = -k;
+			new_p.E0 = E0;
+			return new_p;
+		}
+		else if(k_dot_kn == 1){
+			new_p.k = k;
+			new_p.E0 = E0;
+			return new_p;
+		}
+
+		vec<T> r = k_hat.cross(kn_hat);			//compute the rotation vector
+
+		//std::cout<<"r: "<<r<<std::endl;
+
+		T theta = asin(r.len());				//compute the angle of the rotation about r
+
+		
+
+		//deal with a zero vector (both k and kn point in the same direction)
+		//if(theta == (T)0)
+		//{
+		//	new_p = *this;
+		//	return new_p;
+		//}
+
+		//create a quaternion to capture the rotation
+		quaternion<T> q;
+		q.CreateRotation(theta, r.norm());
+
+		//apply the rotation to E0
+		vec<T> E0n = q.toMatrix3() * E0;
+
+		new_p.k = kn_hat * kmag();
+		new_p.E0 = E0n;
+
+		return new_p;
+	}
+
 public:
-	vec<P> k;	//k = tau / lambda
-	vec<P> E0;		//amplitude
 
 
 	///constructor: create a plane wave propagating along z, polarized along x
-	/*planewave(P lambda = (P)1)
+	/*planewave(T lambda = (T)1)
 	{
-		k = rts::vec<P>(0, 0, 1) * (TAU/lambda);
-		E0 = rts::vec<P>(1, 0, 0);
+		k = rts::vec<T>(0, 0, 1) * (TAU/lambda);
+		E0 = rts::vec<T>(1, 0, 0);
 	}*/
-	///constructor: create a plane wave propagating along _k, polarized along _E0, at frequency _omega
-	planewave(vec<P> _k = rts::vec<P>(0, 0, TAU), vec<P> _E0 = rts::vec<P>(1, 0, 0))
+	///constructor: create a plane wave propagating along k, polarized along _E0, at frequency _omega
+	CUDA_CALLABLE planewave(vec<T> kvec = rts::vec<T>(0, 0, rtsTAU), 
+							vec<T> E = rts::vec<T>(1, 0, 0), 
+							T phase = 0)
 	{
-		k = _k;
-		vec<P> k_hat = _k.norm();
-
-		vec<P> s = (k.cross(_E0)).norm();		//compute an orthogonal side vector
-		vec<P> E0_hat = (s.cross(k)).norm();	//compute a normalized E0 direction vector
-		E0 = E0_hat * E0_hat.dot(_E0);					//compute the projection of _E0 onto E0_hat
+		phi = phase;
+		k = kvec;
+		vec<T> k_hat = k.norm();
+
+		if(E.len() == 0)			//if the plane wave has an amplitude of 0
+			E0 = vec<T>(0);			//just return it
+		else{
+			vec<T> s = (k.cross(E)).norm();		//compute an orthogonal side vector
+			vec<T> E_hat = (s.cross(k)).norm();	//compute a normalized E0 direction vector
+			E0 = E_hat * E_hat.dot(E);					//compute the projection of _E0 onto E0_hat
+		}
 	}
 
 	///multiplication operator: scale E0
-    planewave<P> & operator* (const P & rhs)
+    CUDA_CALLABLE planewave<T> & operator* (const T & rhs)
 	{
 		
 		E0 = E0 * rhs;
 		return *this;
 	}
 
-	P lambda() const
+	CUDA_CALLABLE T lambda() const
 	{
-		return TAU / k.len();
+		return rtsTAU / k.len();
 	}
 
-	P kmag() const
+	CUDA_CALLABLE T kmag() const
 	{
 		return k.len();
 	}
 
+	CUDA_CALLABLE vec<T> E(){
+		return E0;
+	}
+
+	CUDA_CALLABLE vec<T> kvec(){
+		return k;
+	}
+
+	CUDA_CALLABLE T phase(){
+		return phi;
+	}
+
+	CUDA_CALLABLE void phase(T p){
+		phi = p;
+	}
+
+	CUDA_CALLABLE vec< complex<T> > pos(vec<T> p = vec<T>(0, 0, 0)){
+		vec< complex<T> > result;
+
+		T kdp = k.dot(p);
+		complex<T> x = complex<T>(0, kdp);
+		complex<T> expx = exp(x);
+
+		result[0] = E0[0] * expx;
+		result[1] = E0[1] * expx;
+		result[2] = E0[2] * expx;
+
+		return result;
+	}
+
+	//scales k based on a transition from material ni to material nt
+	CUDA_CALLABLE planewave<T> n(T ni, T nt){
+		return planewave<T>(k * (nt / ni), E0);
+	}
 
-	planewave<P> refract(rts::vec<P> kn) const
+	CUDA_CALLABLE planewave<T> refract(rts::vec<T> kn) const
 	{
-		vec<P> kn_hat = kn.norm();				//normalize the new k
-		vec<P> k_hat = k.norm();				//normalize the current k
+		return bend(kn);
+	}
 
-		vec<P> r = k_hat.cross(kn_hat);			//compute the rotation vector
+	void scatter(rts::plane<T> P, T nr, planewave<T> &r, planewave<T> &t){
 
-		P theta = asin(r.len());				//compute the angle of the rotation about r
+		int facing = P.face(k);		//determine which direction the plane wave is coming in
 
-		planewave<P> new_p;						//create a new plane wave
+		//if(facing == 0)				//if the wave is tangent to the plane, return an identical wave
+		//	return *this;
+		//else 
+		if(facing == -1){		//if the wave hits the back of the plane, invert the plane and nr
+			P = P.flip();			//flip the plane
+			nr = 1/nr;				//invert the refractive index (now nr = n0/n1)
+		}
 
-		//deal with a zero vector (both k and kn point in the same direction)
-		if(theta == (P)0)
-		{
-			new_p = *this;
-			return new_p;
+		//use Snell's Law to calculate the transmitted angle
+		T cos_theta_i = k.norm().dot(-P.norm());				//compute the cosine of theta_i
+		T theta_i = acos(cos_theta_i);							//compute theta_i
+		T sin_theta_t = (1/nr) * sin(theta_i);						//compute the sine of theta_t using Snell's law
+		T theta_t = asin(sin_theta_t);							//compute the cosine of theta_t
+
+		bool tir = false;						//flag for total internal reflection
+		if(theta_t != theta_t){
+			tir = true;
+			theta_t = rtsPI / (T)2;
 		}
 
-		//create a quaternion to capture the rotation
-		quaternion<P> q;
-		q.CreateRotation(theta, r.norm());
+		//handle the degenerate case where theta_i is 0 (the plane wave hits head-on)
+		if(theta_i == 0){
+			T rp = (1 - nr) / (1 + nr);		//compute the Fresnel coefficients
+			T tp = 2 / (1 + nr);
+			vec<T> kr = -k;
+			vec<T> kt = k * nr;			//set the k vectors for theta_i = 0
+			vec<T> Er = E0 * rp;		//compute the E vectors
+			vec<T> Et = E0 * tp;
+			T phase_t = P.p().dot(k - kt);	//compute the phase offset
+			T phase_r = P.p().dot(k - kr);
+			//std::cout<<"Degeneracy: Head-On"<<std::endl;
+			//std::cout<<"rs: "<<rp<<"  rp: "<<rp<<"  ts: "<<tp<<"  tp: "<<tp<<std::endl;
+			//std::cout<<"phase r: "<<phase_r<<"  phase t: "<<phase_t<<std::endl;
+
+			//create the plane waves
+			r = planewave<T>(kr, Er, phase_r);
+			t = planewave<T>(kt, Et, phase_t);
+
+			//std::cout<<"i + r: "<<pos()[0] + r.pos()[0]<<pos()[1] + r.pos()[1]<<pos()[2] + r.pos()[2]<<std::endl;
+			//std::cout<<"t:     "<<t.pos()[0]<<t.pos()[1]<<t.pos()[2]<<std::endl;
+			//std::cout<<"--------------------------------"<<std::endl;
+			return;
+		}
 
-		//apply the rotation to E0
-		vec<P> E0n = q.toMatrix3() * E0;
 
-		new_p.k = kn_hat * k.len();
-		new_p.E0 = E0n;
+		//compute the Fresnel coefficients
+		T rp, rs, tp, ts;
+		rp = tan(theta_t - theta_i) / tan(theta_t + theta_i);
+		rs = sin(theta_t - theta_i) / sin(theta_t + theta_i);
+		
+		if(tir){
+			tp = ts = 0;
+		}
+		else{
+			tp = ( 2 * sin(theta_t) * cos(theta_i) ) / ( sin(theta_t + theta_i) * cos(theta_t - theta_i) );
+			ts = ( 2 * sin(theta_t) * cos(theta_i) ) / sin(theta_t + theta_i);
+		}
 
-		return new_p;
+		//compute the coordinate space for the plane of incidence
+		vec<T> z_hat = -P.norm();
+		vec<T> y_hat = P.parallel(k).norm();
+		vec<T> x_hat = y_hat.cross(z_hat).norm();
+
+		//compute the k vectors for r and t
+		vec<T> kr, kt;
+		kr = ( y_hat * sin(theta_i) - z_hat * cos(theta_i) ) * kmag();
+		kt = ( y_hat * sin(theta_t) + z_hat * cos(theta_t) ) * kmag() * nr;
+
+		//compute the magnitude of the p- and s-polarized components of the incident E vector
+		T Ei_s = E0.dot(x_hat);
+		int sgn = (0 < E0.dot(y_hat)) - (E0.dot(y_hat) < 0);
+		T Ei_p = sgn * ( E0 - x_hat * Ei_s ).len();
+		//T Ei_p = ( E0 - x_hat * Ei_s ).len();
+		//compute the magnitude of the p- and s-polarized components of the reflected E vector
+		T Er_s = Ei_s * rs;
+		T Er_p = Ei_p * rp;
+		//compute the magnitude of the p- and s-polarized components of the transmitted E vector
+		T Et_s = Ei_s * ts;
+		T Et_p = Ei_p * tp;
+
+		//std::cout<<"E0: "<<E0<<std::endl;
+		//std::cout<<"E0 dot y_hat: "<<E0.dot(y_hat)<<std::endl;
+		//std::cout<<"theta i: "<<theta_i<<"  theta t: "<<theta_t<<std::endl;
+		//std::cout<<"x_hat: "<<x_hat<<"  y_hat: "<<y_hat<<"  z_hat: "<<z_hat<<std::endl;
+		//std::cout<<"Ei_s: "<<Ei_s<<"  Ei_p: "<<Ei_p<<"  Er_s: "<<Er_s<<"  Er_p: "<<Er_p<<"  Et_s: "<<Et_s<<"  Et_p: "<<Et_p<<std::endl;
+		//std::cout<<"rs: "<<rs<<"  rp: "<<rp<<"  ts: "<<ts<<"  tp: "<<tp<<std::endl;
+		
+
+		//compute the reflected E vector
+		vec<T> Er = (y_hat * cos(theta_i) + z_hat * sin(theta_i)) * Er_p + x_hat * Er_s;
+		//compute the transmitted E vector
+		vec<T> Et = (y_hat * cos(theta_t) - z_hat * sin(theta_t)) * Et_p + x_hat * Et_s;
+
+		T phase_t = P.p().dot(k - kt);
+		T phase_r = P.p().dot(k - kr);
+
+		//std::cout<<"phase r: "<<phase_r<<"  phase t: "<<phase_t<<std::endl;
+
+		//std::cout<<"phase: "<<phase<<std::endl;
+
+		//create the plane waves
+		r.k = kr;
+		r.E0 = Er;
+		r.phi = phase_r;
+
+		//t = bend(kt);
+		//t.k = t.k * nr;
+
+		t.k = kt;
+		t.E0 = Et;
+		t.phi = phase_t;
+		//std::cout<<"i: "<<str()<<std::endl;
+		//std::cout<<"r: "<<r.str()<<std::endl;
+		//std::cout<<"t: "<<t.str()<<std::endl;
+
+		//std::cout<<"i + r: "<<pos()[0] + r.pos()[0]<<pos()[1] + r.pos()[1]<<pos()[2] + r.pos()[2]<<std::endl;
+		//std::cout<<"t:     "<<t.pos()[0]<<t.pos()[1]<<t.pos()[2]<<std::endl;
+		//std::cout<<"--------------------------------"<<std::endl;
 
-		/*vec<P> kn_hat = kn.norm();	//normalize new_k
-		vec<P> k_hat = k.norm();
-
-		//compute the side vector (around which we will be rotating)
-		vec<P> E0_hat = E0.norm();
-		rts::vec<P> s = k_hat.cross(E0_hat);
-
-		//compute the projection of k' onto the k-E plane
-		rts::vec<P> s_prime = s * (kn_hat.dot(s));
-
-		//compute the angle between
-		vec<P> kp_hat = (kn_hat - s_prime).norm();
-		P theta = acos(k_hat.dot( kp_hat ));
-		if(kn_hat.dot(E0_hat) < 0)
-			theta = -theta;
-
-		//rotate E0 around s by theta
-		quaternion<P> q;
-		q.CreateRotation(theta, s);
-		rts::vec<P> E0_prime = q.toMatrix3() * E0;
-
-		//create the refracted plane wave
-		planewave<P> new_p;
-		new_p.E0 = E0_prime;
-		new_p.k = kn_hat * k.len();
-
-		return new_p;*/
-
-		/*vec<P> kn_hat = kn.norm();		//normalize kn
-		vec<P> k_hat = k.norm();
-		vec<P> E0_hat = E0.norm();
-
-		vec<P> B = k_hat.cross(E0_hat);
-		planewave<P> newp;
-		newp.k = kn_hat * k.len();
-		newp.E0 = B.cross(kn_hat).norm();
-		std::cout<<newp.k.norm().dot(newp.E0.norm())<<std::endl;
-		return newp;*/
-/*
-		//compute the side vector (around which we will be rotating)
-		rts::vec<P> s_hat = k_hat.cross(E0_hat);
-		//std::cout<<s.len()<<std::endl;
-
-		//project kn_hat into the k-E0 plane
-		rts::vec<P> sp = s_hat * (kn_hat.dot(s_hat));	//project k_new onto s
-		rts::vec<P> kp = (kn_hat - sp);	//correct k_new so it lies on the E0-k plane
-		rts::vec<P> kp_hat = kp.norm();
-
-		//compute the angle and direction between k_prime and k
-		P theta = acos(k_hat.dot(kp_hat));
-		if(kp_hat.dot(E0_hat) < 0)
-			theta = -theta;
-
-		//rotate E0 around s by theta
-		quaternion<P> q;
-		q.CreateRotation(theta, s_hat);
-		rts::vec<P> E0n = q.toMatrix3() * E0;
-		rts::vec<P> E0n_hat = E0n.norm();
-
-		//std::cout<<s_hat.dot(kp_hat)<<"  "<<s_hat.dot(E0n_hat)<<"  "<<s_hat.dot(E0_hat)<<"  "<<s_hat.dot(k_hat)<<"  "<<
-		//	E0_hat.dot(k_hat)<<"  "<<k_hat.dot(kp_hat)<<"  "<<E0_hat.dot(E0n_hat)<<"  "<<E0n_hat.dot(kp_hat)<<std::endl;
-
-		//create the refracted plane wave
-		//std::cout<<"cos: "<<cos(theta)<<"   k*kp: "<<k_hat.dot(kp_hat)<<"  E0*E0p: "<<E0_hat.dot(E0n_hat)<<"  E0p*kp: "<<E0n_hat.dot(kp_hat)<<std::endl;
-
-		//std::cout<<"kp*s: "<<kp.dot(s_hat)<<"   E0n*s: "<<E0n.dot(s_hat)<<"  |E0n|: "<<E0n.len()<<"  E0n*kp: "<<E0n.dot(kp_hat)<<"  E0n*kn: "<<E0n.dot(kn_hat)<<std::endl;
-
-		planewave<P> new_p(kn_hat * k.len(), E0n);				//create the plane wave
-		std::cout<<"|E0n|: "<<new_p.E0.len()<<"  E0n*kn: "<<(new_p.E0.norm()).dot(new_p.k.norm())<<std::endl;
-
-		return new_p;*/
 	}
 
 	std::string str()
 	{
 		std::stringstream ss;
-		ss<<E0<<" e^i ( "<<k<<" . r )";
+		ss<<"Plane Wave:"<<std::endl;
+		ss<<"	"<<E0<<" e^i ( "<<k<<" . r )";
 		return ss.str();
 	}
 };
 }
 
+template <typename T>
+std::ostream& operator<<(std::ostream& os, rts::planewave<T> p)
+{
+    os<<p.str();
+    return os;
+}
+
 #endif
\ No newline at end of file
--
libgit2 0.21.4