1 files changed, 33 insertions, 21 deletions

diff --git a/fourier_integrator.cpp b/fourier_integrator.cpp
index 86c3b24..4f4eff5 100644
--- a/fourier_integrator.cpp
+++ b/fourier_integrator.cpp
@@ -17,7 +17,7 @@ int main(int argc, char* argv[]) {
 
   unsigned maxIterations = 1000;
   Real ε = 1e-14;
-  Real γ = 0;
+  Real γ = 1;
 
   bool loadData = false;
 
@@ -65,11 +65,11 @@ int main(int argc, char* argv[]) {
 
   unsigned n = pow(2, log2n);
 
-  Real Δτ = τₘₐₓ / M_PI / n;
-  Real Δω = M_PI / τₘₐₓ;
+  Real Δτ = (1 + τ₀ / 2) * τₘₐₓ / M_PI / n;
+  Real Δω = M_PI / ((1 + τ₀ / 2) * τₘₐₓ);
 
-  Real z = (-1+sqrt(1+2*τ₀)) / (2 * τ₀);
-  Real Γ₀ = 1;
+  Real z = 0.5;
+  Real Γ₀ = 1 + τ₀ / 2;
 
   std::vector<Real> C(2 * n);
   std::vector<Real> R(2 * n);
@@ -84,7 +84,11 @@ int main(int argc, char* argv[]) {
     // start from the exact solution for τ₀ = 0
     for (unsigned i = 0; i < n; i++) {
       Real τ = i * Δτ * M_PI;
-      C[i] = Γ₀ / 2 * (exp(-z * τ) - z * τ₀ * exp(-τ / τ₀)) / (z - pow(z, 3) * pow(τ₀, 2));
+      if (τ₀ > 0) {
+        C[i] = Γ₀ / 2 * (exp(-z * τ) - z * τ₀ * exp(-τ / τ₀)) / (z - pow(z, 3) * pow(τ₀, 2));
+      } else {
+        C[i] = Γ₀ / 2 * exp(-z * τ) / z;
+      }
       if (i > 0) {
         C[2 * n - i] = C[i];
       }
@@ -94,24 +98,26 @@ int main(int argc, char* argv[]) {
     Rt = fft.fourier(R);
   } else {
     std::ifstream cfile(fourierFile("C", p, s, λ, τ₀, y, log2n, τₘₐₓ), std::ios::binary);
-    cfile.read((char*)(C.data()), C.size() * sizeof(Real));
+    cfile.read((char*)(C.data()), (C.size() / 2) * sizeof(Real));
     cfile.close();
+    for (unsigned i = 1; i < n; i++) {
+      C[2 * n - i] = C[i];
+    }
     std::ifstream rfile(fourierFile("R", p, s, λ, τ₀, y, log2n, τₘₐₓ), std::ios::binary);
-    rfile.read((char*)(R.data()), R.size() * sizeof(Real));
+    rfile.read((char*)(R.data()), (R.size() / 2) * sizeof(Real));
     rfile.close();
 
     Ct = fft.fourier(C);
     Rt = fft.fourier(R);
 
-    auto [RddfCt, dfCt] = RddfCtdfCt(fft, C, R, p, s, λ);
-
-    z = ((Γ₀ * std::conj(Rt[0]) + pow(y, 2) * (RddfCt[0] * Ct[0] + dfCt[0] * std::conj(Rt[0]))) / Ct[0]).real();
+    z = estimateZ(fft, C, Ct, R, Rt, p, s, λ, τ₀, y);
   }
 
   while (y += Δy, y <= yₘₐₓ) {
-    Real ΔC = 1;;
+    Real ΔC = 1;
+    Real ΔCprev = 1000;
     unsigned it = 0;
-    while (sqrt(ΔC / C.size()) > ε) {
+    while (sqrt(2 * ΔC / C.size()) > ε) {
       it++;
       auto [RddfCt, dfCt] = RddfCtdfCt(fft, C, R, p, s, λ);
 
@@ -132,7 +138,7 @@ int main(int argc, char* argv[]) {
       }
 
       ΔC = 0;
-      for (unsigned i = 0; i < Cnew.size(); i++) {
+      for (unsigned i = 0; i < Cnew.size() / 2; i++) {
         ΔC += pow(Cnew[i] - C[i], 2);
         ΔC += pow(Rnew[i] - R[i], 2);
       }
@@ -141,18 +147,24 @@ int main(int argc, char* argv[]) {
         C[i] += γ * (Cnew[i] - C[i]);
       }
 
-      for (unsigned i = 0; i < Rnew.size(); i++) {
+      for (unsigned i = 0; i < Rnew.size() / 2; i++) {
         R[i] += γ * (Rnew[i] - R[i]);
       }
 
       z *= Cnew[0];
 
-      std::cerr << it << " " << p << " " << s << " " << τ₀ << " " << sqrt(ΔC / C.size()) << std::endl;
+      if (it % maxIterations == 0) {
+        if (ΔC < ΔCprev) {
+          γ = std::min(1.1 * γ, 1.0);
+        } else {
+          γ /= 2;
+        }
 
-      if (it > maxIterations) {
-        it = 0;
-        γ /= 2;
+        ΔCprev = ΔC;
       }
+
+      std::cerr << it << " " << p << " " << s << " " << τ₀ << " " << y << " " << sqrt(2 * ΔC / C.size()) << " " << γ << std::endl;
+
     }
 
     Real e = energy(C, R, p, s, λ, y, Δτ);
@@ -160,11 +172,11 @@ int main(int argc, char* argv[]) {
     std::cerr << "y " << y << " " << e << " " << z << std::endl;
 
     std::ofstream outfile(fourierFile("C", p, s, λ, τ₀, y, log2n, τₘₐₓ), std::ios::out | std::ios::binary);
-    outfile.write((const char*)(C.data()), C.size() * sizeof(Real));
+    outfile.write((const char*)(C.data()), (C.size() / 2) * sizeof(Real));
     outfile.close();
 
     std::ofstream outfileR(fourierFile("R", p, s, λ, τ₀, y, log2n, τₘₐₓ), std::ios::out | std::ios::binary);
-    outfileR.write((const char*)(R.data()), R.size() * sizeof(Real));
+    outfileR.write((const char*)(R.data()), (R.size() / 2) * sizeof(Real));
     outfileR.close();
   }