#ifndef STIM_PLANEWAVE_H
#define STIM_PLANEWAVE_H

#include <string>
#include <sstream>
#include <cmath>

#include "../math/vec3.h"
#include "../math/quaternion.h"
#include "../math/constants.h"
#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 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
		/// @param y is the y-coordinate of the point
		/// @param z is the z-coordinate of the point
		/// @param A is the amplitude of the plane wave, specifically the field at (0, 0, 0)
		/// @param kx is the k-vector component in the x direction
		/// @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>
		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
			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
		}

		/// evaluate the scalar field produced by a plane wave at several positions

		/// @param field is a pre-allocated block of memory that will store the complex field at all points
		/// @param N is the number of field values to be evaluated
		/// @param x is a set of x coordinates defining positions within the field (NULL implies that all values are zero)
		/// @param y is a set of y coordinates defining positions within the field (NULL implies that all values are zero)
		/// @param z is a set of z coordinates defining positions within the field (NULL implies that all values are zero)
		/// @param A is the amplitude of the plane wave, specifically the field at (0, 0, 0)
		/// @param kx is the k-vector component in the x direction
		/// @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(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
				(y == NULL) ? py = 0 : py = y[i];
				(z == NULL) ? pz = 0 : pz = z[i];

				field[i] = planewave_scalar(px, py, pz, A, kx, ky, kz);			// call the single-value plane wave function
			}
		}


/*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:

	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::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)

		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

		T k_dot_kn = k_hat.dot(kn_hat);					//if kn is equal to k or -k, handle the degenerate case

		//if k . n < 0, then the bend is a reflection
		if(k_dot_kn < 0) k_hat = -k_hat;				//flip k_hat

		if(k_dot_kn == -1){
			new_p.m_k = -m_k;
			new_p.m_E = m_E;
			return new_p;
		}
		else if(k_dot_kn == 1){
			new_p.m_k = m_k;
			new_p.m_E = m_E;
			return new_p;
		}

		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.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))
	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();
	}

	CUDA_CALLABLE planewave() : planewave(0, 0, 1, 1, 0, 0) {}

	//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;

		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) {
		m_E = m_E * rhs;
		return *this;
	}

	/*CUDA_CALLABLE T lambda() const{
		return stim::TAU / k.len();
	}

	CUDA_CALLABLE T kmag() const{
		return k.len();
	}

	CUDA_CALLABLE vec< complex<T> > E(){
		return E0;
	}

	CUDA_CALLABLE vec<T> kvec(){
		return k;
	}

	/// calculate the value of the field produced by the plane wave given a three-dimensional position
	CUDA_CALLABLE vec< complex<T> > pos(T x, T y, T z){
		return pos( stim::vec<T>(x, y, z) );
	}

	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);
	}

	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

	/// @param P is a plane representing the position and orientation of the surface
	/// @param n0 is the refractive index outside of the surface (in the direction of the normal)
	/// @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 n0, T n1, planewave<T> &r, planewave<T> &t){
		scatter(P, n1/n0, r, t);
	}*/

	/// 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 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) {
		
		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;

		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

		if (facing == -1) {									//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
		}

		//use Snell's Law to calculate the transmitted angle
		T cos_theta_i = k_dir.dot(-plane_normal);					//compute the cosine of theta_i
		T sin_theta_i = std::sqrt(1 - cos_theta_i * cos_theta_i);
		T theta_i = acos(cos_theta_i);								//compute theta_i

		//handle the degenerate case where theta_i is 0 (the plane wave hits head-on)
		if (theta_i == 0) {
			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;
			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]) +
				plane_position[1] * (ki[1] - kt[1]) +
				plane_position[2] * (ki[2] - kt[2]);
		}
		else {
			T cos_theta_r = cos_theta_i;
			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> 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);

			//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;
			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;

			//calculate the k vectors
			ki = ki_dir * ki_mag;
			kr = kr_dir * kr_mag;
			kt = kt_dir * kt_mag;

			//calculate the Fresnel coefficients
			std::complex<T> rs = std::sin(theta_t - theta_i) / std::sin(theta_t + theta_i);
			std::complex<T> rp = std::tan(theta_t - theta_i) / std::tan(theta_t + theta_i);
			std::complex<T> ts = (2.0 * (sin_theta_t * cos_theta_i)) / std::sin(theta_t + theta_i);
			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;
			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;

			Et[0] = Etp_dir[0] * Et_p + x_hat[0] * Et_s;
			Et[1] = Etp_dir[1] * Et_p + x_hat[1] * Et_s;
			Et[2] = Etp_dir[2] * Et_p + x_hat[2] * Et_s;
		}


		//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));
		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]);

		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()
	{
		std::stringstream ss;
		ss << "k: " << m_k << std::endl;
		ss << "E: " << m_E << std::endl;
		return ss.str();
	}
};					//end planewave class
}					//end namespace optics
}					//end namespace stim

template <typename T>
std::ostream& operator<<(std::ostream& os, stim::optics::planewave<T> p)
{
    os<<p.str();
    return os;
}

#endif