codebase / stimlib

  #ifndef RTS_PLANEWAVE
  #define RTS_PLANEWAVE
  
  #include <string>
  #include <sstream>
  

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

  #include "../math/plane.h"

  #include "../cuda/callable.h"

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

  namespace stim{

  template<typename T>

  class planewave{
  

  protected:
  
  	vec<T> k;	//k = tau / lambda

  	vec< complex<T> > E0;		//amplitude
  	//T phi;

  	CUDA_CALLABLE planewave<T> bend(rts::vec<T> kn) const{

  
  		vec<T> kn_hat = kn.norm();				//normalize the new k
  		vec<T> k_hat = k.norm();				//normalize the current k
  
  		//std::cout<<"PLANE WAVE BENDING------------------"<<std::endl;
  		//std::cout<<"kn_hat: "<<kn_hat<<"     k_hat: "<<k_hat<<std::endl;
  
  		planewave<T> new_p;						//create a new plane wave
  
  		//if kn is equal to k or -k, handle the degenerate case
  		T k_dot_kn = k_hat.dot(kn_hat);
  
  		//if k . n < 0, then the bend is a reflection
  			//flip k_hat
  		if(k_dot_kn < 0) k_hat = -k_hat;
  
  		//std::cout<<"k dot kn: "<<k_dot_kn<<std::endl;
  
  		//std::cout<<"k_dot_kn: "<<k_dot_kn<<std::endl;
  		if(k_dot_kn == -1){
  			new_p.k = -k;
  			new_p.E0 = E0;
  			return new_p;
  		}
  		else if(k_dot_kn == 1){
  			new_p.k = k;
  			new_p.E0 = E0;
  			return new_p;
  		}
  
  		vec<T> r = k_hat.cross(kn_hat);			//compute the rotation vector
  
  		//std::cout<<"r: "<<r<<std::endl;
  
  		T theta = asin(r.len());				//compute the angle of the rotation about r
  
  		
  
  		//deal with a zero vector (both k and kn point in the same direction)
  		//if(theta == (T)0)
  		//{
  		//	new_p = *this;
  		//	return new_p;
  		//}
  
  		//create a quaternion to capture the rotation
  		quaternion<T> q;
  		q.CreateRotation(theta, r.norm());
  
  		//apply the rotation to E0

  		vec< complex<T> > E0n = q.toMatrix3() * E0;

  
  		new_p.k = kn_hat * kmag();
  		new_p.E0 = E0n;
  
  		return new_p;
  	}
  

  public:

  	///constructor: create a plane wave propagating along z, polarized along x

  	/*planewave(T lambda = (T)1)

  	{

  		k = rts::vec<T>(0, 0, 1) * (TAU/lambda);
  		E0 = rts::vec<T>(1, 0, 0);

  	}*/

  	///constructor: create a plane wave propagating along k, polarized along _E0, at frequency _omega
  	CUDA_CALLABLE planewave(vec<T> kvec = rts::vec<T>(0, 0, rtsTAU), 

  							vec< complex<T> > E = rts::vec<T>(1, 0, 0), T phase = 0)

  	{

  		//phi = phase;
  

  		k = kvec;

  		vec< complex<T> > k_hat = k.norm();

  
  		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
  		}

  
  		E0 = E0 * exp( complex<T>(0, phase) );

  	}
  
  	///multiplication operator: scale E0

      CUDA_CALLABLE planewave<T> & operator* (const T & rhs)

  	{
  		
  		E0 = E0 * rhs;
  		return *this;
  	}
  

  	CUDA_CALLABLE T lambda() const

  	{

  		return rtsTAU / 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;
  	}
  

  	/*CUDA_CALLABLE T phase(){

  		return phi;
  	}
  
  	CUDA_CALLABLE void phase(T p){
  		phi = p;

  	}*/

  
  	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(rts::vec<T> kn) const

  	{

  		return bend(kn);
  	}

  	void scatter(rts::plane<T> P, T nr, planewave<T> &r, planewave<T> &t){

  		int facing = P.face(k);		//determine which direction the plane wave is coming in

  		//if(facing == 0)				//if the wave is tangent to the plane, return an identical wave
  		//	return *this;
  		//else 
  		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)
  		}

  		//use Snell's Law to calculate the transmitted angle
  		T cos_theta_i = k.norm().dot(-P.norm());				//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
  
  		bool tir = false;						//flag for total internal reflection
  		if(theta_t != theta_t){
  			tir = true;
  			theta_t = rtsPI / (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::cout<<"Degeneracy: Head-On"<<std::endl;
  			//std::cout<<"rs: "<<rp<<"  rp: "<<rp<<"  ts: "<<tp<<"  tp: "<<tp<<std::endl;
  			//std::cout<<"phase r: "<<phase_r<<"  phase t: "<<phase_t<<std::endl;
  
  			//create the plane waves
  			r = planewave<T>(kr, Er, phase_r);
  			t = planewave<T>(kt, Et, phase_t);
  
  			//std::cout<<"i + r: "<<pos()[0] + r.pos()[0]<<pos()[1] + r.pos()[1]<<pos()[2] + r.pos()[2]<<std::endl;
  			//std::cout<<"t:     "<<t.pos()[0]<<t.pos()[1]<<t.pos()[2]<<std::endl;
  			//std::cout<<"--------------------------------"<<std::endl;
  			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
  		vec<T> z_hat = -P.norm();
  		vec<T> y_hat = P.parallel(k).norm();
  		vec<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;
  
  		//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 = (0 < E0.dot(y_hat)) - (E0.dot(y_hat) < 0);
  		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;

  		//T Ei_p = ( E0 - x_hat * Ei_s ).len();
  		//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;

  		//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::cout<<"E0: "<<E0<<std::endl;
  		//std::cout<<"E0 dot y_hat: "<<E0.dot(y_hat)<<std::endl;
  		//std::cout<<"theta i: "<<theta_i<<"  theta t: "<<theta_t<<std::endl;
  		//std::cout<<"x_hat: "<<x_hat<<"  y_hat: "<<y_hat<<"  z_hat: "<<z_hat<<std::endl;
  		//std::cout<<"Ei_s: "<<Ei_s<<"  Ei_p: "<<Ei_p<<"  Er_s: "<<Er_s<<"  Er_p: "<<Er_p<<"  Et_s: "<<Et_s<<"  Et_p: "<<Et_p<<std::endl;
  		//std::cout<<"rs: "<<rs<<"  rp: "<<rp<<"  ts: "<<ts<<"  tp: "<<tp<<std::endl;
  		
  
  		//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;

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

  
  		T phase_t = P.p().dot(k - kt);
  		T phase_r = P.p().dot(k - kr);
  
  		//std::cout<<"phase r: "<<phase_r<<"  phase t: "<<phase_t<<std::endl;
  
  		//std::cout<<"phase: "<<phase<<std::endl;
  
  		//create the plane waves
  		r.k = kr;

  		r.E0 = Er * exp( complex<T>(0, phase_r) );
  		//r.phi = phase_r;

  
  		//t = bend(kt);
  		//t.k = t.k * nr;
  
  		t.k = kt;

  		t.E0 = Et * exp( complex<T>(0, phase_t) );
  		//t.phi = phase_t;

  		//std::cout<<"i: "<<str()<<std::endl;
  		//std::cout<<"r: "<<r.str()<<std::endl;
  		//std::cout<<"t: "<<t.str()<<std::endl;
  
  		//std::cout<<"i + r: "<<pos()[0] + r.pos()[0]<<pos()[1] + r.pos()[1]<<pos()[2] + r.pos()[2]<<std::endl;
  		//std::cout<<"t:     "<<t.pos()[0]<<t.pos()[1]<<t.pos()[2]<<std::endl;
  		//std::cout<<"--------------------------------"<<std::endl;

  	}
  
  	std::string str()
  	{
  		std::stringstream ss;

  		ss<<"Plane Wave:"<<std::endl;
  		ss<<"	"<<E0<<" e^i ( "<<k<<" . r )";

  		return ss.str();
  	}
  };
  }
  

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

  #endif