From 734dea28279ce7272e6ed23f821c69c225e2d947 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Fri, 14 Feb 2020 14:22:59 -0500
Subject: Split up Ising-related functions and added an executable for
 measuring Ising energy autocorrelation time.

---
 quantity.hpp | 47 ++++++++++++++++++++++-------------------------
 1 file changed, 22 insertions(+), 25 deletions(-)

(limited to 'quantity.hpp')

diff --git a/quantity.hpp b/quantity.hpp
index 55d5a28..d2a7dea 100644
--- a/quantity.hpp
+++ b/quantity.hpp
@@ -1,57 +1,54 @@
 
 #pragma once
 
-#include <vector>
-#include <list>
-#include <cmath>
 #include <array>
+#include <cmath>
+#include <list>
+#include <vector>
 
 class Quantity {
 private:
   double total;
   double total2;
   std::vector<double> C;
-  unsigned wait;
 
 public:
   unsigned n;
   std::list<double> hist;
 
-  Quantity(unsigned lag, unsigned wait) : C(lag), wait(wait) {
+  Quantity(unsigned lag) : C(lag) {
     n = 0;
     total = 0;
     total2 = 0;
   }
 
   void add(double x) {
-    if (n > wait) {
-      hist.push_front(x);
-      if (hist.size() > C.size()) {
-        hist.pop_back();
-        unsigned t = 0;
-        for (double a : hist) {
-          C[t] += a * x;
-          t++;
-        }
-        total += x;
-        total2 += pow(x, 2);
+    hist.push_front(x);
+    if (hist.size() > C.size()) {
+      hist.pop_back();
+      unsigned t = 0;
+      for (double a : hist) {
+        C[t] += a * x;
+        t++;
       }
+      total += x;
+      total2 += pow(x, 2);
+      n++;
     }
-    n++;
   }
 
-  double avg() const { return total / (n - wait); }
+  double avg() const { return total / n; }
 
-  double avg2() const { return total2 / (n - wait); }
+  double avg2() const { return total2 / n; }
 
   std::vector<double> ρ() const {
-    double C0 = C.front() / (n - wait);
-    double avg2 = pow(total / (n - wait), 2);
+    double C0 = C.front() / n;
+    double avg2 = pow(total / n, 2);
 
     std::vector<double> ρtmp;
 
     for (double Ct : C) {
-      ρtmp.push_back((Ct / (n - wait) - avg2) / (C0 - avg2));
+      ρtmp.push_back((Ct / n - avg2) / (C0 - avg2));
     }
 
     return ρtmp;
@@ -69,15 +66,15 @@ public:
       M++;
     }
 
-    return {τtmp, 2.0 * (2.0 * M + 1) * pow(τtmp, 2) / (n - wait)};
+    return {τtmp, 2.0 * (2.0 * M + 1) * pow(τtmp, 2) / n};
   }
 
   double σ() const {
-    return 2.0 / (n - wait) * this->τ()[0] * (C[0] / (n - wait) - pow(this->avg(), 2));
+    return 2.0 / n * this->τ()[0] * (C[0] / n - pow(this->avg(), 2));
   }
 
   double serr() const { return sqrt(this->σ()); }
 
-  unsigned num_added() const { return n - wait; }
+  unsigned num_added() const { return n; }
 };
 
-- 
cgit v1.2.3-54-g00ecf