codebase / stimlib

  #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