diff options
| -rw-r--r-- | fourier_integrator.cpp | 25 |
1 files changed, 14 insertions, 11 deletions
diff --git a/fourier_integrator.cpp b/fourier_integrator.cpp index e2446fd..7fc8ab4 100644 --- a/fourier_integrator.cpp +++ b/fourier_integrator.cpp @@ -161,26 +161,32 @@ int main(int argc, char* argv[]) { } for (unsigned i = 0; i < Rt.size(); i++) { - Rt[i] += γ * (Rtnew[i] - Rt[i]); - } - - for (unsigned i = 0; i < Rt.size(); i++) { Real ω = i * Δω; Ctnew[i] = (Γ₀ * std::conj(Rt[i]) / (1 + pow(τ₀ * ω, 2)) + pow(y, 2) * (RddfCt[i] * Ct[i] + dfCt[i] * std::conj(Rt[i]))) / (z + 1i * ω); } - for (unsigned i = 0; i < Ct.size(); i++) { - Ct[i] += γ * (Ctnew[i] - Ct[i]); - } - std::vector<Real> Cnew = fft.inverse(Ctnew); std::vector<Real> Rnew = fft.inverse(Rtnew); ΔC = 0; for (unsigned i = 0; i < Cnew.size(); i++) { ΔC += pow(Cnew[i] - C[i], 2); + ΔC += pow(Rnew[i] - R[i], 2); + } + + /* + for (unsigned i = 0; i < Cnew.size(); i++) { + C[i] += γ * (Cnew[i] - C[i]); } + for (unsigned i = 0; i < Rnew.size(); i++) { + R[i] += γ * (Rnew[i] - R[i]); + } + */ + + C = Cnew; + R = Rnew; + z *= Cnew[0]; Real energy = 0; @@ -189,9 +195,6 @@ int main(int argc, char* argv[]) { energy += y * R[i] * df(p, C[i]) * M_PI * Δτ; } - C = Cnew; - R = Rnew; - std::cerr << "y " << y << " iteration " << it << ": " << sqrt(ΔC / C.size()) << " " << pow(y, 2) / z << " " << energy << std::endl; if (it > maxIterations) { |