1 files changed, 76 insertions, 0 deletions
diff --git a/correlations.cpp b/correlations.cpp
new file mode 100644
index 0000000..0d48844
--- /dev/null
+++ b/correlations.cpp
@@ -0,0 +1,76 @@
+#include <getopt.h>
+#include <fstream>
+#include <iostream>
+
+#include "eigen/Eigen/Dense"
+
+using Real = float;
+using Vector = Eigen::Matrix<Real, Eigen::Dynamic, 1>;
+
+int main(int argc, char* argv[]) {
+  unsigned N = 10;
+  Real E = 0;
+  Real Δt = 1e-4;
+  Real Δw = 1e-2;
+  Real T0 = 0;
+  Real T = 100;
+
+  int opt;
+
+  while ((opt = getopt(argc, argv, "N:E:t:w:0:T:")) != -1) {
+    switch (opt) {
+    case 'N':
+      N = (unsigned)atof(optarg);
+      break;
+    case 'E':
+      E = atof(optarg);
+      break;
+    case 't':
+      Δt = atof(optarg);
+      break;
+    case 'w':
+      Δw = atof(optarg);
+      break;
+    case '0':
+      T0 = atof(optarg);
+      break;
+    case 'T':
+      T = atof(optarg);
+      break;
+    default:
+      exit(1);
+    }
+  }
+
+  std::string filebase = std::to_string(N) + "_" + std::to_string(E) + "_" + std::to_string(-std::log10(Δt)) + "_" + std::to_string(Δw);
+
+  std::ifstream file(filebase + ".dat", std::ios::binary|std::ios::in|std::ios::ate);
+
+  unsigned size = file.tellg() / (N * sizeof(float)) - 1;
+  unsigned i0 = T0 / Δw;
+
+  if (i0 > size) {
+    return 0;
+  }
+
+  file.seekg((i0 + 1) * sizeof(float) * N);
+
+  Vector x0(N);
+  file.read(reinterpret_cast<char*>(x0.data()), N * sizeof(float));
+
+  file.seekg((i0 + 1) * sizeof(float) * N);
+
+  for (Real t = 0; t <= T; t += Δw) {
+    if (file.peek() == EOF) {
+      break;
+    }
+
+    Vector x(N);
+    file.read(reinterpret_cast<char*>(x.data()), N * sizeof(float));
+
+    std::cout << t << "\t" << x0.dot(x) / N << std::endl;
+  }
+
+
+  return 0;
+}