Changed interative scheme and many variable names

1 files changed, 75 insertions, 73 deletions

diff --git a/log-fourier_integrator.cpp b/log-fourier_integrator.cpp
index fd8b220..bf83da5 100644
--- a/log-fourier_integrator.cpp
+++ b/log-fourier_integrator.cpp
@@ -64,98 +64,100 @@ int main(int argc, char* argv[]) {
     μ = (sqrt(1+4*Γ₀*τ₀)-1)/(2*τ₀);
   }
 
-  std::vector<Real> C(N);
-  std::vector<Real> R(N);
-  std::vector<Complex> Ct(N);
-  std::vector<Complex> Rt(N);
+  std::vector<Real> Cₜ(N);
+  std::vector<Real> Rₜ(N);
+  std::vector<Complex> Ĉₜ(N);
+  std::vector<Complex> Ȓₜ(N);
 
   // start from the exact solution for β = 0
   for (unsigned n = 0; n < N; n++) {
     if (τ₀ > 0) {
-      C[n] = Γ₀ * (exp(-μ * fft.t(n)) - μ * τ₀ * exp(-fft.t(n) / τ₀)) / (μ - pow(μ, 3) * pow(τ₀, 2));
+      Cₜ[n] = Γ₀ * (exp(-μ * fft.t(n)) - μ * τ₀ * exp(-fft.t(n) / τ₀)) / (μ - pow(μ, 3) * pow(τ₀, 2));
     } else {
-      C[n] = Γ₀ * exp(-μ * fft.t(n)) / μ;
+      Cₜ[n] = Γ₀ * exp(-μ * fft.t(n)) / μ;
     }
-    R[n] = exp(-μ * fft.t(n));
+    Rₜ[n] = exp(-μ * fft.t(n));
 
-    Ct[n] = 2 * Γ₀ / (pow(μ, 2) + pow(fft.ν(n), 2)) / (1 + pow(τ₀ * fft.ν(n), 2));
-    Rt[n] = 1.0 / (μ + 1i * fft.ν(n));
+    Ĉₜ[n] = 2 * Γ₀ / (pow(μ, 2) + pow(fft.ν(n), 2)) / (1 + pow(τ₀ * fft.ν(n), 2));
+    Ȓₜ[n] = 1.0 / (μ + 1i * fft.ν(n));
   }
 
   Real β = 0;
   while (β < βₘₐₓ) {
-  Real ΔC = 100;
-  while (ΔC > ε) {
-    std::vector<Real> RddfC(N);
-    std::vector<Real> dfC(N);
-    for (unsigned n = 0; n < N; n++) {
-      RddfC[n] = R[n] * ddf(λ, p, s, C[n]);
-      dfC[n] = df(λ, p, s, C[n]);
-    }
-    std::vector<Complex> RddfCt = fft.fourier(RddfC, false);
-    std::vector<Complex> dfCt = fft.fourier(dfC, true);
-
-    std::vector<Complex> Rtnew(N);
-    std::vector<Complex> Ctnew(N);
-
-    for (unsigned n = 0; n < N; n++) {
-      Rtnew[n] = (1.0 + pow(β, 2) * RddfCt[n] * Rt[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));
-    }
-
+    Real ΔC = 100;
+    while (ΔC > ε) {
+      std::vector<Real> RddfC(N);
+//      std::vector<Real> dfC(N);
+      for (unsigned n = 0; n < N; n++) {
+        RddfC[n] = Rₜ[n] * ddf(λ, p, s, Cₜ[n]);
+//        dfC[n] = df(λ, p, s, Cₜ[n]);
+      }
+      std::vector<Complex> RddfCt = fft.fourier(RddfC, false);
+//      std::vector<Complex> dfCt = fft.fourier(dfC, true);
+
+      std::vector<Complex> Ȓₜ₊₁(N);
+      std::vector<Complex> Ĉₜ₊₁(N);
+
+      for (unsigned n = 0; n < N; n++) {
+        Ȓₜ₊₁[n] = (1.0 + pow(β, 2) * RddfCt[n] * Ȓₜ[n]) / (μ + 1i * fft.ν(n));
+        Ĉₜ₊₁[n] = - 2 * Γ₀ * Ȓₜ₊₁[n].imag() / (1 + pow(τ₀ * fft.ν(n), 2)) / fft.ν(n);
+      }
+
+      std::vector<Real> Rₜ₊₁ = fft.inverse(Ȓₜ₊₁);
+      /*
+      for (unsigned n = 0; n < N; n++) {
+        RddfC[n] = Rₜ₊₁[n] * ddf(λ, p, s, Cₜ[n]);
+      }
+      RddfCt = fft.fourier(RddfC, false);
+
+      for (unsigned n = 0; n < N; n++) {
+        Ĉₜ₊₁[n] = (2 * Γ₀ * std::conj(Ȓₜ₊₁[n]) / (1 + pow(τ₀ * fft.ν(n), 2)) + pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ₊₁[n]))) / (μ + 1i * fft.ν(n));
+      }
+      */
+
+      std::vector<Real> Cₜ₊₁ = fft.inverse(Ĉₜ₊₁);
+
+      ΔC = 0;
+      for (unsigned i = 0; i < N; i++) {
+        ΔC += std::norm(Ĉₜ[i] - Ĉₜ₊₁[i]);
+        ΔC += std::norm(Ȓₜ[i] - Ȓₜ₊₁[i]);
+      }
+      ΔC = sqrt(ΔC) / (2*N);
+
+      for (unsigned i = 0; i < N; i++) {
+        Cₜ[i] += γ * (Cₜ₊₁[i] - Cₜ[i]);
+        Rₜ[i] += γ * (Rₜ₊₁[i] - Rₜ[i]);
+        Ĉₜ[i] += γ * (Ĉₜ₊₁[i] - Ĉₜ[i]);
+        Ȓₜ[i] += γ * (Ȓₜ₊₁[i] - Ȓₜ[i]);
+      }
+
+      μ *= Cₜ[0];
 
-    std::vector<Real> Cnew = fft.inverse(Ctnew);
-
-    ΔC = 0;
-    for (unsigned i = 0; i < N; i++) {
-      ΔC += std::norm(Ct[i] - Ctnew[i]);
-      ΔC += std::norm(Rt[i] - Rtnew[i]);
+//    std::cerr << ΔC << std::endl;
     }
-    ΔC = sqrt(ΔC) / (2*N);
 
-    for (unsigned i = 0; i < N; i++) {
-      C[i] += γ * (Cnew[i] - C[i]);
-      R[i] += γ * (Rnew[i] - R[i]);
-      Ct[i] += γ * (Ctnew[i] - Ct[i]);
-      Rt[i] += γ * (Rtnew[i] - Rt[i]);
+    /* Integrate the energy using Simpson's rule */
+    Real E = 0;
+    for (unsigned n = 0; n < N/2-1; n++) {
+      Real h₂ₙ   = fft.t(2*n+1) - fft.t(2*n);
+      Real h₂ₙ₊₁ = fft.t(2*n+2) - fft.t(2*n+1);
+      Real f₂ₙ   = Rₜ[2*n]   * df(λ, p, s, Cₜ[2*n]);
+      Real f₂ₙ₊₁ = Rₜ[2*n+1] * df(λ, p, s, Cₜ[2*n+1]);
+      Real f₂ₙ₊₂ = Rₜ[2*n+2] * df(λ, p, s, Cₜ[2*n+2]);
+      E += (h₂ₙ + h₂ₙ₊₁) / 6 * (
+            (2 - h₂ₙ₊₁ / h₂ₙ) * f₂ₙ
+            + pow(h₂ₙ + h₂ₙ₊₁, 2) / (h₂ₙ * h₂ₙ₊₁) * f₂ₙ₊₁
+            + (2 - h₂ₙ / h₂ₙ₊₁) * f₂ₙ₊₂
+          );
     }
+    E *= β;
 
-    μ *= C[0];
-
-//    std::cerr << ΔC << std::endl;
-  }
-
-  /* Integrate the energy using Simpson's rule */
-  Real E = 0;
-  for (unsigned i = 0; i < N/2-1; i++) {
-    Real h₂ᵢ   = fft.t(2*i+1) - fft.t(2*i);
-    Real h₂ᵢ₊₁ = fft.t(2*i+2) - fft.t(2*i+1);
-    Real f₂ᵢ   = R[2*i]   * df(λ, p, s, C[2*i]);
-    Real f₂ᵢ₊₁ = R[2*i+1] * df(λ, p, s, C[2*i+1]);
-    Real f₂ᵢ₊₂ = R[2*i+2] * df(λ, p, s, C[2*i+2]);
-    E += (h₂ᵢ + h₂ᵢ₊₁) / 6 * (
-          (2 - h₂ᵢ₊₁ / h₂ᵢ) * f₂ᵢ
-          + pow(h₂ᵢ + h₂ᵢ₊₁, 2) / (h₂ᵢ * h₂ᵢ₊₁) * f₂ᵢ₊₁
-          + (2 - h₂ᵢ / h₂ᵢ₊₁) * f₂ᵢ₊₂
-        );
-  }
-  E *= β;
-
-  std::cerr << β << " " << μ << " " << Ct[0].real() << " " << E << std::endl;
-   β += Δβ;
+    std::cerr << β << " " << μ << " " << Ĉₜ[0].real() << " " << E << std::endl;
+    β += Δβ;
   }
 
   for (unsigned i = 0; i < N; i++) {
-    std::cout << fft.t(i) << " " << C[i] << std::endl;
+    std::cout << fft.t(i) << " " << Cₜ[i] << std::endl;
   }
 
   return 0;