1 files changed, 18 insertions, 3 deletions

diff --git a/hadamard.cpp b/hadamard.cpp
index d9e5bfd..86c27d8 100644
--- a/hadamard.cpp
+++ b/hadamard.cpp
@@ -1,5 +1,6 @@
 
 #include "hadamard_pt.hpp"
+#include "matrices.hpp"
 
 #include <fstream>
 #include <iostream>
@@ -9,6 +10,7 @@ class MeasureEnergy : public Measurement {
 public:
   unsigned N;
   double totalE;
+  double totalE2;
   unsigned n;
   std::vector<unsigned> ρ_dist;
 
@@ -16,11 +18,13 @@ public:
     n = n_bins;
     N = 0;
     totalE = 0;
+    totalE2 = 0;
   }
 
   void after_sweep(double, double E, const Orthogonal& M) override {
     N++;
     totalE += E;
+    totalE2 += pow(E, 2);
     double max = sqrt(M.size());
     for (unsigned i = 0; i < M.size(); i++) {
       for (unsigned j = 0; j < M.size(); j++) {
@@ -30,6 +34,7 @@ public:
   }
 
   double energy() const { return totalE / N; }
+  double specific_heat() const { return totalE2 / N - pow(totalE / N, 2); }
 };
 
 class MeasureTransitionRates : public ParallelMeasurement {
@@ -65,6 +70,10 @@ int main(int argc, char* argv[]) {
     switch (opt) {
     case 'k':
       k = atoi(optarg);
+      if (k == 0 || k > 8) {
+        std::cout << "The size k must be an integer from 1 to 8!" << std::endl;
+        exit(1);
+      }
       break;
     case 'b':
       β₀ = atof(optarg);
@@ -95,7 +104,7 @@ int main(int argc, char* argv[]) {
     }
   }
 
-  unsigned n = pow(2, k);
+  unsigned n = 4 * k;
 
   std::vector<Measurement*> As(N);
   for (Measurement*& A : As) {
@@ -106,7 +115,7 @@ int main(int argc, char* argv[]) {
   PT p(β₀, β₁, N, n, B, As);
 
   for (MCMC& sim : p.Ms) {
-    sim.M = walsh(k);
+    sim.M = hadamards[k - 1];
     sim.E = sim.M.energy();
   }
 
@@ -141,7 +150,13 @@ int main(int argc, char* argv[]) {
   file << std::endl;
 
   for (unsigned i = 0; i < As.size(); i++) {
-    file << std::fixed << ((MeasureEnergy*)As[i])->totalE / ((MeasureEnergy*)As[i])->N << " ";
+    file << std::fixed << ((MeasureEnergy*)As[i])->energy() << " ";
+  }
+
+  file << std::endl;
+
+  for (unsigned i = 0; i < As.size(); i++) {
+    file << std::fixed << ((MeasureEnergy*)As[i])->specific_heat() << " ";
   }
 
   file << std::endl;