diff --git a/stim/math/matrix_sq.h b/stim/math/matrix_sq.h
index 3d56b99..3c642d0 100644
--- a/stim/math/matrix_sq.h
+++ b/stim/math/matrix_sq.h
@@ -8,7 +8,7 @@
 #include <stim/math/vec3.h>
 #include <stim/cuda/cudatools/callable.h>
 
-namespace stim{
+namespace tira {
 
 template <class T, int N>
 struct matrix_sq
@@ -100,7 +100,7 @@ struct matrix_sq
 
 	template<typename Y>
 	CUDA_CALLABLE vec3<Y> operator*(vec3<Y> rhs){
-		vec3<Y> result;
+		vec3<Y> result(0, 0, 0);
 		for(int r=0; r<3; r++)
 			for(int c=0; c<3; c++)
 				result[r] += (*this)(r, c) * rhs[c];
@@ -152,7 +152,7 @@ struct matrix_sq
 }	//end namespace rts
 
 template <typename T, int N>
-std::ostream& operator<<(std::ostream& os, stim::matrix_sq<T, N> M)
+std::ostream& operator<<(std::ostream& os, tira::matrix_sq<T, N> M)
 {
     os<<M.toStr();
     return os;
diff --git a/stim/math/plane.h b/stim/math/plane.h
index 6dd54f6..88ef37a 100644
--- a/stim/math/plane.h
+++ b/stim/math/plane.h
@@ -1,5 +1,5 @@
-#ifndef STIM_PLANE_H
-#define STIM_PLANE_H
+#ifndef TIRA_PLANE_H
+#define TIRA_PLANE_H
 
 #include <iostream>
 #include <stim/math/vec3.h>
@@ -7,32 +7,32 @@
 #include <stim/math/quaternion.h>
 
 
-namespace stim
+namespace tira
 {
 template<typename T> class plane;
 }
 
 template<typename T>
-CUDA_CALLABLE stim::plane<T> operator-(stim::plane<T> v);
+CUDA_CALLABLE tira::plane<T> operator-(tira::plane<T> v);
 
-namespace stim
+namespace tira
 {
 
 template <typename T>
 class plane
 {
 	protected:
-		stim::vec3<T> P;
-		stim::vec3<T> N;
-		stim::vec3<T> U;
+		vec3<T> P;
+		vec3<T> N;
+		vec3<T> U;
 
 		///Initializes the plane with standard coordinates.
 		///
 		CUDA_CALLABLE void init()
 		{
-			P = stim::vec3<T>(0, 0, 0);
-			N = stim::vec3<T>(0, 0, 1);
-			U = stim::vec3<T>(1, 0, 0);
+			P = vec3<T>(0, 0, 0);
+			N = vec3<T>(0, 0, 1);
+			U = vec3<T>(1, 0, 0);
 		}
 
 	public:
@@ -60,7 +60,7 @@ class plane
 		{
 			init();
 			P = c;
-			stim::vec3<T> n = (c - a).cross(b - a);
+			vec3<T> n = (c - a).cross(b - a);
 			try
 			{
 				if(n.len() != 0)
@@ -165,9 +165,9 @@ class plane
 
 		}
 
-		CUDA_CALLABLE stim::plane<T> operator-()
+		CUDA_CALLABLE plane<T> operator-()
 		{
-			stim::plane<T> p = *this;
+			plane<T> p = *this;
 
 			//negate the normal vector
 			p.N = -p.N;
diff --git a/stim/math/quaternion.h b/stim/math/quaternion.h
index bbb6cb3..65b0fb8 100644
--- a/stim/math/quaternion.h
+++ b/stim/math/quaternion.h
@@ -1,10 +1,10 @@
-#ifndef RTS_QUATERNION_H
-#define RTS_QUATERNION_H
+#ifndef TIRA_QUATERNION_H
+#define TIRA_QUATERNION_H
 
 #include <stim/math/matrix_sq.h>
 #include <stim/cuda/cudatools/callable.h>
 
-namespace stim{
+namespace tira {
 
 template<typename T>
 class quaternion
diff --git a/stim/math/vec3.h b/stim/math/vec3.h
index e8f71d5..1a52ea0 100644
--- a/stim/math/vec3.h
+++ b/stim/math/vec3.h
@@ -1,5 +1,5 @@
-#ifndef STIM_VEC3_H
-#define STIM_VEC3_H
+#ifndef TIRA_VEC3_H
+#define TIRA_VEC3_H
 
 
 #include <stim/cuda/cudatools/callable.h>
@@ -9,7 +9,7 @@
 #include <sstream>
 
 
-namespace stim{
+namespace tira{
 
 
 /// A class designed to act as a 3D vector with CUDA compatibility
@@ -313,20 +313,20 @@ std::string str() const{
 			return result;
 		}
 	};						//end class cvec3
-}							//end namespace stim
+}							//end namespace tira
 
 
 
 
 /// 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){
+tira::vec3<T> operator*(T lhs, tira::vec3<T> rhs){
     return rhs * lhs;
 }
 
 //stream operator
 template<typename T>
-std::ostream& operator<<(std::ostream& os, stim::vec3<T> const& rhs){
+std::ostream& operator<<(std::ostream& os, tira::vec3<T> const& rhs){
 	os<<rhs.str();
 	return os;
 }
diff --git a/stim/math/vector.h b/stim/math/vector.h
index 0af7cf8..e9c6035 100644
--- a/stim/math/vector.h
+++ b/stim/math/vector.h
@@ -1,5 +1,5 @@
-#ifndef STIM_VECTOR_H
-#define STIM_VECTOR_H
+#ifndef TIRA_VECTOR_H
+#define TIRA_VECTOR_H
 
 #include <iostream>
 #include <cmath>
@@ -11,7 +11,7 @@
 #include <stim/cuda/cudatools/callable.h>
 #include <stim/math/vec3.h>
 
-namespace stim
+namespace tira
 {
 
 template <class T>
@@ -338,8 +338,8 @@ struct vec : public std::vector<T>
 	}
 
 	/// Cast to a vec3
-	operator stim::vec3<T>(){
-		stim::vec3<T> r;
+	operator tira::vec3<T>(){
+		tira::vec3<T> r;
 		size_t N = std::min(size(), (size_t)3);
 		for(size_t i = 0; i < N; i++)
 			r[i] = at(i);
@@ -405,10 +405,10 @@ struct vec : public std::vector<T>
 };
 
 
-}	//end namespace rts
+}	//end namespace tira
 
 template <typename T>
-std::ostream& operator<<(std::ostream& os, stim::vec<T> v)
+std::ostream& operator<<(std::ostream& os, tira::vec<T> v)
 {
     os<<v.str();
     return os;
@@ -417,9 +417,9 @@ std::ostream& operator<<(std::ostream& os, stim::vec<T> v)
 
 /// Multiply a vector by a constant when the vector is on the right hand side
 template <typename T>
-stim::vec<T> operator*(T lhs, stim::vec<T> rhs)
+tira::vec<T> operator*(T lhs, tira::vec<T> rhs)
 {
-	stim::vec<T> r;
+	tira::vec<T> r;
 
     return rhs * lhs;
 }
diff --git a/stim/optics/beam.h b/stim/optics/beam.h
new file mode 100644
index 0000000..b82c158
--- /dev/null
+++ b/stim/optics/beam.h
@@ -0,0 +1,98 @@
+#ifndef TIRA_BEAM
+#define TIRA_BEAM
+
+#include "../math/vec3.h"
+#include "../math/function.h"
+#include "../optics/planewave.h"
+#include <vector>
+#include <numbers>
+#include <random>
+
+namespace tira{
+	namespace optics{
+template<typename T>
+class beam {
+
+	vec3<T> m_direction;
+	vec3<T> m_focus;
+	cvec3<T> m_E;
+	T m_k;
+
+
+
+public:
+	/// <summary>
+	/// Beam constructor
+	/// </summary>
+	/// <param name="d">Propagation direction of the beam (automatically normalized).</param>
+	/// <param name="f">Focal point for the beam.</param>
+	/// <param name="pol">Beam polarization (automatically normalized and orthogonalized). Only works for linear polarization, though.</param>
+	/// <param name="E">Complex beam amplitude (applied to the polarization to create a linearly polarized beam).</param>
+	/// <param name="lambda">Wavelength</param>
+	beam(
+		vec3<T> d = vec3<T>(0, 0, 1), 
+		vec3<T> f = vec3<T>(0, 0, 0),
+		vec3<T> pol = vec3<T>(1, 0, 0),
+		std::complex<T> E = std::complex<T>(1.0, 0.0),
+		T lambda = 1.0
+	) {
+
+		m_direction = d.direction();				//normalize and set the plane wave direction
+		pol = (pol - d.dot(pol)).direction();			//orthogonalize and normalize the polarization vector
+		m_focus = f;								//set the focal point
+		m_E[0] = pol[0] * E;						//convert the polarization and complex amplitude to a linearly-polarized E vector
+		m_E[1] = pol[1] * E;
+		m_E[2] = pol[2] * E;
+		m_k = 2.0 * std::numbers::pi * lambda;		//calculate and store the wavenumber
+	}
+
+	/// <summary>
+	/// Generate a focal spot for the beam using Monte-Carlo sampling.
+	/// </summary>
+	/// <param name="NA">Numerical aperture of the focusing optics.</param>
+	/// <param name="samples">Number of Monte-Carlo samples.</param>
+	/// <param name="NA_in">NA of the internal obscuration.</param>
+	/// <returns></returns>
+	std::vector< planewave<T> > mcfocus(T NA, size_t samples, T NA_in = 0) {
+
+		//create a rotation matrix to orient the beam along the specified direction vector
+		tira::matrix_sq<double, 3> R;// = tira::matrix_sq<double, 3>::identity();		//initialize the rotation matrix to the identity matrix
+		tira::vec3<double> Z(0, 0, 1);												//create a vector along the Z direction
+		tira::quaternion<double> q;
+		q.CreateRotation(Z, m_direction);
+		R = q.toMatrix3();
+
+		std::vector< planewave<T> > result(samples);								//generate an array of N plane waves
+		double alpha = std::asin(NA);												//calculate the angle of the band-pass
+		double alpha_in = std::asin(NA_in);											//calculate the angle of the central obscuration
+		double cos_alpha = std::cos(alpha);											//calculate the cosine of the band-pass and central obscuration
+		double cos_alpha_in = std::cos(alpha_in);
+
+		//random sampling is done using Archimede's method - uniform samples are taken in cylindrical space and projected onto a sphere
+		std::random_device rd;														//create a random number generator
+		std::default_random_engine eng(rd());
+		std::uniform_real_distribution<> theta_dist(0.0, 2.0 * std::numbers::pi);	//create a distribution for theta (in cylindrical coordinates)
+		std::uniform_real_distribution<> z_dist(cos_alpha, cos_alpha_in);			//create a uniform distribution for sampling the z-axis
+
+		//generate a guide plane wave
+		planewave<T> pw(m_direction[0] * m_k, m_direction[1] * m_k, m_direction[2] * m_k, m_E[0], m_E[1], m_E[2]);
+
+		double theta, phi;															//variables to store the coordinates for theta/phi on a sphere
+		tira::vec3<T> dir;
+		T w = 1.0 / samples;														//integration weight
+		for (size_t n = 0; n < samples; n++) {
+			theta = theta_dist(eng);												//randomly generate a theta coordinate between 0 and 2pi
+			phi = std::acos(z_dist(eng));											//randomly generate a z coordinate based on the NA and project to phi
+			dir = R * tira::vec3<double>(1.0, theta, phi).sph2cart();				//convert the value to Cartesian coordinates and store the sample vector
+			result[n] = pw.refract(dir);											//refract the plane wave to align with the sample direction
+			result[n] = result[n].translate(m_focus[0], m_focus[1], m_focus[2]);	//refocus the plane wave to the focal position
+			result[n] = result[n] * w;												//apply the integration weight (1/N)
+		}
+
+		return result;																//return the sample array
+	}
+};
+}								//end namespace optics
+}								//end namespace tira
+
+#endif
diff --git a/stim/optics/planewave.h b/stim/optics/planewave.h
index e63354a..1832458 100644
--- a/stim/optics/planewave.h
+++ b/stim/optics/planewave.h
@@ -1,5 +1,5 @@
-#ifndef STIM_PLANEWAVE_H
-#define STIM_PLANEWAVE_H
+#ifndef TIRA_PLANEWAVE_H
+#define TIRA_PLANEWAVE_H
 
 #include <string>
 #include <sstream>
@@ -11,29 +11,8 @@
 #include "../math/plane.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 tira{
 	namespace optics{
 
 
@@ -77,121 +56,6 @@ 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{
 
@@ -201,12 +65,12 @@ protected:
 	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::vec3<T> v) const {
+	CUDA_CALLABLE planewave<T> bend(vec3<T> v) const {
 
-		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)
+		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
+		vec3<T> kn_hat = v.direction();					//normalize the new k
+		vec3<T> k_hat = k_real.direction();				//normalize the current k
 
 		planewave<T> new_p;								//create a new plane wave
 
@@ -227,10 +91,10 @@ protected:
 		}
 
 		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
+		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.direction());	
-		stim::matrix_sq<T, 3> R = q.toMatrix3();
+		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();
@@ -252,8 +116,6 @@ 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))
 	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) {
 
@@ -313,22 +175,10 @@ public:
 		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;
@@ -350,6 +200,17 @@ public:
 		return *this;
 	}
 
+	///return a plane wave with the applied refractive index (scales the k-vector by the input)
+	CUDA_CALLABLE planewave<T> ri(T n) {
+		planewave<T> result;
+		result.m_E = m_E;
+		result.m_k = m_k * n;
+		return result;
+	}
+	CUDA_CALLABLE planewave<T> refract(vec3<T> kn) const {
+		return bend(kn);
+	}
+
 	/*CUDA_CALLABLE T lambda() const{
 		return stim::TAU / k.len();
 	}
@@ -390,9 +251,7 @@ public:
 		return planewave<T>(k * (nt / ni), E0);
 	}
 
-	CUDA_CALLABLE planewave<T> refract(stim::vec<T> kn) const{
-		return bend(kn);
-	}
+	
 
 	/// Calculate the result of a plane wave hitting an interface between two refractive indices
 
@@ -408,29 +267,28 @@ public:
 	/// Calculate the scattering result when nr = n1/n0
 
 	/// @param P is a plane representing the position and orientation of the surface
-	/// @param r is the ration n1/n0
+	/// @param r is the ratio n1/n0
 	/// @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(vec3<T> plane_normal, vec3<T> plane_position, std::complex<T> nr, planewave<T>& r, planewave<T>& t) {
+	int scatter(vec3<T> plane_normal, vec3<T> plane_position, std::complex<T> nr, planewave<T>& r, planewave<T>& t) {
 		
 		if (m_k[0].imag() != 0.0 || m_k[1].imag() != 0.0 || m_k[2].imag() != 0) {
 			std::cout << "ERROR: cannot scatter a plane wave with an imaginary k-vector." << std::endl;
 		}
 
-		stim::vec3<T> ki(m_k[0].real(), m_k[1].real(), m_k[2].real());	//force the current k vector to be real
-		stim::vec3<T> kr;
-		stim::cvec3<T> kt, Ei, Er, Et;
+		vec3<T> ki(m_k[0].real(), m_k[1].real(), m_k[2].real());	//force the current k vector to be real
+		vec3<T> kr;
+		cvec3<T> kt, Ei, Er, Et;
 
 		plane_normal = plane_normal.direction();
-		stim::vec3<T> k_dir = ki.direction();				// calculate the direction of the incident plane wave
-
-		int facing = plane_face(k_dir, plane_normal);		//determine which direction the plane wave is coming in
+		vec3<T> k_dir = ki.direction();								//calculate the direction of the incident plane wave
 
-		if (facing == -1) {									//if the wave hits the back of the plane, invert the plane and nr
+		//int facing = plane_face(k_dir, plane_normal);				//determine which direction the plane wave is coming in
+		if (k_dir.dot(plane_normal) > 0) {							//if the wave hits the back of the plane, invert the plane and nr
 			std::cout << "ERROR: k-vector intersects the wrong side of the boundary." << std::endl;
-			return -1;										//the plane wave is impacting the wrong side of the surface
+			return -1;												//the plane wave is impacting the wrong side of the surface
 		}
 
 		//use Snell's Law to calculate the transmitted angle
@@ -443,17 +301,15 @@ public:
 			std::complex<T> rp = (1.0 - nr) / (1.0 + nr);			//compute the Fresnel coefficients
 			std::complex<T> tp = 2.0 / (1.0 + nr);
 
-			//ki = kv * ni;											//calculate the incident k-vector (scaled by the incident refractive index)
 			kr = -ki;												//the reflection vector is the inverse of the incident vector
 			kt[0] = ki[0] * nr;
 			kt[1] = ki[1] * nr;
 			kt[2] = ki[2] * nr;
-			//Ei = Ev;												//compute the E vectors for all three plane waves based on the Fresnel coefficients
-			Er = m_E * rp;
+			
+			Er = m_E * rp;											//compute the E vectors based on the Fresnel coefficients
 			Et = m_E * tp;
 
 			//calculate the phase offset based on the plane positions
-			//T phase_i = plane_position.dot(kv - ki);			//compute the phase offset for each plane wave
 			T phase_r = plane_position.dot(ki - kr);
 			std::complex<T> phase_t =
 				plane_position[0] * (ki[0] - kt[0]) +
@@ -465,25 +321,25 @@ public:
 			T sin_theta_r = sin_theta_i;
 			T theta_r = theta_i;
 
-			std::complex<T> sin_theta_t = nr * sin(theta_i);		//compute the sine of theta_t using Snell's law
+			std::complex<T> sin_theta_t = (1.0/nr) * sin(theta_i);		//compute the sine of theta_t using Snell's law
 			std::complex<T> cos_theta_t = std::sqrt(1.0 - sin_theta_t * sin_theta_t);
 			std::complex<T> theta_t = asin(sin_theta_t);				//compute the cosine of theta_t
 
 			//Define the basis vectors for the calculation (plane of incidence)
-			stim::vec3<T> z_hat = -plane_normal;
-			stim::vec3<T> y_hat = plane_parallel(k_dir, plane_normal);
-			stim::vec3<T> x_hat = y_hat.cross(z_hat);
+			vec3<T> z_hat = -plane_normal;
+			vec3<T> plane_perpendicular = plane_normal * k_dir.dot(plane_normal);
+			vec3<T> y_hat = (k_dir - plane_perpendicular).direction();
+			vec3<T> x_hat = y_hat.cross(z_hat);
 
 			//calculate the k-vector magnitudes
-			//T kv_mag = kv.norm2();
 			T ki_mag = ki.norm2();
 			T kr_mag = ki_mag;
 			std::complex<T> kt_mag = ki_mag * nr;
 
 			//calculate the k vector directions
-			stim::vec3<T> ki_dir = y_hat * sin_theta_i + z_hat * cos_theta_i;
-			stim::vec3<T> kr_dir = y_hat * sin_theta_r - z_hat * cos_theta_r;
-			stim::cvec3<T> kt_dir;
+			vec3<T> ki_dir = y_hat * sin_theta_i + z_hat * cos_theta_i;
+			vec3<T> kr_dir = y_hat * sin_theta_r - z_hat * cos_theta_r;
+			cvec3<T> kt_dir;
 			kt_dir[0] = y_hat[0] * sin_theta_t + z_hat[0] * cos_theta_t;
 			kt_dir[1] = y_hat[1] * sin_theta_t + z_hat[1] * cos_theta_t;
 			kt_dir[2] = y_hat[2] * sin_theta_t + z_hat[2] * cos_theta_t;
@@ -500,34 +356,22 @@ public:
 			std::complex<T> tp = ((2.0 * sin_theta_t * cos_theta_i) / (std::sin(theta_t + theta_i) * std::cos(theta_t - theta_i)));
 
 			//calculate the p component directions for each E vector
-			stim::vec3<T> Eip_dir = y_hat * cos_theta_i - z_hat * sin_theta_i;
-			stim::vec3<T> Erp_dir = y_hat * cos_theta_r + z_hat * sin_theta_r;
-			stim::cvec3<T> Etp_dir;
+			vec3<T> Eip_dir = y_hat * cos_theta_i - z_hat * sin_theta_i;
+			vec3<T> Erp_dir = y_hat * cos_theta_r + z_hat * sin_theta_r;
+			cvec3<T> Etp_dir;
 			Etp_dir[0] = y_hat[0] * cos_theta_t - z_hat[0] * sin_theta_t;
 			Etp_dir[1] = y_hat[1] * cos_theta_t - z_hat[1] * sin_theta_t;
 			Etp_dir[2] = y_hat[2] * cos_theta_t - z_hat[2] * sin_theta_t;
 
 			//calculate the s and t components of each E vector
-			//std::complex<T> E_mag = Ev.norm2();
 			std::complex<T> Ei_s = m_E.dot(x_hat);
-			//std::complex<T> Ei_p = std::sqrt(E_mag * E_mag - Ei_s * Ei_s);
 			std::complex<T> Ei_p = m_E.dot(Eip_dir);
 			std::complex<T> Er_s = rs * Ei_s;
 			std::complex<T> Er_p = rp * Ei_p;
 			std::complex<T> Et_s = ts * Ei_s;
 			std::complex<T> Et_p = tp * Ei_p;
 
-
-
 			//calculate the E vector for each plane wave
-			//Ei[0] = Eip_dir[0] * Ei_p + x_hat[0] * Ei_s;
-			//Ei[1] = Eip_dir[1] * Ei_p + x_hat[1] * Ei_s;
-			//Ei[2] = Eip_dir[2] * Ei_p + x_hat[2] * Ei_s;
-
-			//std::cout << Ev << std::endl;
-			//std::cout << Ei << std::endl;
-			//std::cout << std::endl;
-
 			Er[0] = Erp_dir[0] * Er_p + x_hat[0] * Er_s;
 			Er[1] = Erp_dir[1] * Er_p + x_hat[1] * Er_s;
 			Er[2] = Erp_dir[2] * Er_p + x_hat[2] * Er_s;
@@ -539,134 +383,21 @@ public:
 
 
 		//calculate the phase offset based on the plane positions
-		//T phase_i = plane_position.dot(kv - ki);			//compute the phase offset for each plane wave
 		T phase_r = plane_position.dot(ki - kr);
 		std::complex<T> phase_t =
 			plane_position[0] * (ki[0] - kt[0]) +
 			plane_position[1] * (ki[1] - kt[1]) +
 			plane_position[2] * (ki[2] - kt[2]);
 
-
-		//Ei = Ei * std::exp(std::complex<T>(0, 1));
+		//apply the phase offset
 		Er = Er * std::exp(std::complex<T>(0, 1) * phase_r);
 		Et = Et * std::exp(std::complex<T>(0, 1) * phase_t);
 
-		//Pi = stim::optics::planewave<T>(ki[0], ki[1], ki[2], Ei[0], Ei[1], Ei[2]);
-		r = stim::optics::planewave<T>(kr[0], kr[1], kr[2], Er[0], Er[1], Er[2]);
-		t = stim::optics::planewave<T>(kt[0], kt[1], kt[2], Et[0], Et[1], Et[2]);
+		//generate the reflected and transmitted waves
+		r = planewave<T>(kr[0], kr[1], kr[2], Er[0], Er[1], Er[2]);
+		t = planewave<T>(kt[0], kt[1], kt[2], Et[0], Et[1], Et[2]);
 
 		return 0;
-		/*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 = 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
-			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
-		//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
-		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;
-		//}
-
-		//handle the degenerate case where theta_i is 0 (the plane wave hits head-on)
-		if(theta_i == 0){
-			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(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;
-		}
-
-
-		//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);
-		}
-
-		//compute the coordinate space for the plane of incidence
-		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
-		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
-		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
-		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
-		std::complex<T> Et_s = Ei_s * ts;
-		std::complex<T> Et_p = Ei_p * tp;
-
-		//compute the reflected E vector
-		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
-		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 = plane_position.dot(k - kt);
-		T phase_r = plane_position.dot(k - kr);
-
-		//create the plane waves
-		//r.m_k = kr;
-		//r.m_E = Er * exp( complex<T>(0, phase_r) );
-
-		//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()
@@ -678,10 +409,10 @@ public:
 	}
 };					//end planewave class
 }					//end namespace optics
-}					//end namespace stim
+}					//end namespace tira
 
 template <typename T>
-std::ostream& operator<<(std::ostream& os, stim::optics::planewave<T> p)
+std::ostream& operator<<(std::ostream& os, tira::optics::planewave<T> p)
 {
     os<<p.str();
     return os;
--
libgit2 0.21.4