1 files changed, 18 insertions, 11 deletions

diff --git a/fourier_integrator.cpp b/fourier_integrator.cpp
index d83fb9d..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,16 +98,19 @@ 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);
 
-    z = estimateZ(fft, C, Ct, R, Rt, p, s, λ, y);
+    z = estimateZ(fft, C, Ct, R, Rt, p, s, λ, τ₀, y);
   }
 
   while (y += Δy, y <= yₘₐₓ) {
@@ -165,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();
   }