#ifndef RTS_PLANEWAVE #define RTS_PLANEWAVE #include #include #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{ namespace optics{ template class planewave{ protected: vec k; //k = tau / lambda vec< complex > E0; //amplitude //T phi; CUDA_CALLABLE planewave bend(rts::vec kn) const{ vec kn_hat = kn.norm(); //normalize the new k vec k_hat = k.norm(); //normalize the current k //std::cout<<"PLANE WAVE BENDING------------------"< 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: "< r = k_hat.cross(kn_hat); //compute the rotation vector //std::cout<<"r: "< q; q.CreateRotation(theta, r.norm()); //apply the rotation to E0 vec< complex > 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(0, 0, 1) * (TAU/lambda); E0 = rts::vec(1, 0, 0); }*/ ///constructor: create a plane wave propagating along k, polarized along _E0, at frequency _omega CUDA_CALLABLE planewave(vec kvec = rts::vec(0, 0, rtsTAU), vec< complex > E = rts::vec(1, 0, 0), T phase = 0) { //phi = phase; k = kvec; vec< complex > k_hat = k.norm(); if(E.len() == 0) //if the plane wave has an amplitude of 0 E0 = vec(0); //just return it else{ vec< complex > s = (k_hat.cross(E)).norm(); //compute an orthogonal side vector vec< complex > 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(0, phase) ); } ///multiplication operator: scale E0 CUDA_CALLABLE planewave & 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 > E(){ return E0; } CUDA_CALLABLE vec kvec(){ return k; } /*CUDA_CALLABLE T phase(){ return phi; } CUDA_CALLABLE void phase(T p){ phi = p; }*/ CUDA_CALLABLE vec< complex > pos(vec p = vec(0, 0, 0)){ vec< complex > result; T kdp = k.dot(p); complex x = complex(0, kdp); complex 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 n(T ni, T nt){ return planewave(k * (nt / ni), E0); } CUDA_CALLABLE planewave refract(rts::vec kn) const { return bend(kn); } void scatter(rts::plane P, T nr, planewave &r, planewave &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 kr = -k; vec kt = k * nr; //set the k vectors for theta_i = 0 vec< complex > Er = E0 * rp; //compute the E vectors vec< complex > 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"<(kr, Er, phase_r); t = planewave(kt, Et, phase_t); //std::cout<<"i + r: "< z_hat = -P.norm(); vec y_hat = P.parallel(k).norm(); vec x_hat = y_hat.cross(z_hat).norm(); //compute the k vectors for r and t vec 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 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 > cx_hat = x_hat; complex 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 Er_s = Ei_s * rs; complex Er_p = Ei_p * rp; //compute the magnitude of the p- and s-polarized components of the transmitted E vector complex Et_s = Ei_s * ts; complex Et_p = Ei_p * tp; //std::cout<<"E0: "< > Er = vec< complex >(y_hat * cos(theta_i) + z_hat * sin(theta_i)) * Er_p + cx_hat * Er_s; //compute the transmitted E vector vec< complex > Et = vec< complex >(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: "<(0, phase_r) ); //r.phi = phase_r; //t = bend(kt); //t.k = t.k * nr; t.k = kt; t.E0 = Et * exp( complex(0, phase_t) ); //t.phi = phase_t; //std::cout<<"i: "< std::ostream& operator<<(std::ostream& os, rts::planewave p) { os<