integrate boundary scattering into the plane wave class

David Mayerich
1 parent 14634fca
Showing 1 changed file with 148 additions and 7 deletions Show diff stats
stim/optics/planewave.h
@@ -345,12 +345,12 @@ public:
 	}
 	///multiplication operator: scale E0
-    /*CUDA_CALLABLE planewave<T>& operator* (const T& rhs) {
-		E0 = E0 * rhs;
+    CUDA_CALLABLE planewave<T>& operator* (const T& rhs) {
+		m_E = m_E * rhs;
 		return *this;
 	}
-	CUDA_CALLABLE T lambda() const{
+	/*CUDA_CALLABLE T lambda() const{
 		return stim::TAU / k.len();
 	}
@@ -413,9 +413,150 @@ public:
 	/// @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) {
+	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;
-		//TODO: Generate a real vector from the current K vector - this will not address complex k-vectors
+		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];
@@ -525,8 +666,8 @@ public:
 		//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]);
-	}*/
+		t = planewave(kt[0], kt[1], kt[2], Et[0], Et[1], Et[2]);*/
+	}
 	std::string str()
 	{