1 files changed, 100 insertions, 0 deletions
diff --git a/log_get_energy.cpp b/log_get_energy.cpp
new file mode 100644
index 0000000..d156fd4
--- /dev/null
+++ b/log_get_energy.cpp
@@ -0,0 +1,100 @@
+#include "log-fourier.hpp"
+
+#include <getopt.h>
+#include <iostream>
+#include <iomanip>
+#include <filesystem>
+
+int main(int argc, char* argv[]) {
+  /* Model parameters */
+  unsigned p = 2;
+  unsigned s = 2;
+  Real λ = 0.5;
+  Real τ₀ = 0;
+
+  /* Log-Fourier parameters */
+  unsigned log2n = 8;
+  Real Δτ = 0.1;
+  Real k = 0.1;
+  unsigned pad = 2;
+  Real logShift = 0;
+  bool shiftSquare = false;
+
+  /* Iteration parameters */
+  Real β₀ = 0;
+  Real βₘₐₓ = 0.7;
+  Real Δβ = 0.01;
+
+  int opt;
+
+  while ((opt = getopt(argc, argv, "p:s:2:T:t:b:d:k:h:D:0:S")) != -1) {
+    switch (opt) {
+    case 'p':
+      p = atoi(optarg);
+      break;
+    case 's':
+      s = atoi(optarg);
+      break;
+    case '2':
+      log2n = atoi(optarg);
+      break;
+    case 't':
+      τ₀ = atof(optarg);
+      break;
+    case 'b':
+      βₘₐₓ = atof(optarg);
+      break;
+    case 'd':
+      Δβ = atof(optarg);
+      break;
+    case 'k':
+      k = atof(optarg);
+      break;
+    case 'h':
+      logShift = atof(optarg);
+      break;
+    case 'D':
+      Δτ = atof(optarg);
+      break;
+    case '0':
+      β₀ = atof(optarg);
+      break;
+    case 'S':
+      shiftSquare = true;
+      break;
+    default:
+      exit(1);
+    }
+  }
+
+  unsigned N = pow(2, log2n);
+  Real Γ₀ = 1;
+  Real μ₀ = τ₀ > 0 ? (sqrt(1+4*Γ₀*τ₀)-1)/(2*τ₀) : Γ₀;
+
+  Real shift = μ₀ * pow(10, logShift);
+  if (shiftSquare) shift *= μ₀;
+  LogarithmicFourierTransform fft(N, k, Δτ, pad, shift);
+
+  std::vector<Real> C(N);
+  std::vector<Real> R(N);
+  std::vector<Complex> Ĉ(N);
+  std::vector<Complex> Ȓ(N);
+
+  Real β = β₀;
+
+  std::cout << std::setprecision(16);
+
+  while (β = std::round(1e6 * (β + Δβ)) / 1e6, β <= βₘₐₓ) {
+    if (std::filesystem::exists(logFourierFile("C", p, s, λ, τ₀, β, log2n, Δτ, logShift))) {
+      logFourierLoad(C, R, Ĉ, Ȓ, p, s, λ, τ₀, β, log2n, Δτ, logShift);
+
+      Real e = energy(fft, C, R, p, s, λ, β);
+
+      Real μ = estimateZ(fft, C, Ĉ, R, Ȓ, p, s, λ, τ₀, β);
+
+      std::cout << p << " " << s << " " << λ << " " << τ₀ << " " << β << " " << μ << " " << Ĉ[0].real() << " " << e << " " << Ȓ[0].real() << std::endl;
+    }
+  }
+
+  return 0;
+}