### Blame view

stim/optics/planewave.h 14.4 KB
 ```1 2``` `````` #ifndef TIRA_PLANEWAVE_H #define TIRA_PLANEWAVE_H `````` ```3 4 5``` `````` #include #include `````` `6` `````` #include `````` `7` `````` `````` `8` `````` #include "../math/vec3.h" `````` ```9 10``` `````` #include "../math/quaternion.h" #include "../math/constants.h" `````` `11` `````` #include "../math/plane.h" `````` ```12 13``` `````` #include `````` `14` `````` `````` `15` `````` namespace tira{ `````` ```16 17``` `````` namespace optics{ `````` `18` `````` `````` ```19 20 21 22 23 24 25 26 27 28``` `````` /// evaluate the scalar field produced by a plane wave at a point (x, y, z) /// @param x is the x-coordinate of the point /// @param y is the y-coordinate of the point /// @param z is the z-coordinate of the point /// @param A is the amplitude of the plane wave, specifically the field at (0, 0, 0) /// @param kx is the k-vector component in the x direction /// @param ky is the k-vector component in the y direction /// @param kz is the k-vector component in the z direction template `````` `29` `````` std::complex planewave_scalar(T x, T y, T z, std::complex A, T kx, T ky, T kz){ `````` `30` `````` T d = x * kx + y * ky + z * kz; //calculate the dot product between k and p = (x, y, z) to find the distance p is along the propagation direction `````` `31` `````` std::complex di = std::complex(0, d); //calculate the phase shift that will have to be applied to propagate the wave distance d `````` ```32 33 34 35 36 37 38 39 40 41 42 43 44 45 46``` `````` return A * exp(di); //multiply the phase term by the amplitude at (0, 0, 0) to propagate the wave to p } /// evaluate the scalar field produced by a plane wave at several positions /// @param field is a pre-allocated block of memory that will store the complex field at all points /// @param N is the number of field values to be evaluated /// @param x is a set of x coordinates defining positions within the field (NULL implies that all values are zero) /// @param y is a set of y coordinates defining positions within the field (NULL implies that all values are zero) /// @param z is a set of z coordinates defining positions within the field (NULL implies that all values are zero) /// @param A is the amplitude of the plane wave, specifically the field at (0, 0, 0) /// @param kx is the k-vector component in the x direction /// @param ky is the k-vector component in the y direction /// @param kz is the k-vector component in the z direction template `````` `47` `````` void cpu_planewave_scalar(std::complex* field, size_t N, T* x, T* y = NULL, T* z = NULL, std::complex A = 1.0, T kx = 0.0, T ky = 0.0, T kz = 0.0){ `````` ```48 49 50 51 52 53 54 55 56``` `````` T px, py, pz; for(size_t i = 0; i < N; i++){ // for each element in the array (x == NULL) ? px = 0 : px = x[i]; // test for NULL values (y == NULL) ? py = 0 : py = y[i]; (z == NULL) ? pz = 0 : pz = z[i]; field[i] = planewave_scalar(px, py, pz, A, kx, ky, kz); // call the single-value plane wave function } } `````` `57` `````` `````` `58` `````` `````` `59` `````` template `````` ```60 61``` `````` class planewave{ `````` ```62 63``` `````` protected: `````` ```64 65``` `````` cvec3 m_k; //k-vector, pointed in propagation direction with magnitude |k| = tau / lambda = 2pi / lambda cvec3 m_E; //amplitude (for a scalar plane wave, only E0[0] is used) `````` `66` `````` `````` `67` `````` /// Bend a plane wave via refraction, given that the new propagation direction is known `````` `68` `````` CUDA_CALLABLE planewave bend(vec3 v) const { `````` `69` `````` `````` `70` `````` vec3 k_real(m_k.get(0).real(), m_k.get(1).real(), m_k.get(2).real()); //force the vector to be real (can only refract real directions) `````` `71` `````` `````` ```72 73``` `````` vec3 kn_hat = v.direction(); //normalize the new k vec3 k_hat = k_real.direction(); //normalize the current k `````` `74` `````` `````` `75` `````` planewave new_p; //create a new plane wave `````` `76` `````` `````` `77` `````` T k_dot_kn = k_hat.dot(kn_hat); //if kn is equal to k or -k, handle the degenerate case `````` ```78 79``` `````` //if k . n < 0, then the bend is a reflection `````` `80` `````` if(k_dot_kn < 0) k_hat = -k_hat; //flip k_hat `````` `81` `````` `````` `82` `````` if(k_dot_kn == -1){ `````` ```83 84``` `````` new_p.m_k = -m_k; new_p.m_E = m_E; `````` ```85 86 87``` `````` return new_p; } else if(k_dot_kn == 1){ `````` ```88 89``` `````` new_p.m_k = m_k; new_p.m_E = m_E; `````` ```90 91 92``` `````` return new_p; } `````` `93` `````` vec3 r = k_hat.cross(kn_hat); //compute the rotation vector `````` ```94 95``` `````` T theta = asin(r.len()); //compute the angle of the rotation about r quaternion q; //create a quaternion to capture the rotation `````` `96` `````` q.CreateRotation(theta, r.direction()); `````` `97` `````` matrix_sq R = q.toMatrix3(); `````` ```98 99 100 101 102 103 104 105 106 107 108 109``` `````` vec3< std::complex > E(m_E.get(0), m_E.get(1), m_E.get(2)); vec3< std::complex > E0n = R * E; //apply the rotation to E0 //new_p.m_k = kn_hat * kmag(); //new_p.m_E = E0n; new_p.m_k[0] = kn_hat[0] * kmag(); new_p.m_k[1] = kn_hat[1] * kmag(); new_p.m_k[2] = kn_hat[2] * kmag(); new_p.m_E[0] = E0n[0]; new_p.m_E[1] = E0n[1]; new_p.m_E[2] = E0n[2]; `````` ```110 111 112 113``` `````` return new_p; } `````` `114` `````` public: `````` `115` `````` `````` ```116 117``` `````` `````` `118` `````` ///constructor: create a plane wave propagating along k `````` ```119 120 121 122 123 124 125``` `````` CUDA_CALLABLE planewave(std::complex kx, std::complex ky, std::complex kz, std::complex Ex, std::complex Ey, std::complex Ez) { m_k = cvec3(kx, ky, kz); m_E = cvec3(Ex, Ey, Ez); force_orthogonal(); } `````` `126` `````` `````` `127` `````` CUDA_CALLABLE planewave() : planewave(0, 0, 1, 1, 0, 0) {} `````` `128` `````` `````` ```129 130 131 132 133 134 135 136 137 138``` `````` //copy constructor CUDA_CALLABLE planewave(const planewave& other) { m_k = other.m_k; m_E = other.m_E; } /// Assignment operator CUDA_CALLABLE planewave& operator=(const planewave& rhs) { m_k = rhs.m_k; m_E = rhs.m_E; `````` `139` `````` `````` ```140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177``` `````` return *this; } /// Forces the k and E vectors to be orthogonal CUDA_CALLABLE void force_orthogonal() { /*if (m_E.norm2() == 0) return; cvec3 k_dir = m_k.direction(); //calculate the normalized direction vectors for k and E cvec3 E_dir = m_E.direction(); cvec3 side = k_dir.cross(E_dir); //calculate a side vector for projection cvec3 side_dir = side.direction(); //normalize the side vector E_dir = side_dir.cross(k_dir); //calculate the new E vector direction T E_norm = m_E.norm2(); m_E = E_dir * E_norm; //apply the new direction to the existing E vector */ } CUDA_CALLABLE cvec3 k() { return m_k; } CUDA_CALLABLE cvec3 E() { return m_E; } CUDA_CALLABLE cvec3 evaluate(T x, T y, T z) { std::complex k_dot_r = m_k[0] * x + m_k[1] * y + m_k[2] * z; std::complex e_k_dot_r = std::exp(std::complex(0, 1) * k_dot_r); cvec3 result; result[0] = m_E[0] * e_k_dot_r; result[1] = m_E[1] * e_k_dot_r; result[2] = m_E[2] * e_k_dot_r; return result; } `````` ```178 179 180 181``` `````` CUDA_CALLABLE T kmag() const { return std::sqrt(std::real(m_k.get(0) * std::conj(m_k.get(0)) + m_k.get(1) * std::conj(m_k.get(1)) + m_k.get(2) * std::conj(m_k.get(2)))); } `````` ```182 183 184 185 186 187 188 189 190 191 192 193 194``` `````` /// Return a plane wave with the origin translated by (x, y, z) CUDA_CALLABLE planewave translate(T x, T y, T z) const { planewave result; cvec3 k = m_k; result.m_k = k; std::complex k_dot_r = k[0] * (-x) + k[1] * (-y) + k[2] * (-z); std::complex exp_k_dot_r = std::exp(std::complex(0.0, 1.0) * k_dot_r); cvec3 E = m_E; result.m_E[0] = E[0] * exp_k_dot_r; result.m_E[1] = E[1] * exp_k_dot_r; result.m_E[2] = E[2] * exp_k_dot_r; return result; `````` ```195 196 197``` `````` } ///multiplication operator: scale E0 `````` ```198 199``` `````` CUDA_CALLABLE planewave& operator* (const T& rhs) { m_E = m_E * rhs; `````` ```200 201 202``` `````` return *this; } `````` ```203 204 205 206 207 208 209 210 211 212 213``` `````` ///return a plane wave with the applied refractive index (scales the k-vector by the input) CUDA_CALLABLE planewave ri(T n) { planewave result; result.m_E = m_E; result.m_k = m_k * n; return result; } CUDA_CALLABLE planewave refract(vec3 kn) const { return bend(kn); } `````` `214` `````` /*CUDA_CALLABLE T lambda() const{ `````` `215` `````` return stim::TAU / k.len(); `````` ```216 217``` `````` } `````` `218` `````` CUDA_CALLABLE T kmag() const{ `````` ```219 220 221``` `````` return k.len(); } `````` `222` `````` CUDA_CALLABLE vec< complex > E(){ `````` ```223 224 225 226 227 228 229``` `````` return E0; } CUDA_CALLABLE vec kvec(){ return k; } `````` ```230 231 232``` `````` /// calculate the value of the field produced by the plane wave given a three-dimensional position CUDA_CALLABLE vec< complex > pos(T x, T y, T z){ return pos( stim::vec(x, y, z) ); `````` ```233 234``` `````` } `````` ```235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252``` `````` 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); } `````` `253` `````` `````` `254` `````` `````` `255` `````` `````` ```256 257 258 259 260 261 262 263 264``` `````` /// Calculate the result of a plane wave hitting an interface between two refractive indices /// @param P is a plane representing the position and orientation of the surface /// @param n0 is the refractive index outside of the surface (in the direction of the normal) /// @param n1 is the refractive index inside the surface (in the direction away from the normal) /// @param r is the reflected component of the plane wave /// @param t is the transmitted component of the plane wave void scatter(stim::plane P, T n0, T n1, planewave &r, planewave &t){ scatter(P, n1/n0, r, t); `````` `265` `````` }*/ `````` ```266 267 268 269``` `````` /// Calculate the scattering result when nr = n1/n0 /// @param P is a plane representing the position and orientation of the surface `````` `270` `````` /// @param r is the ratio n1/n0 `````` ```271 272 273``` `````` /// @param n1 is the refractive index inside the surface (in the direction away from the normal) /// @param r is the reflected component of the plane wave /// @param t is the transmitted component of the plane wave `````` `274` `````` `````` `275` `````` int scatter(vec3 plane_normal, vec3 plane_position, std::complex nr, planewave& r, planewave& t) { `````` ```276 277 278 279 280``` `````` 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; } `````` ```281 282 283``` `````` vec3 ki(m_k[0].real(), m_k[1].real(), m_k[2].real()); //force the current k vector to be real vec3 kr; cvec3 kt, Ei, Er, Et; `````` `284` `````` `````` `285` `````` plane_normal = plane_normal.direction(); `````` `286` `````` vec3 k_dir = ki.direction(); //calculate the direction of the incident plane wave `````` `287` `````` `````` ```288 289``` `````` //int facing = plane_face(k_dir, plane_normal); //determine which direction the plane wave is coming in if (k_dir.dot(plane_normal) > 0) { //if the wave hits the back of the plane, invert the plane and nr `````` `290` `````` std::cout << "ERROR: k-vector intersects the wrong side of the boundary." << std::endl; `````` `291` `````` return -1; //the plane wave is impacting the wrong side of the surface `````` ```292 293 294 295 296 297 298 299 300 301 302 303``` `````` } //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 rp = (1.0 - nr) / (1.0 + nr); //compute the Fresnel coefficients std::complex tp = 2.0 / (1.0 + nr); `````` ```304 305 306 307``` `````` 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; `````` ```308 309``` `````` Er = m_E * rp; //compute the E vectors based on the Fresnel coefficients `````` ```310 311 312``` `````` Et = m_E * tp; //calculate the phase offset based on the plane positions `````` ```313 314 315 316 317 318 319 320 321 322 323``` `````` T phase_r = plane_position.dot(ki - kr); std::complex 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; `````` `324` `````` std::complex sin_theta_t = (1.0/nr) * sin(theta_i); //compute the sine of theta_t using Snell's law `````` ```325 326 327 328``` `````` std::complex cos_theta_t = std::sqrt(1.0 - sin_theta_t * sin_theta_t); std::complex theta_t = asin(sin_theta_t); //compute the cosine of theta_t //Define the basis vectors for the calculation (plane of incidence) `````` ```329 330 331 332``` `````` vec3 z_hat = -plane_normal; vec3 plane_perpendicular = plane_normal * k_dir.dot(plane_normal); vec3 y_hat = (k_dir - plane_perpendicular).direction(); vec3 x_hat = y_hat.cross(z_hat); `````` ```333 334``` `````` //calculate the k-vector magnitudes `````` ```335 336 337 338 339``` `````` T ki_mag = ki.norm2(); T kr_mag = ki_mag; std::complex kt_mag = ki_mag * nr; //calculate the k vector directions `````` ```340 341 342``` `````` vec3 ki_dir = y_hat * sin_theta_i + z_hat * cos_theta_i; vec3 kr_dir = y_hat * sin_theta_r - z_hat * cos_theta_r; cvec3 kt_dir; `````` ```343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358``` `````` 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 rs = std::sin(theta_t - theta_i) / std::sin(theta_t + theta_i); std::complex rp = std::tan(theta_t - theta_i) / std::tan(theta_t + theta_i); std::complex ts = (2.0 * (sin_theta_t * cos_theta_i)) / std::sin(theta_t + theta_i); std::complex 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 `````` ```359 360 361``` `````` vec3 Eip_dir = y_hat * cos_theta_i - z_hat * sin_theta_i; vec3 Erp_dir = y_hat * cos_theta_r + z_hat * sin_theta_r; cvec3 Etp_dir; `````` ```362 363 364 365 366``` `````` 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 `````` `367` `````` std::complex Ei_s = m_E.dot(x_hat); `````` ```368 369 370 371 372 373``` `````` std::complex Ei_p = m_E.dot(Eip_dir); std::complex Er_s = rs * Ei_s; std::complex Er_p = rp * Ei_p; std::complex Et_s = ts * Ei_s; std::complex Et_p = tp * Ei_p; `````` `374` `````` //calculate the E vector for each plane wave `````` ```375 376 377 378 379 380 381 382 383 384 385``` `````` 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 `````` ```386 387 388 389 390 391``` `````` T phase_r = plane_position.dot(ki - kr); std::complex phase_t = plane_position[0] * (ki[0] - kt[0]) + plane_position[1] * (ki[1] - kt[1]) + plane_position[2] * (ki[2] - kt[2]); `````` `392` `````` //apply the phase offset `````` ```393 394 395``` `````` Er = Er * std::exp(std::complex(0, 1) * phase_r); Et = Et * std::exp(std::complex(0, 1) * phase_t); `````` ```396 397 398``` `````` //generate the reflected and transmitted waves 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]); `````` ```399 400``` `````` return 0; `````` `401` `````` } `````` ```402 403 404 405``` `````` std::string str() { std::stringstream ss; `````` ```406 407``` `````` ss << "k: " << m_k << std::endl; ss << "E: " << m_E << std::endl; `````` ```408 409``` `````` return ss.str(); } `````` ```410 411``` `````` }; //end planewave class } //end namespace optics `````` `412` `````` } //end namespace tira `````` `413` `````` `````` `414` `````` template `````` `415` `````` std::ostream& operator<<(std::ostream& os, tira::optics::planewave p) `````` ```416 417 418 419 420``` `````` { os<