edits to update scattering code

David Mayerich
1 parent 4fcb3a86
Showing 6 changed files with 420 additions and 101 deletions Show diff stats
cmake/FindCNPY.cmake
python/optics.py
stim/math/matrix_sq.h
stim/math/vec3.h
stim/math/vector.h
stim/optics/planewave.h
+# finds the STIM library (downloads it if it isn't present)
+# set STIMLIB_PATH to the directory containing the stim subdirectory (the stim repository)
+
+include(FindPackageHandleStandardArgs)
+
+
+
+IF(NOT UNIX)
+    set(CNPY_ROOT $ENV{PROGRAMFILES}/CNPY)
+    find_path(CNPY_INCLUDE_DIR PATHS ${CNPY_ROOT} NAMES cnpy.h)
+    find_library(CNPY_LIBRARY NAMES cnpy PATHS ${CNPY_ROOT})
+ENDIF(NOT UNIX)
+
+
+find_package_handle_standard_args(CNPY DEFAULT_MSG CNPY_LIBRARY CNPY_INCLUDE_DIR)
+
+set(CNPY_LIBRARIES ${CNPY_LIBRARY})
+set(CNPY_INCLUDE_DIRS ${CNPY_INCLUDE_DIR})
 \ No newline at end of file
@@ -40,12 +40,12 @@ class planewave:
         self.k = np.array(k)              
         self.E = E
-        if ( np.linalg.norm(k) > 1e-15 and np.linalg.norm(E) >1e-15):
-            s = np.cross(k, E)              #compute an orthogonal side vector
-            s = s / np.linalg.norm(s)       #normalize it
-            Edir = np.cross(s, k)              #compute new E vector which is orthogonal
-            self.k = np.array(k)
-            self.E = Edir / np.linalg.norm(Edir) * np.linalg.norm(E)
+        #if ( np.linalg.norm(k) > 1e-15 and np.linalg.norm(E) >1e-15):
+        #    s = np.cross(k, E)              #compute an orthogonal side vector
+        #    s = s / np.linalg.norm(s)       #normalize it
+        #    Edir = np.cross(s, k)              #compute new E vector which is orthogonal
+        #    self.k = np.array(k)
+        #    self.E = Edir / np.linalg.norm(Edir) * np.linalg.norm(E)
     def __str__(self):
         return str(self.k) + "\n" + str(self.E)     #for verify field vectors use print command
@@ -100,7 +100,7 @@ struct matrix_sq
 	template<typename Y>
 	CUDA_CALLABLE vec3<Y> operator*(vec3<Y> rhs){
-		vec3<Y> result = 0;
+		vec3<Y> result;
 		for(int r=0; r<3; r++)
 			for(int c=0; c<3; c++)
 				result[r] += (*this)(r, c) * rhs[c];
@@ -3,6 +3,7 @@
 #include <stim/cuda/cudatools/callable.h>
+#include <complex>
 #include <cmath>
 #include <sstream>
@@ -43,21 +44,25 @@ public:
 	CUDA_CALLABLE T& operator[](size_t idx){
 		return ptr[idx];
 	}
+	//read only accessor method
+	CUDA_CALLABLE T get(size_t idx) const {
+		return ptr[idx];
+	}
 	CUDA_CALLABLE T* data(){
 		return ptr;
 	}
 /// Casting operator. Creates a new vector with a new type U.
-	template< typename U >
+/*	template< typename U >
 	CUDA_CALLABLE operator vec3<U>(){
-		vec3<U> result;
-		result.ptr[0] = (U)ptr[0];
-		result.ptr[1] = (U)ptr[1];
-		result.ptr[2] = (U)ptr[2];
+		vec3<U> result((U)ptr[0], (U)ptr[1], (U)ptr[2]);
+		//result.ptr[0] = (U)ptr[0];
+		//result.ptr[1] = (U)ptr[1];
+		//result.ptr[2] = (U)ptr[2];
 		return result;
-	}
+	}*/
 	// computes the squared Euclidean length (useful for several operations where only >, =, or < matter)
 	CUDA_CALLABLE T len_sq() const{
@@ -68,7 +73,22 @@ public:
 	CUDA_CALLABLE T len() const{
 		return sqrt(len_sq());
 	}
-	
+
+	/// Calculate the L2 norm of the vector, accounting for the possibility that it is complex
+	CUDA_CALLABLE T norm2() const{
+		T conj_dot = std::real(ptr[0] * std::conj(ptr[0]) + ptr[1] * std::conj(ptr[1]) + ptr[2] * std::conj(ptr[2]));
+		return sqrt(conj_dot);
+	}
+
+	/// Calculate the normalized direction vector
+	CUDA_CALLABLE vec3<T> direction() const {
+		vec3<T> result;
+		T length = norm2();
+		result[0] = ptr[0] / length;
+		result[1] = ptr[1] / length;
+		result[2] = ptr[2] / length;
+		return result;
+	}
 	/// Convert the vector from cartesian to spherical coordinates (x, y, z -> r, theta, phi where theta = [-PI, PI])
 	CUDA_CALLABLE vec3<T> cart2sph() const{
@@ -268,8 +288,36 @@ std::string str() const{
 	size_t size(){ return 3; }
 	};						//end class vec3
+
+	/// Start Class cvec3 (complex vector class)
+	template<typename T>
+	class cvec3 : public vec3< std::complex<T> > {
+
+	public:
+		CUDA_CALLABLE cvec3() : vec3< std::complex<T> >() {}
+		CUDA_CALLABLE cvec3(T x, T y, T z) : vec3< std::complex<T> >(std::complex<T>(x), std::complex<T>(y), std::complex<T>(z)) {}
+		CUDA_CALLABLE cvec3(std::complex<T> x, std::complex<T> y, std::complex<T> z) : vec3< std::complex<T> >(x, y, z) {}
+
+		CUDA_CALLABLE cvec3& operator=(const vec3< std::complex<T> > rhs) {
+			vec3< std::complex<T> >::ptr[0] = rhs.get(0);
+			vec3< std::complex<T> >::ptr[1] = rhs.get(1);
+			vec3< std::complex<T> >::ptr[2] = rhs.get(2);
+			return *this;
+		}
+		CUDA_CALLABLE std::complex<T> dot(const vec3<T> rhs) {
+			std::complex<T> result =
+				vec3< std::complex<T> >::ptr[0] * rhs.get(0) +
+				vec3< std::complex<T> >::ptr[1] * rhs.get(1) +
+				vec3< std::complex<T> >::ptr[2] * rhs.get(2);
+
+			return result;
+		}
+	};						//end class cvec3
 }							//end namespace stim
+
+
+
 /// Multiply a vector by a constant when the vector is on the right hand side
 template <typename T>
 stim::vec3<T> operator*(T lhs, stim::vec3<T> rhs){
@@ -6,6 +6,7 @@
 #include <sstream>
 #include <vector>
 #include <algorithm>
+#include <complex>
 #include <stim/cuda/cudatools/callable.h>
 #include <stim/math/vec3.h>
@@ -5,15 +5,38 @@
 #include <sstream>
 #include <cmath>
-#include "../math/vector.h"
+#include "../math/vec3.h"
 #include "../math/quaternion.h"
 #include "../math/constants.h"
 #include "../math/plane.h"
-#include "../math/complex.h"
+#include <complex>
+
+/// Calculate whether or not a vector v intersects the front (1) or back (-1) of a plane.
+/// This function returns -1 if the vector lies within the plane (is orthogonal to the surface normal)
+/*template <typename T>
+int plane_face(stim::vec3<T> v, stim::vec3<T> plane_normal) {
+	T dotprod = v.dot(plane_normal);						//calculate the dot product
+
+	if (dotprod < 0) return 1;
+	if (dotprod > 0) return -1;
+	return 0;
+}
+
+/// Calculate the component of a vector v that is perpendicular to a plane defined by a normal.
+template <typename T>
+stim::vec3<T> plane_perpendicular(stim::vec3<T> v, stim::vec3<T> plane_normal) {
+	return plane_normal * v.dot(plane_normal);
+}
+
+template <typename T>
+stim::vec3<T> plane_parallel(stim::vec3<T> v, stim::vec3<T> plane_normal) {
+	return v - plane_perpendicular(v, plane_normal);
+}*/
 namespace stim{
 	namespace optics{
+
 		/// evaluate the scalar field produced by a plane wave at a point (x, y, z)
 		/// @param x is the x-coordinate of the point
@@ -24,9 +47,9 @@ namespace stim{
 		/// @param ky is the k-vector component in the y direction
 		/// @param kz is the k-vector component in the z direction
 		template<typename T>
-		stim::complex<T> planewave_scalar(T x, T y, T z, stim::complex<T> A, T kx, T ky, T kz){
+		std::complex<T> planewave_scalar(T x, T y, T z, std::complex<T> A, T kx, T ky, T kz){
 			T d = x * kx + y * ky + z * kz;						//calculate the dot product between k and p = (x, y, z) to find the distance p is along the propagation direction
-			stim::complex<T> di = stim::complex<T>(0, d);		//calculate the phase shift that will have to be applied to propagate the wave distance d
+			std::complex<T> di = std::complex<T>(0, d);		//calculate the phase shift that will have to be applied to propagate the wave distance d
 			return A * exp(di);									//multiply the phase term by the amplitude at (0, 0, 0) to propagate the wave to p
 		}
@@ -42,7 +65,7 @@ namespace stim{
 		/// @param ky is the k-vector component in the y direction
 		/// @param kz is the k-vector component in the z direction
 		template<typename T>
-		void cpu_planewave_scalar(stim::complex<T>* field, size_t N, T* x, T* y = NULL, T* z = NULL, stim::complex<T> A = 1.0, T kx = 0.0, T ky = 0.0, T kz = 0.0){
+		void cpu_planewave_scalar(std::complex<T>* field, size_t N, T* x, T* y = NULL, T* z = NULL, std::complex<T> A = 1.0, T kx = 0.0, T ky = 0.0, T kz = 0.0){
 			T px, py, pz;
 			for(size_t i = 0; i < N; i++){										// for each element in the array
 				(x == NULL) ? px = 0 : px = x[i];								// test for NULL values
@@ -53,19 +76,137 @@ namespace stim{
 			}
 		}
+
+/*template<typename T>
+class cvec3 {
+public:
+	std::complex<T> x;
+	std::complex<T> y;
+	std::complex<T> z;
+
+	cvec3(std::complex<T> _x, std::complex<T> _y, std::complex<T> _z) {
+		x = _x;
+		y = _y;
+		z = _z;
+	}
+};*/
+
+/*template<typename T>
+class cvec3 {
+public:
+	std::complex<T> m_v[3];
+
+	cvec3(std::complex<T> x, std::complex<T> y, std::complex<T> z) {
+		m_v[0] = x;
+		m_v[1] = y;
+		m_v[2] = z;
+	}
+
+	cvec3() : cvec3(0, 0, 0) {}
+
+	void operator()(std::complex<T> x, std::complex<T> y, std::complex<T> z) {
+		m_v[0] = x;
+		m_v[1] = y;
+		m_v[2] = z;
+	}
+
+	std::complex<T> operator[](int i) const { return m_v[i]; }
+	std::complex<T>& operator[](int i) { return m_v[i]; }
+
+	/// Calculate the 2-norm of the complex vector
+	T norm2() {
+		T xx = std::real(m_v[0] * std::conj(m_v[0]));
+		T yy = std::real(m_v[1] * std::conj(m_v[1]));
+		T zz = std::real(m_v[2] * std::conj(m_v[2]));
+		return std::sqrt(xx + yy + zz);
+	}
+
+	/// Returns the normalized direction vector
+	cvec3 direction() {
+		cvec3 result;
+		
+		std::complex<T> length = norm2();
+		result[0] = m_v[0] / length;
+		result[1] = m_v[1] / length;
+		result[2] = m_v[2] / length;
+
+		return result;
+	}
+
+	std::string str() {
+		std::stringstream ss;
+		ss << m_v[0] << ", " << m_v[1] << ", " << m_v[2];
+		return ss.str();
+	}
+
+	//copy constructor
+	cvec3(const cvec3 &other) {
+		m_v[0] = other.m_v[0];
+		m_v[1] = other.m_v[1];
+		m_v[2] = other.m_v[2];
+	}
+
+	/// Assignment operator
+	cvec3& operator=(const cvec3 &rhs) {
+		m_v[0] = rhs.m_v[0];
+		m_v[1] = rhs.m_v[1];
+		m_v[2] = rhs.m_v[2];
+
+		return *this;
+	}
+
+	/// Calculate and return the cross product between this vector and another
+	cvec3 cross(cvec3 rhs) {
+		cvec3 result;
+
+		//compute the cross product (only valid for 3D vectors)
+		result[0] = (m_v[1] * rhs[2] - m_v[2] * rhs[1]);
+		result[1] = (m_v[2] * rhs[0] - m_v[0] * rhs[2]);
+		result[2] = (m_v[1] * rhs[1] - m_v[1] * rhs[0]);
+
+		return result;
+	}
+
+	/// Calculate and return the dot product between this vector and another
+	std::complex<T> dot(cvec3 rhs) {
+		return m_v[0] * rhs[0] + m_v[1] * rhs[1] + m_v[2] * rhs[2];
+	}
+
+	/// Arithmetic multiplication: returns the vector scaled by a constant value
+	template<typename R>
+	cvec3 operator*(R rhs) const
+	{
+		cvec3 result;
+		result[0] = m_v[0] * rhs;
+		result[1] = m_v[1] * rhs;
+		result[2] = m_v[2] * rhs;
+
+		return result;
+	}
+
+};
+template <typename T>
+std::ostream& operator<<(std::ostream& os, stim::optics::cvec3<T> p) {
+	os << p.str();
+	return os;
+}
+*/
+
 template<typename T>
 class planewave{
 protected:
-	stim::vec<T> k;							//k-vector, pointed in propagation direction with magnitude |k| = tau / lambda = 2pi / lambda
-	stim::vec< stim::complex<T> > E0;		//amplitude (for a scalar plane wave, only E0[0] is used)
+	cvec3<T> m_k;					//k-vector, pointed in propagation direction with magnitude |k| = tau / lambda = 2pi / lambda
+	cvec3<T> m_E;					//amplitude (for a scalar plane wave, only E0[0] is used)
 	/// Bend a plane wave via refraction, given that the new propagation direction is known
-	CUDA_CALLABLE planewave<T> bend(stim::vec<T> kn) const{
+	CUDA_CALLABLE planewave<T> bend(stim::vec3<T> v) const {
-		stim::vec<T> kn_hat = kn.norm();				//normalize the new k
-		stim::vec<T> k_hat = k.norm();					//normalize the current k
+		stim::vec3<T> k_real(m_k.get(0).real(), m_k.get(1).real(), m_k.get(2).real());			//force the vector to be real (can only refract real directions)
+
+		stim::vec3<T> kn_hat = v.direction();				//normalize the new k
+		stim::vec3<T> k_hat = k_real.direction();			//normalize the current k
 		planewave<T> new_p;								//create a new plane wave
@@ -75,51 +216,136 @@ protected:
 		if(k_dot_kn < 0) k_hat = -k_hat;				//flip k_hat
 		if(k_dot_kn == -1){
-			new_p.k = -k;
-			new_p.E0 = E0;
+			new_p.m_k = -m_k;
+			new_p.m_E = m_E;
 			return new_p;
 		}
 		else if(k_dot_kn == 1){
-			new_p.k = k;
-			new_p.E0 = E0;
+			new_p.m_k = m_k;
+			new_p.m_E = m_E;
 			return new_p;
 		}
-		vec<T> r = k_hat.cross(kn_hat);					//compute the rotation vector
+		vec3<T> r = k_hat.cross(kn_hat);					//compute the rotation vector
 		T theta = asin(r.len());						//compute the angle of the rotation about r
 		quaternion<T> q;								//create a quaternion to capture the rotation
-		q.CreateRotation(theta, r.norm());		
-		vec< complex<T> > E0n = q.toMatrix3() * E0;		//apply the rotation to E0
-		new_p.k = kn_hat * kmag();
-		new_p.E0 = E0n;
+		q.CreateRotation(theta, r.direction());	
+		stim::matrix_sq<T, 3> R = q.toMatrix3();
+		vec3< std::complex<T> > E(m_E.get(0), m_E.get(1), m_E.get(2));
+		vec3< std::complex<T> > E0n = R * E;					//apply the rotation to E0
+		//new_p.m_k = kn_hat * kmag();
+		//new_p.m_E = E0n;
+		new_p.m_k[0] = kn_hat[0] * kmag();
+		new_p.m_k[1] = kn_hat[1] * kmag();
+		new_p.m_k[2] = kn_hat[2] * kmag();
+
+		new_p.m_E[0] = E0n[0];
+		new_p.m_E[1] = E0n[1];
+		new_p.m_E[2] = E0n[2];
+
 		return new_p;
 	}
 public:
+	
+	
 	///constructor: create a plane wave propagating along k
-	CUDA_CALLABLE planewave(vec<T> kvec = stim::vec<T>(0, 0, stim::TAU), 
-							vec< complex<T> > E = stim::vec<T>(1, 0, 0))
-	{
-		//phi = phase;
+	//CUDA_CALLABLE planewave(vec<T> kvec = stim::vec<T>(0, 0, stim::TAU), 
+	//						vec< complex<T> > E = stim::vec<T>(1, 0, 0))
+	CUDA_CALLABLE planewave(std::complex<T> kx, std::complex<T> ky, std::complex<T> kz,
+		std::complex<T> Ex, std::complex<T> Ey, std::complex<T> Ez) {
+
+		m_k = cvec3<T>(kx, ky, kz);
+		m_E = cvec3<T>(Ex, Ey, Ez);
+		force_orthogonal();
+	}
-		k = kvec;
-		vec< complex<T> > k_hat = k.norm();
+	CUDA_CALLABLE planewave() : planewave(0, 0, 1, 1, 0, 0) {}
-		if(E.len() == 0)			//if the plane wave has an amplitude of 0
-			E0 = vec<T>(0);			//just return it
-		else{
-			vec< complex<T> > s = (k_hat.cross(E)).norm();		//compute an orthogonal side vector
-			vec< complex<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
-		}
+	//copy constructor
+	CUDA_CALLABLE planewave(const planewave& other) {
+		m_k = other.m_k;
+		m_E = other.m_E;
+	}
+
+	/// Assignment operator
+	CUDA_CALLABLE planewave& operator=(const planewave& rhs) {
+		m_k = rhs.m_k;
+		m_E = rhs.m_E;
-		E0 = E0 * exp( complex<T>(0, phase) );
+		return *this;
+	}
+
+	/// Forces the k and E vectors to be orthogonal
+	CUDA_CALLABLE void force_orthogonal() {
+
+		/*if (m_E.norm2() == 0) return;
+
+		cvec3<T> k_dir = m_k.direction();							//calculate the normalized direction vectors for k and E
+		cvec3<T> E_dir = m_E.direction();
+		cvec3<T> side = k_dir.cross(E_dir);						//calculate a side vector for projection
+		cvec3<T> side_dir = side.direction();					//normalize the side vector
+		E_dir = side_dir.cross(k_dir);								//calculate the new E vector direction
+		T E_norm = m_E.norm2();
+		m_E = E_dir * E_norm;								//apply the new direction to the existing E vector
+		*/
+	}
+
+	CUDA_CALLABLE cvec3<T> k() {
+		return m_k;
+	}
+
+	CUDA_CALLABLE cvec3<T> E() {
+		return m_E;
+	}
+
+	CUDA_CALLABLE cvec3<T> evaluate(T x, T y, T z) {
+		
+		std::complex<T> k_dot_r = m_k[0] * x + m_k[1] * y + m_k[2] * z;
+		std::complex<T> e_k_dot_r = std::exp(std::complex<T>(0, 1) * k_dot_r);
+
+		cvec3<T> result;
+		result[0] = m_E[0] * e_k_dot_r;
+		result[1] = m_E[1] * e_k_dot_r;
+		result[2] = m_E[2] * e_k_dot_r;
+		return result;
+	}
+
+	/*int scatter(vec3<T> surface_normal, vec3<T> surface_point, T ni, std::complex<T> nt,
+		planewave& Pi, planewave& Pr, planewave& Pt) {
+
+		
+		return 0;
+
+	}*/
+
+	CUDA_CALLABLE T kmag() const {
+		return std::sqrt(std::real(m_k.get(0) * std::conj(m_k.get(0)) + m_k.get(1) * std::conj(m_k.get(1)) + m_k.get(2) * std::conj(m_k.get(2))));
+	}
+
+	CUDA_CALLABLE planewave<T> refract(stim::vec3<T> kn) const {
+		return bend(kn);
+	}
+
+	/// Return a plane wave with the origin translated by (x, y, z)
+	CUDA_CALLABLE planewave<T> translate(T x, T y, T z) const {
+		planewave<T> result;
+		cvec3<T> k = m_k;
+		result.m_k = k;
+		std::complex<T> k_dot_r = k[0] * (-x) + k[1] * (-y) + k[2] * (-z);
+		std::complex<T> exp_k_dot_r = std::exp(std::complex<T>(0.0, 1.0) * k_dot_r);
+
+		cvec3<T> E = m_E;
+		result.m_E[0] = E[0] * exp_k_dot_r;
+		result.m_E[1] = E[1] * exp_k_dot_r;
+		result.m_E[2] = E[2] * exp_k_dot_r;
+		return result;
 	}
 	///multiplication operator: scale E0
-    CUDA_CALLABLE planewave<T> & operator* (const T & rhs){		
+    /*CUDA_CALLABLE planewave<T>& operator* (const T& rhs) {
 		E0 = E0 * rhs;
 		return *this;
 	}
@@ -177,7 +403,7 @@ public:
 	/// @param t is the transmitted component of the plane wave
 	void scatter(stim::plane<T> P, T n0, T n1, planewave<T> &r, planewave<T> &t){
 		scatter(P, n1/n0, r, t);
-	}
+	}*/
 	/// Calculate the scattering result when nr = n1/n0
@@ -186,41 +412,63 @@ public:
 	/// @param n1 is the refractive index inside the surface (in the direction away from the normal)
 	/// @param r is the reflected component of the plane wave
 	/// @param t is the transmitted component of the plane wave
-	void scatter(stim::plane<T> P, T nr, planewave<T> &r, planewave<T> &t){
+	
+	/*void scatter(vec3<T> plane_normal, vec3<T> plane_position, std::complex<T> nr, planewave<T>& r, planewave<T>& t) {
+
+		//TODO: Generate a real vector from the current K vector - this will not address complex k-vectors
+		vec3< std::complex<T> > k(3);
+		k[0] = m_k[0];
+		k[1] = m_k[1];
+		k[2] = m_k[2];
+
+		vec3< std::complex<T> > E(3);
+		E[0] = m_E[0];// std::complex<T>(m_E[0].real(), m_E[0].imag());
+		E[1] = m_E[1];// std::complex<T>(m_E[1].real(), m_E[1].imag());
+		E[2] = m_E[2];// std::complex<T>(m_E[2].real(), m_E[2].imag());
+
+		std::complex<T> cOne(1, 0);
+		std::complex<T> cTwo(2, 0);
+
+		T kmag = m_k.norm2();
-		int facing = P.face(k);		//determine which direction the plane wave is coming in
+		int facing = plane_face(k, plane_normal);		//determine which direction the plane wave is coming in
-		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)
+		if(facing == -1){								//if the wave hits the back of the plane, invert the plane and nr
+			plane_normal = -plane_normal;				//flip the plane
+			nr = cOne / nr;									//invert the refractive index (now nr = n0/n1)
 		}
 		//use Snell's Law to calculate the transmitted angle
-		T cos_theta_i = k.norm().dot(-P.norm());				//compute the cosine of theta_i
+		//vec3<T> plane_normal = -P.n();
+		T cos_theta_i = k.norm().dot(-plane_normal);				//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
+		std::complex<T> sin_theta_t = (cOne / nr) * sin(theta_i);						//compute the sine of theta_t using Snell's law
+		std::complex<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 = stim::PI / (T)2;
-		}
+		//bool tir = false;						//flag for total internal reflection
+		//if(theta_t != theta_t){
+		//	tir = true;
+		//	theta_t = stim::PI / (T)2;
+		//}
 		//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< complex<T> > Er = E0 * rp;		//compute the E vectors
-			vec< complex<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::complex<T> rp = (cOne - nr) / (cOne + nr);		//compute the Fresnel coefficients
+			std::complex<T> tp = (cOne * cTwo) / (cOne + nr);
+			vec3<T> kr = -k;
+			vec3<T> kt = k * nr;			//set the k vectors for theta_i = 0
+			vec3< std::complex<T> > Er = E * rp;		//compute the E vectors
+			vec3< std::complex<T> > Et = E * tp;
+			T phase_t = plane_position.dot(k - kt);	//compute the phase offset
+			T phase_r = plane_position.dot(k - kr);
 			//create the plane waves
-			r = planewave<T>(kr, Er, phase_r);
-			t = planewave<T>(kt, Et, phase_t);
+			//r = planewave<T>(kr, Er, phase_r);
+			//t = planewave<T>(kt, Et, phase_t);
+
+			r = planewave(kr[0], kr[1], kr[2], Er[0], Er[1], Er[2]);
+			t = planewave(kt[0], kt[1], kt[2], Et[0], Et[1], Et[2]);
+
 			return;
 		}
@@ -230,57 +478,61 @@ public:
 		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{
+		//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);
 		}
 		//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();
+		vec3<T> z_hat = -plane_normal;
+		vec3<T> y_hat = plane_parallel(k, plane_normal).norm();
+		vec3<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;
+		vec3<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
-		complex<T> Ei_s = E0.dot(x_hat);
-		int sgn = E0.dot(y_hat).sgn();
-		vec< complex<T> > cx_hat = x_hat;
-		complex<T> Ei_p = ( E0 - cx_hat * Ei_s ).len() * sgn;
+		std::complex<T> Ei_s = E.dot(x_hat);
+		//int sign = sgn(E0.dot(y_hat));
+		vec3< std::complex<T> > cx_hat = x_hat;
+		std::complex<T> Ei_p = (E - cx_hat * Ei_s).len();// *(T)sign;
 		//compute the magnitude of the p- and s-polarized components of the reflected E vector
-		complex<T> Er_s = Ei_s * rs;
-		complex<T> Er_p = Ei_p * rp;
+		std::complex<T> Er_s = Ei_s * rs;
+		std::complex<T> Er_p = Ei_p * rp;
 		//compute the magnitude of the p- and s-polarized components of the transmitted E vector
-		complex<T> Et_s = Ei_s * ts;
-		complex<T> Et_p = Ei_p * tp;
+		std::complex<T> Et_s = Ei_s * ts;
+		std::complex<T> Et_p = Ei_p * tp;
 		//compute the reflected E vector
-		vec< complex<T> > Er = vec< complex<T> >(y_hat * cos(theta_i) + z_hat * sin(theta_i)) * Er_p + cx_hat * Er_s;
+		vec3< std::complex<T> > Er = vec3< std::complex<T> >(y_hat * cos(theta_i) + z_hat * sin(theta_i)) * Er_p + cx_hat * Er_s;
 		//compute the transmitted E vector
-		vec< complex<T> > Et = vec< complex<T> >(y_hat * cos(theta_t) - z_hat * sin(theta_t)) * Et_p + cx_hat * Et_s;
+		vec3< std::complex<T> > Et = vec3< std::complex<T> >(y_hat * cos(theta_t) - z_hat * sin(theta_t)) * Et_p + cx_hat * Et_s;
-		T phase_t = P.p().dot(k - kt);
-		T phase_r = P.p().dot(k - kr);
+		T phase_t = plane_position.dot(k - kt);
+		T phase_r = plane_position.dot(k - kr);
 		//create the plane waves
-		r.k = kr;
-		r.E0 = Er * exp( complex<T>(0, phase_r) );
+		//r.m_k = kr;
+		//r.m_E = Er * exp( complex<T>(0, phase_r) );
-		t.k = kt;
-		t.E0 = Et * exp( complex<T>(0, phase_t) );
-	}
+		//t.m_k = kt;
+		//t.m_E = Et * exp( complex<T>(0, phase_t) );
+
+		r = planewave(kr[0], kr[1], kr[2], Er[0], Er[1], Er[2]);
+		t = planewave(kt[0], kt[1], kt[2], Et[0], Et[1], Et[2]);
+	}*/
 	std::string str()
 	{
 		std::stringstream ss;
-		ss<<"Plane Wave:"<<std::endl;
-		ss<<"	"<<E0<<" e^i ( "<<k<<" . r )";
+		ss << "k: " << m_k << std::endl;
+		ss << "E: " << m_E << std::endl;
 		return ss.str();
 	}
 };					//end planewave class