| 264 | } |
| 265 | |
| 266 | void dft_chunk::update_dft(double time) { |
| 267 | if (!fc->f[c][0]) return; |
| 268 | |
| 269 | const int Nomega = omega.size(); |
| 270 | for (int i = 0; i < Nomega; ++i) |
| 271 | dft_phase[i] = polar(1.0, omega[i] * time) * scale; |
| 272 | |
| 273 | int numcmp = fc->f[c][1] ? 2 : 1; |
| 274 | |
| 275 | PLOOP_OVER_IVECS(fc->gv, is, ie, idx) { |
| 276 | size_t idx_dft = IVEC_LOOP_COUNTER; |
| 277 | double w; |
| 278 | if (include_dV_and_interp_weights) { |
| 279 | w = IVEC_LOOP_WEIGHT(s0, s1, e0, e1, dV0 + dV1 * loop_i2); |
| 280 | if (sqrt_dV_and_interp_weights) w = sqrt(w); |
| 281 | } |
| 282 | else |
| 283 | w = 1.0; |
| 284 | realnum f[2]; // real/imag field value at epsilon point |
| 285 | if (avg2) |
| 286 | for (int cmp = 0; cmp < numcmp; ++cmp) |
| 287 | f[cmp] = (w * 0.25) * (fc->f[c][cmp][idx] + fc->f[c][cmp][idx + avg1] + |
| 288 | fc->f[c][cmp][idx + avg2] + fc->f[c][cmp][idx + (avg1 + avg2)]); |
| 289 | else if (avg1) |
| 290 | for (int cmp = 0; cmp < numcmp; ++cmp) |
| 291 | f[cmp] = (w * 0.5) * (fc->f[c][cmp][idx] + fc->f[c][cmp][idx + avg1]); |
| 292 | else |
| 293 | for (int cmp = 0; cmp < numcmp; ++cmp) |
| 294 | f[cmp] = w * fc->f[c][cmp][idx]; |
| 295 | |
| 296 | if (numcmp == 2) { |
| 297 | complex<realnum> fc(f[0], f[1]); |
| 298 | for (int i = 0; i < Nomega; ++i) |
| 299 | dft[Nomega * idx_dft + i] += dft_phase[i] * fc; |
| 300 | } |
| 301 | else { |
| 302 | realnum fr = f[0]; |
| 303 | for (int i = 0; i < Nomega; ++i) |
| 304 | dft[Nomega * idx_dft + i] += |
| 305 | std::complex<realnum>{fr * dft_phase[i].real(), fr * dft_phase[i].imag()}; |
| 306 | } |
| 307 | } |
| 308 | } |
| 309 | |
| 310 | /* Return the L2 norm of the DFTs themselves. This is useful |
| 311 | to check whether the simulation is finished (whether all relevant fields have decayed). |
no test coverage detected