1 files changed, 55 insertions, 20 deletions

diff --git a/log-fourier_integrator.cpp b/log-fourier_integrator.cpp
index 81e7634..a032401 100644
--- a/log-fourier_integrator.cpp
+++ b/log-fourier_integrator.cpp
@@ -8,6 +8,26 @@ std::string logFourierFile(std::string prefix, unsigned p, unsigned s, Real λ,
   return prefix + "_" + std::to_string(p) + "_" + std::to_string(s) + "_" + std::to_string(λ) + "_" + std::to_string(τ₀) + "_" + std::to_string(β) + "_" + std::to_string(log2n) + "_" + std::to_string(Δτ)  + "_" + std::to_string(k) + ".dat";
 }
 
+std::tuple<std::vector<Complex>, std::vector<Complex>> RddfCtdfCt(LogarithmicFourierTransform& fft, const std::vector<Real>& C, const std::vector<Real>& R, unsigned p, unsigned s, Real λ) {
+  std::vector<Real> dfC(C.size());
+  std::vector<Real> RddfC(C.size());
+  for (unsigned n = 0; n < C.size(); 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);
+
+  return {RddfCt, dfCt};
+}
+
+Real estimateZ(LogarithmicFourierTransform& fft, const std::vector<Real>& C, const std::vector<Complex>& Ct, const std::vector<Real>& R, const std::vector<Complex>& Rt, unsigned p, unsigned s, Real λ, Real τ₀, Real β) {
+  auto [RddfCt, dfCt] = RddfCtdfCt(fft, C, R, p, s, λ);
+  Real Γ₀ = 1 + τ₀ / 2;
+
+  return ((2 * Γ₀ * std::conj(Rt[0]) + pow(β, 2) * (RddfCt[0] * Ct[0] + dfCt[0] * std::conj(Rt[0]))) / Ct[0]).real();
+}
+
 int main(int argc, char* argv[]) {
   /* Model parameters */
   unsigned p = 2;
@@ -23,12 +43,14 @@ int main(int argc, char* argv[]) {
   /* Iteration parameters */
   Real ε = 1e-14;
   Real γ = 1;
+  Real β₀ = 0;
   Real βₘₐₓ = 0.7;
   Real Δβ = 0.01;
+  bool loadData = false;
 
   int opt;
 
-  while ((opt = getopt(argc, argv, "p:s:2:T:t:b:d:g:k:D:e:")) != -1) {
+  while ((opt = getopt(argc, argv, "p:s:2:T:t:b:d:g:k:D:e:0:l")) != -1) {
     switch (opt) {
     case 'p':
       p = atoi(optarg);
@@ -60,6 +82,12 @@ int main(int argc, char* argv[]) {
     case 'e':
       ε = atof(optarg);
       break;
+    case '0':
+      β₀ = atof(optarg);
+      break;
+    case 'l':
+      loadData = true;
+      break;
     default:
       exit(1);
     }
@@ -77,17 +105,31 @@ int main(int argc, char* argv[]) {
   std::vector<Complex> Ĉₜ₋₁(N);
   std::vector<Complex> Ȓₜ₋₁(N);
 
-  /* Start from the exact solution for β = 0 */
-  for (unsigned n = 0; n < N; n++) {
-    if (τ₀ != 1) {
-      Cₜ₋₁[n] = Γ₀ * (exp(-μₜ₋₁ * fft.t(n)) - μₜ₋₁ * τ₀ * exp(-fft.t(n) / τ₀)) / (μₜ₋₁ - pow(μₜ₋₁, 3) * pow(τ₀, 2));
-    } else {
-      Cₜ₋₁[n] = Γ₀ * exp(-fft.t(n)) * (1 + fft.t(n));
-    }
-    Rₜ₋₁[n] = exp(-μₜ₋₁ * fft.t(n));
+  if (!loadData) {
+    /* Start from the exact solution for β = 0 */
+    for (unsigned n = 0; n < N; n++) {
+      if (τ₀ != 1) {
+        Cₜ₋₁[n] = Γ₀ * (exp(-μₜ₋₁ * fft.t(n)) - μₜ₋₁ * τ₀ * exp(-fft.t(n) / τ₀)) / (μₜ₋₁ - pow(μₜ₋₁, 3) * pow(τ₀, 2));
+      } else {
+        Cₜ₋₁[n] = Γ₀ * exp(-fft.t(n)) * (1 + fft.t(n));
+      }
+      Rₜ₋₁[n] = exp(-μₜ₋₁ * fft.t(n));
 
-    Ĉₜ₋₁[n] = 2 * Γ₀ / (pow(μₜ₋₁, 2) + pow(fft.ν(n), 2)) / (1 + pow(τ₀ * fft.ν(n), 2));
-    Ȓₜ₋₁[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));
+    }
+  } else {
+    std::ifstream cfile(logFourierFile("C", p, s, λ, τ₀, β₀, log2n, Δτ, k), std::ios::binary);
+    cfile.read((char*)(Cₜ₋₁.data()), N * sizeof(Real));
+    cfile.close();
+    std::ifstream rfile(logFourierFile("R", p, s, λ, τ₀, β₀, log2n, Δτ, k), std::ios::binary);
+    rfile.read((char*)(Rₜ₋₁.data()), N * sizeof(Real));
+    rfile.close();
+
+    Ĉₜ₋₁ = fft.fourier(Cₜ₋₁, true);
+    Ȓₜ₋₁ = fft.fourier(Rₜ₋₁, false);
+
+    μₜ₋₁ = estimateZ(fft, Cₜ₋₁, Ĉₜ₋₁, Rₜ₋₁, Ȓₜ₋₁, p, s, λ, τ₀, β₀);
   }
 
   std::vector<Real> Cₜ = Cₜ₋₁;
@@ -96,18 +138,11 @@ int main(int argc, char* argv[]) {
   std::vector<Complex> Ȓₜ = Ȓₜ₋₁;
   Real μₜ = μₜ₋₁;
 
-  Real β = 0;
+  Real β = β₀ + Δβ;
   while (β < βₘₐₓ) {
     Real ΔC = 100;
     while (ΔC > ε) {
-      std::vector<Real> dfC(N);
-      std::vector<Real> RddfC(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);
+      auto [RddfCt, dfCt] = RddfCtdfCt(fft, Cₜ, Rₜ, p, s, λ);
 
       std::vector<Complex> Ĉₜ₊₁(N);
       std::vector<Complex> Ȓₜ₊₁(N);