1 files changed, 20 insertions, 19 deletions

diff --git a/log-fourier_integrator.cpp b/log-fourier_integrator.cpp
index 7b24a55..9db21f0 100644
--- a/log-fourier_integrator.cpp
+++ b/log-fourier_integrator.cpp
@@ -13,7 +13,7 @@ int main(int argc, char* argv[]) {
   unsigned log2n = 8;
   Real Δτ = 0.1;
   Real k = -0.01;
-  Real shift = 0.5;
+  Real logShift = 0;
 
   /* Iteration parameters */
   Real ε = 1e-15;
@@ -54,12 +54,12 @@ int main(int argc, char* argv[]) {
     case 'k':
       k = atof(optarg);
       break;
+    case 'h':
+      logShift = atof(optarg);
+      break;
     case 'D':
       Δτ = atof(optarg);
       break;
-    case 'h':
-      shift = atof(optarg);
-      break;
     case 'e':
       ε = atof(optarg);
       break;
@@ -85,13 +85,14 @@ int main(int argc, char* argv[]) {
 
   unsigned N = pow(2, log2n);
 
-  LogarithmicFourierTransform fft(N, k, Δτ, pad, shift);
-
   Real Γ₀ = 1;
-  Real μₜ₋₁ = Γ₀;
-  if (τ₀ > 0) {
-    μₜ₋₁ = (sqrt(1+4*Γ₀*τ₀)-1)/(2*τ₀);
-  }
+  Real μ₀ = τ₀ > 0 ? (sqrt(1+4*Γ₀*τ₀)-1)/(2*τ₀) : Γ₀;
+
+  LogarithmicFourierTransform fft(N, k, Δτ, pad, μ₀ * pow(10, logShift));
+
+  std::cerr << "Starting, μ₀ = " << μ₀ << ", range " << fft.t(0) << " " << fft.t(N-1) << std::endl;
+
+  Real μₜ₋₁ = μ₀;
 
   std::vector<Real> Cₜ₋₁(N);
   std::vector<Real> Rₜ₋₁(N);
@@ -103,20 +104,20 @@ int main(int argc, char* argv[]) {
     for (unsigned n = 0; n < N; n++) {
       if (τ₀ > 0) {
         if (τ₀ == 2) {
-          Cₜ₋₁[n] = Γ₀ * exp(-fft.t(n) / 2) * (1 + fft.t(n) / 2);
+          Cₜ₋₁[n] = Γ₀ * std::exp(-fft.t(n) / 2) * (1 + fft.t(n) / 2);
         } else {
-          Cₜ₋₁[n] = Γ₀ * (exp(-μₜ₋₁ * fft.t(n)) - μₜ₋₁ * τ₀ * exp(-fft.t(n) / τ₀)) / (μₜ₋₁ - pow(μₜ₋₁, 3) * pow(τ₀, 2));
+          Cₜ₋₁[n] = Γ₀ * (std::exp(-μ₀ * fft.t(n)) - μ₀ * τ₀ * std::exp(-fft.t(n) / τ₀)) / (μ₀ - pow(μ₀, 3) * pow(τ₀, 2));
         }
       } else {
-        Cₜ₋₁[n] = Γ₀ * exp(-μₜ₋₁ * fft.t(n)) / μₜ₋₁;
+        Cₜ₋₁[n] = Γ₀ * std::exp(-μ₀ * fft.t(n)) / μ₀;
       }
-      Rₜ₋₁[n] = exp(-μₜ₋₁ * fft.t(n));
+      Rₜ₋₁[n] = std::exp(-μ₀ * fft.t(n));
 
-      Ĉₜ₋₁[n] = 2 * Γ₀ / (pow(μₜ₋₁, 2) + pow(fft.ν(n), 2)) / (1 + pow(τ₀ * fft.ν(n), 2));
-      Ȓₜ₋₁[n] = (Real)1.0 / (μₜ₋₁ + II * fft.ν(n));
+      Ĉₜ₋₁[n] = 2 * Γ₀ / (pow(μ₀, 2) + pow(fft.ν(n), 2)) / (1 + pow(τ₀ * fft.ν(n), 2));
+      Ȓₜ₋₁[n] = (Real)1.0 / (μ₀ + II * fft.ν(n));
     }
   } else {
-    logFourierLoad(Cₜ₋₁, Rₜ₋₁, Ĉₜ₋₁, Ȓₜ₋₁, p, s, λ, τ₀, β₀, log2n, Δτ, shift);
+    logFourierLoad(Cₜ₋₁, Rₜ₋₁, Ĉₜ₋₁, Ȓₜ₋₁, p, s, λ, τ₀, β₀, log2n, Δτ, logShift);
     μₜ₋₁ = estimateZ(fft, Cₜ₋₁, Ĉₜ₋₁, Rₜ₋₁, Ȓₜ₋₁, p, s, λ, τ₀, β₀);
   }
 
@@ -139,7 +140,7 @@ int main(int argc, char* argv[]) {
       std::vector<Complex> Ȓₜ₊₁(N);
       for (unsigned n = 0; n < N; n++) {
         Ȓₜ₊₁[n] = ((Real)1.0 + std::pow(β, 2) * RddfCt[n] * Ȓₜ[n]) / (μₜ + II * fft.ν(n));
-        Ĉₜ₊₁[n] = (2 * Γ₀ * std::conj(Ȓₜ[n]) / (1 + std::pow(τ₀ * fft.ν(n), 2)) + std::pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ[n]))) / (μₜ + II * fft.ν(n));
+        Ĉₜ₊₁[n] = ((2 * Γ₀ * std::conj(Ȓₜ[n]) / (1 + std::pow(τ₀ * fft.ν(n), 2)) + std::pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ[n]))) / (μₜ + II * fft.ν(n))).real();
       }
       std::vector<Real> Rₜ₊₁ = fft.inverse(Ȓₜ₊₁);
       std::vector<Real> Cₜ₊₁ = fft.inverse(Ĉₜ₊₁);
@@ -190,7 +191,7 @@ int main(int argc, char* argv[]) {
       std::cerr << "\x1b[2K" << "\r";
       std::cerr << β << " " << μₜ << " " << Ĉₜ[0].real() << " " << E << " " << γ << std::endl;
 
-      logFourierSave(Cₜ, Rₜ, Ĉₜ, Ȓₜ, p, s, λ, τ₀, β, log2n, Δτ, shift);
+      logFourierSave(Cₜ, Rₜ, Ĉₜ, Ȓₜ, p, s, λ, τ₀, β, log2n, Δτ, logShift);
 
       β += Δβ;
       Cₜ₋₁ = Cₜ;