1 files changed, 106 insertions, 0 deletions
diff --git a/ising.hpp b/ising.hpp
new file mode 100644
index 0000000..6f140be
--- /dev/null
+++ b/ising.hpp
@@ -0,0 +1,106 @@
+
+#include "smiley.hpp"
+#include "space_wolff.hpp"
+#include "torus_symmetries.hpp"
+
+const unsigned D = 2;
+typedef Model<signed, D, TorusGroup<signed, D>, signed> isingModel;
+typedef Spin<signed, D, signed> isingSpin;
+
+std::function<double(const isingSpin&, const isingSpin&)> isingZ(unsigned L) {
+  return [L](const isingSpin& s1, const isingSpin& s2) -> double {
+    bool old_one_one = false;
+    bool old_many_ones = false;
+    bool old_any_two = false;
+
+    Vector<signed, D> old_diff = diff<signed, D>(L, s1.x, s2.x);
+
+    for (unsigned i = 0; i < D; i++) {
+      if (old_diff(i) == 1 && !old_one_one) {
+        old_one_one = true;
+      } else if (old_diff(i) == 1 && old_one_one) {
+        old_many_ones = true;
+      } else if (old_diff(i) > 1) {
+        old_any_two = true;
+      }
+    }
+
+    bool were_on_someone = !old_one_one && !old_any_two;
+    bool were_nearest_neighbors = old_one_one && !old_many_ones && !old_any_two;
+
+    if (were_on_someone) {
+      return -std::numeric_limits<double>::infinity();
+      ;
+    } else if (were_nearest_neighbors) {
+      return s1.s * s2.s;
+    } else {
+      return 0.0;
+    }
+  };
+}
+
+std::function<double(isingSpin)> isingBFace(unsigned L, double H) {
+  return [L, H](const isingSpin& s) -> double {
+    return H * s.s * smiley[15 - s.x(1) * 16 / L][s.x(0) * 16 / L];
+  };
+}
+
+std::function<double(isingSpin)> isingBMod(unsigned L, unsigned mod, double H) {
+  return [mod, L, H](const isingSpin& s) -> double {
+    return H * s.s * sin(s.x(0) * mod * 2 * M_PI / L);
+  };
+}
+
+Gen<signed, D, TorusGroup<signed, D>, signed> isingGen(unsigned L) {
+  std::vector<Vector<signed, 2>> torusVectors = torus_vecs<signed, 2>(L);
+  std::vector<Matrix<signed, 2>> torusMatrices = torus_mats<signed, 2>();
+  return [L, torusVectors,
+          torusMatrices](Model<signed, D, TorusGroup<signed, D>, signed>& M,
+                         Rng& r) -> Transformation<signed, D, TorusGroup<signed, D>, signed>* {
+    Matrix<signed, 2> m;
+    Vector<signed, 2> t;
+
+    m = r.pick(torusMatrices);
+    t(0) = r.uniform<signed>(0, L);
+    t(1) = r.uniform<signed>(0, L);
+    t = t - m * t;
+
+    TorusGroup<signed, 2> g = TorusGroup<signed, 2>({(signed)L, t, m});
+
+    Spin<signed, 2, signed>* ss = r.pick(M.s);
+    Spin<signed, 2, signed> s2 = g.act(*ss);
+
+    while (ss->x(0) == s2.x(0) && ss->x(1) == s2.x(1)) {
+      ss = r.pick(M.s);
+      s2 = g.act(*ss);
+    }
+
+    std::set<Spin<signed, 2, signed>*> s2s = M.dict.at(s2.x);
+
+    if (s2s.empty()) {
+      return new SpinFlip<signed, 2, TorusGroup<signed, 2>, signed>(M, g, ss);
+    } else {
+      return new PairFlip<signed, 2, TorusGroup<signed, 2>, signed>(M, g, ss, *s2s.begin());
+    }
+  };
+}
+
+void isingPopulate(isingModel& m, unsigned L, double pop, double mag) {
+  Rng rng;
+
+  for (unsigned i = 0; i < L; i++) {
+    for (unsigned j = 0; j < L; j++) {
+      if (rng.uniform<double>(0, 1) < pop) {
+      Spin<signed, D, signed>* ss = new Spin<signed, D, signed>();
+      ss->x = {i, j};
+      if (rng.uniform<double>(0, 1) < mag) {
+        ss->s = 1;
+      } else {
+        ss->s = -1;
+      }
+      m.s.push_back(ss);
+      m.dict.insert(ss);
+      }
+    }
+  }
+}