#ifndef RTS_PLANEWAVE
#define RTS_PLANEWAVE

#include <string>
#include <sstream>

#include "../math/vector.h"
#include "../math/quaternion.h"
#include "../math/constants.h"

/*Basic conversions used here (assuming a vacuum)
	lambda =
*/

namespace rts{

template<typename P>
class planewave{

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)
	{
		k = rts::vec<P>(0, 0, 1) * (TAU/lambda);
		E0 = rts::vec<P>(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))
	{
		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
	}

	///multiplication operator: scale E0
    planewave<P> & operator* (const P & rhs)
	{
		
		E0 = E0 * rhs;
		return *this;
	}

	P lambda() const
	{
		return TAU / k.len();
	}

	P kmag() const
	{
		return k.len();
	}


	planewave<P> refract(rts::vec<P> kn) const
	{
		vec<P> kn_hat = kn.norm();				//normalize the new k
		vec<P> k_hat = k.norm();				//normalize the current k

		vec<P> r = k_hat.cross(kn_hat);			//compute the rotation vector

		P theta = asin(r.len());				//compute the angle of the rotation about r

		planewave<P> new_p;						//create a new plane wave

		//deal with a zero vector (both k and kn point in the same direction)
		if(theta == (P)0)
		{
			new_p = *this;
			return new_p;
		}

		//create a quaternion to capture the rotation
		quaternion<P> q;
		q.CreateRotation(theta, r.norm());

		//apply the rotation to E0
		vec<P> E0n = q.toMatrix3() * E0;

		new_p.k = kn_hat * k.len();
		new_p.E0 = E0n;

		return new_p;

		/*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 )";
		return ss.str();
	}
};
}

#endif