Tweaking form of iterations

2 files changed, 11 insertions, 2 deletions

diff --git a/log-fourier.cpp b/log-fourier.cpp
index 1e2e7d9..6149e6f 100644
--- a/log-fourier.cpp
+++ b/log-fourier.cpp
@@ -52,6 +52,7 @@ Complex gamma(Complex z) {
 std::vector<Complex> LogarithmicFourierTransform::fourier(const std::vector<Real>& c, bool symmetric) {
   std::vector<Complex> ĉ(N);
   std::vector<Real> σs = {1};
+  /* c is either even or zero for negative arguments */
   if (symmetric){
     σs.push_back(-1);
   }
diff --git a/log-fourier_integrator.cpp b/log-fourier_integrator.cpp
index 061cb08..fd8b220 100644
--- a/log-fourier_integrator.cpp
+++ b/log-fourier_integrator.cpp
@@ -100,12 +100,20 @@ int main(int argc, char* argv[]) {
 
     for (unsigned n = 0; n < N; n++) {
       Rtnew[n] = (1.0 + pow(β, 2) * RddfCt[n] * Rt[n]) / (μ + 1i * fft.ν(n));
-      Ctnew[n] = (2 * Γ₀ * std::conj(Rtnew[n]) / (1 + pow(τ₀ * fft.ν(n), 2)) + pow(β, 2) * (RddfCt[n] * Ct[n] + dfCt[n] * std::conj(Rtnew[n]))) / (μ + 1i * fft.ν(n));
 //      Ctnew[n] = - 2 * Γ₀ * Rtnew[n].imag() / (1 + pow(τ₀ * fft.ν(n), 2)) / fft.ν(n);
     }
+    std::vector<Real> Rnew = fft.inverse(Rtnew);
+    for (unsigned n = 0; n < N; n++) {
+      RddfC[n] = Rnew[n] * ddf(λ, p, s, C[n]);
+    }
+    RddfCt = fft.fourier(RddfC, false);
+
+    for (unsigned n = 0; n < N; n++) {
+      Ctnew[n] = (2 * Γ₀ * std::conj(Rtnew[n]) / (1 + pow(τ₀ * fft.ν(n), 2)) + pow(β, 2) * (RddfCt[n] * Ct[n] + dfCt[n] * std::conj(Rtnew[n]))) / (μ + 1i * fft.ν(n));
+    }
+
 
     std::vector<Real> Cnew = fft.inverse(Ctnew);
-    std::vector<Real> Rnew = fft.inverse(Rtnew);
 
     ΔC = 0;
     for (unsigned i = 0; i < N; i++) {