Moved some repeated code into the pSpinModel class.

1 files changed, 6 insertions, 18 deletions

diff --git a/tensor.hpp b/tensor.hpp
index fc99042..1212bda 100644
--- a/tensor.hpp
+++ b/tensor.hpp
@@ -5,6 +5,8 @@
 
 #include <eigen3/unsupported/Eigen/CXX11/Tensor>
 
+#include "factorial.hpp"
+
 template <class Scalar>
 using Vector = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>;
 
@@ -89,25 +91,11 @@ Tensor<Scalar, p> generateCouplings(unsigned N, Distribution d, Generator& r) {
   return J;
 }
 
-template <class Scalar, int p, class Distribution, class Generator>
-Tensor<Scalar, p> plantState(const Tensor<Scalar, p>& J, const Vector<Scalar>& z, double β) {
-  Tensor<Scalar, p> JPlanted = J;
-
-  std::function<void(Tensor<Scalar, p>&, std::array<unsigned, p>)> plant =
-    [&z, β] (Tensor<Scalar, p>& JJ, std::array<unsigned, p> ind) {
-      Scalar Ji = getJ<Scalar, p>(JJ, ind);
-      Scalar zzz = 1;
-
-      for (unsigned i : ind) {
-        zzz *= z(i);
-      }
-
-      setJ<Scalar, p>(JJ, ind, Ji - β * zzz / pow(zzz.size(), p));
-    };
-
-  iterateOver<Scalar, p>(JPlanted, plant);
+template <class Real, int p, class Generator>
+Tensor<Real, p> generateRealPSpinCouplings(unsigned N, Generator& r) {
+  Real σp = sqrt(factorial(p) / ((Real)2 * pow(N, p - 1)));
 
-  return JPlanted;
+  return generateCouplings<Real, p>(N, std::normal_distribution<Real>(0, σp), r);
 }
 
 template <class Scalar>