Added file for calculating energy later.

2 files changed, 114 insertions, 1 deletions

diff --git a/Makefile b/Makefile
index 21f3def..c013825 100644
--- a/Makefile
+++ b/Makefile
@@ -1,4 +1,4 @@
-all: walk correlations integrator fourier_integrator
+all: walk correlations integrator fourier_integrator get_energy
 
 CC := clang++ -std=c++17 -Wno-mathematical-notation-identifier-extension -O3 -march=native -mtune=native -fopenmp -lfftw3 -g
 
@@ -13,3 +13,6 @@ integrator: integrator.cpp
 
 fourier_integrator: fourier_integrator.cpp
 	$(CC) fourier_integrator.cpp -o fourier_integrator
+
+get_energy: get_energy.cpp
+	$(CC) get_energy.cpp -o get_energy
diff --git a/get_energy.cpp b/get_energy.cpp
new file mode 100644
index 0000000..503823d
--- /dev/null
+++ b/get_energy.cpp
@@ -0,0 +1,110 @@
+#include <getopt.h>
+#include <vector>
+#include <cmath>
+#include <iostream>
+#include <fftw3.h>
+#include <complex>
+#include <fstream>
+
+using Real = double;
+using Complex = std::complex<Real>;
+using namespace std::complex_literals;
+
+inline Real fP(unsigned p, Real q) {
+  return 0.5 * pow(q, p);
+}
+
+inline Real dfP(unsigned p, Real q) {
+  return 0.5 * p * pow(q, p - 1);
+}
+
+inline Real ddfP(unsigned p, Real q) {
+  return 0.5 * p * (p - 1) * pow(q, p - 2);
+}
+
+inline Real f(Real λ, unsigned p, unsigned s, Real q) {
+  return (1 - λ) * fP(p, q) + λ * fP(s, q);
+}
+
+inline Real df(Real λ, unsigned p, unsigned s, Real q) {
+  return (1 - λ) * dfP(p, q) + λ * dfP(s, q);
+}
+
+inline Real ddf(Real λ, unsigned p, unsigned s, Real q) {
+  return (1 - λ) * ddfP(p, q) + λ * ddfP(s, q);
+}
+
+int main(int argc, char* argv[]) {
+  unsigned p = 2;
+  unsigned s = 2;
+  Real λ = 0.5;
+  Real τ₀ = 0;
+  Real y₀ = 0;
+  Real yₘₐₓ = 0.5;
+  Real Δy = 0.05;
+
+  unsigned log2n = 8;
+  Real τₘₐₓ = 20;
+
+  int opt;
+
+  while ((opt = getopt(argc, argv, "p:s:2:T:t:0:y:d:")) != -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 't':
+      τ₀ = atof(optarg);
+      break;
+    case '0':
+      y₀ = atof(optarg);
+      break;
+    case 'y':
+      yₘₐₓ = atof(optarg);
+      break;
+    case 'd':
+      Δy = atof(optarg);
+      break;
+    default:
+      exit(1);
+    }
+  }
+
+  unsigned n = pow(2, log2n);
+
+  Real Δτ = τₘₐₓ / M_PI / n;
+
+  Real y = y₀;
+
+  while (y += Δy, y <= yₘₐₓ) {
+    std::vector<Real> C(2 * n);
+    std::vector<Real> R(2 * n);
+
+    std::string file_end = std::to_string(p) + "_" + std::to_string(s) + "_" + std::to_string(λ) + "_" + std::to_string(τ₀) + "_" + std::to_string(y) + "_" + std::to_string(log2n) + "_" + std::to_string(τₘₐₓ) + ".dat";
+    std::ifstream cfile("C_"+file_end, std::ios::binary);
+    cfile.read((char*)(C.data()), C.size() * sizeof(Real));
+    cfile.close();
+    std::ifstream rfile("R_"+file_end, std::ios::binary);
+    rfile.read((char*)(R.data()), R.size() * sizeof(Real));
+    rfile.close();
+
+    Real energy = 0;
+
+    for (unsigned i = 0; i < n; i++) {
+      energy += y * R[i] * df(λ, p, s, C[i]) * M_PI * Δτ;
+    }
+
+    std::cerr << y << " " << energy << std::endl;
+  }
+
+  return 0;
+}