Set Vinv = inverse of V, where both V and Vinv are complex matrices.*/
| 244 | |
| 245 | /* Set Vinv = inverse of V, where both V and Vinv are complex matrices.*/ |
| 246 | void matrix_invert(std::complex<double> (&Vinv)[9], std::complex<double> (&V)[9]) { |
| 247 | |
| 248 | std::complex<double> det = |
| 249 | (V[0 + 3 * 0] * (V[1 + 3 * 1] * V[2 + 3 * 2] - V[1 + 3 * 2] * V[2 + 3 * 1]) - |
| 250 | V[0 + 3 * 1] * (V[0 + 3 * 1] * V[2 + 3 * 2] - V[1 + 3 * 2] * V[0 + 3 * 2]) + |
| 251 | V[0 + 3 * 2] * (V[0 + 3 * 1] * V[1 + 3 * 2] - V[1 + 3 * 1] * V[0 + 3 * 2])); |
| 252 | |
| 253 | if (det == 0.0) meep::abort("meep: Matrix is singular, aborting.\n"); |
| 254 | |
| 255 | Vinv[0 + 3 * 0] = 1.0 / det * (V[1 + 3 * 1] * V[2 + 3 * 2] - V[1 + 3 * 2] * V[2 + 3 * 1]); |
| 256 | Vinv[0 + 3 * 1] = 1.0 / det * (V[0 + 3 * 2] * V[2 + 3 * 1] - V[0 + 3 * 1] * V[2 + 3 * 2]); |
| 257 | Vinv[0 + 3 * 2] = 1.0 / det * (V[0 + 3 * 1] * V[1 + 3 * 2] - V[0 + 3 * 2] * V[1 + 3 * 1]); |
| 258 | Vinv[1 + 3 * 0] = 1.0 / det * (V[1 + 3 * 2] * V[2 + 3 * 0] - V[1 + 3 * 0] * V[2 + 3 * 2]); |
| 259 | Vinv[1 + 3 * 1] = 1.0 / det * (V[0 + 3 * 0] * V[2 + 3 * 2] - V[0 + 3 * 2] * V[2 + 3 * 0]); |
| 260 | Vinv[1 + 3 * 2] = 1.0 / det * (V[0 + 3 * 2] * V[1 + 3 * 0] - V[0 + 3 * 0] * V[1 + 3 * 2]); |
| 261 | Vinv[2 + 3 * 0] = 1.0 / det * (V[1 + 3 * 0] * V[2 + 3 * 1] - V[1 + 3 * 1] * V[2 + 3 * 0]); |
| 262 | Vinv[2 + 3 * 1] = 1.0 / det * (V[0 + 3 * 1] * V[2 + 3 * 0] - V[0 + 3 * 0] * V[2 + 3 * 1]); |
| 263 | Vinv[2 + 3 * 2] = 1.0 / det * (V[0 + 3 * 0] * V[1 + 3 * 1] - V[0 + 3 * 1] * V[1 + 3 * 0]); |
| 264 | } |
| 265 | |
| 266 | complex<double> structure_chunk::get_chi1inv_at_pt(component c, direction d, int idx, |
| 267 | double frequency, |
no test coverage detected