Minor refactoring of sphere commands and animation class, and introduction of Dimers

8 files changed, 400 insertions, 317 deletions

diff --git a/animation.hpp b/animation.hpp
new file mode 100644
index 0000000..857354e
--- /dev/null
+++ b/animation.hpp
@@ -0,0 +1,111 @@
+
+#include "space_wolff.hpp"
+#include <GL/glut.h>
+
+template <typename T, int D>
+void drawSphere(const Vector<T, D>& x, Radius r, std::array<double, 3> color = {0, 0, 0},
+                unsigned nPoints = 20) {
+  glColor3d(color[0], color[1], color[2]);
+  glBegin(GL_POLYGON);
+  for (unsigned i = 0; i < nPoints; i++) {
+    glVertex2d(x(0) + r * cos(2 * i * M_PI / nPoints), x(1) + r * sin(2 * i * M_PI / nPoints));
+  }
+  glEnd();
+}
+
+template <typename T, int D>
+void draw(const Spin<T, D, Radius>& s, std::array<double, 3> color = {0, 0, 0}) {
+  drawSphere<T, D>(s.x, s.s, color);
+}
+
+template <typename T, int D>
+void draw(const Spin<T, D, Dimer<T, D>>& s, std::array<double, 3> color = {0, 0, 0}) {
+  drawSphere<T, D>(s.x + s.s.relativePosition, s.s.radius, color);
+  drawSphere<T, D>(s.x - s.s.relativePosition, s.s.radius, color);
+}
+
+template <int D>
+void draw(const Spin<signed, D, IsingSpin>& s, std::array<double, 3> color = {0, 0, 0}) {
+  if (s.s == 1) {
+    glColor3d(color[0], color[1], color[2]);
+  } else if (s.s == -1) {
+    glColor3d(1 - color[0], 1 - color[1], 1 - color[2]);
+  }
+  glRecti(s.x(0), s.x(1), s.x(0) + 1, s.x(1) + 1);
+}
+
+template <typename T, int D, class R, class S> class Animation : public measurement<T, D, R, S> {
+private:
+  unsigned n;
+  unsigned wait;
+  double L;
+  R s₀⁻¹;
+
+public:
+  Animation(T L, unsigned w, int argc, char* argv[], unsigned wait, bool periodic = false)
+      : s₀⁻¹(0), wait(wait), L(1000 * L) {
+    n = 0;
+
+    glutInit(&argc, argv);
+    glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);
+    glutInitWindowSize(w, w);
+    glutCreateWindow("wolffWindow");
+    glClearColor(0.0, 0.0, 0.0, 0.0);
+    glMatrixMode(GL_PROJECTION);
+    glLoadIdentity();
+    if (periodic) {
+      gluOrtho2D(0, L, 0, L);
+    } else {
+      gluOrtho2D(-L, L, -L, L);
+    }
+  }
+
+  void pre_cluster(const Model<T, D, R, S>& m, unsigned,
+                   const Transformation<T, D, R, S>* t) override {
+    s₀⁻¹ = m.s0.inverse();
+    glClearColor(1.0f, 1.0f, 1.0f, 1.0f);
+    glClear(GL_COLOR_BUFFER_BIT);
+    for (const Spin<T, D, S>* s : m.s) {
+      draw(s₀⁻¹.act(*s));
+    }
+    if (n > wait) {
+      glColor3d(1.0, 0.0, 0.0);
+      glLineWidth(3);
+      glBegin(GL_LINES);
+      Vector<double, D> r = (t->r.t).template cast<double>() / 2.0;
+      double θ = atan2(r(1), r(0));
+      Vector<double, D> v1 = s₀⁻¹.template act<double>({r(0) + L * sin(θ), r(1) - L * cos(θ)});
+      Vector<double, D> v2 = s₀⁻¹.template act<double>({r(0) - L * sin(θ), r(1) + L * cos(θ)});
+      glVertex2d(v1(0), v1(1));
+      glVertex2d(v2(0), v2(1));
+      glEnd();
+    }
+    glFlush();
+  }
+
+  void plain_site_transformed(const Model<T, D, R, S>& m,
+                              const Transformation<T, D, R, S>& t) override {
+    if (n > wait) {
+      std::array<double, 3> color;
+      if (t.current().size() > 1) {
+        color = {0, 0, 1};
+      } else {
+        color = {0, 1, 0};
+      }
+      for (const Spin<T, D, S>* s : t.current()) {
+        if (s != NULL) {
+          draw(s₀⁻¹.act(*s), color);
+        } else {
+        }
+      }
+    }
+  }
+
+  void post_cluster(const Model<T, D, R, S>& m) override {
+    if (n > wait) {
+      glFlush();
+      sleep(2);
+    }
+    n++;
+  }
+};
diff --git a/dimers.cpp b/dimers.cpp
new file mode 100644
index 0000000..4416333
--- /dev/null
+++ b/dimers.cpp
@@ -0,0 +1,122 @@
+
+#include <chrono>
+#include <fstream>
+#include <iostream>
+
+#include "animation.hpp"
+#include "space_wolff.hpp"
+#include "spheres.hpp"
+
+const unsigned D = 2;
+typedef Model<double, D, Euclidean<double, D>, Dimer<double, 2>> model;
+
+int main(int argc, char* argv[]) {
+  const unsigned D = 2;
+
+  double L = 32;
+  unsigned N = 1000;
+  double T = 2.0 / log(1.0 + sqrt(2.0));
+  double H = 1.0;
+  unsigned n = 25;
+  unsigned wait = 1000;
+
+  double k = 1e2;
+  double a = 0.1;
+
+  int opt;
+
+  while ((opt = getopt(argc, argv, "n:N:L:T:H:a:k:w:")) != -1) {
+    switch (opt) {
+    case 'n':
+      n = (unsigned)atof(optarg);
+      break;
+    case 'N':
+      N = (unsigned)atof(optarg);
+      break;
+    case 'L':
+      L = atof(optarg);
+      break;
+    case 'T':
+      T = atof(optarg);
+      break;
+    case 'H':
+      H = atof(optarg);
+      break;
+    case 'a':
+      a = atof(optarg);
+      break;
+    case 'k':
+      k = atof(optarg);
+      break;
+    case 'w':
+      wait = atoi(optarg);
+      break;
+    default:
+      exit(1);
+    }
+  }
+
+  auto zSingle = zSpheres<D>(a, k);
+
+  std::function<double(const Spin<double, D, Dimer<double, D>>&,
+                       const Spin<double, D, Dimer<double, D>>&)>
+      Z = [zSingle, a, k](const Spin<double, D, Dimer<double, D>>& s1,
+                          const Spin<double, D, Dimer<double, D>>& s2) -> double {
+    Spin<double, D, Radius> s11 = {.x = s1.x + s1.s.relativePosition, .s = s1.s.radius};
+    Spin<double, D, Radius> s12 = {.x = s1.x - s1.s.relativePosition, .s = s1.s.radius};
+    Spin<double, D, Radius> s21 = {.x = s2.x + s2.s.relativePosition, .s = s2.s.radius};
+    Spin<double, D, Radius> s22 = {.x = s2.x - s2.s.relativePosition, .s = s2.s.radius};
+
+    return zSingle(s11, s21) + zSingle(s12, s21) + zSingle(s11, s22) + zSingle(s12, s22);
+  };
+
+  auto g1 = nudgeGen<D, Dimer<double, D>>(1);
+  auto g2 = swapGen<D, Dimer<double, D>>(0.1);
+  auto g3 = accrossGen<D, Dimer<double, D>>(0.1);
+  auto g4 = centerGen<D, Dimer<double, D>>(0);
+
+  auto tag = std::chrono::high_resolution_clock::now();
+
+  std::string filename = "flips_" + std::to_string(n) + "_" + std::to_string(T) + "_" +
+                         std::to_string(H) + "_" + std::to_string(a) + "_" + std::to_string(k) +
+                         "_" + std::to_string(tag.time_since_epoch().count()) + ".dat";
+
+  std::ofstream file;
+  file.open(filename);
+
+  Animation<double, D, Euclidean<double, D>, Dimer<double, D>> A(L, 750, argc, argv, wait);
+  model sphere(1.0, Z, bCenter<D, Dimer<double, D>>(H));
+
+  Rng rng;
+
+  sphere.s.resize(n);
+
+  unsigned nx = floor(sqrt(n));
+  for (unsigned i = 0; i < sphere.s.size(); i++) {
+    Spin<double, 2, Dimer<double, D>>* ss = new Spin<double, 2, Dimer<double, D>>();
+    ss->x = {(i / nx) * L / nx - L / 2, (i % nx) * L / nx - L / 2};
+    ss->s = Dimer<double, D>{.relativePosition = {0.2, 0}, .radius = 0.25};
+    sphere.s[i] = ss;
+    sphere.dict.insert(ss);
+  }
+
+  sphere.wolff(T, {g1, g2, g3, g4}, A, N);
+
+  file.close();
+
+  /*
+  std::ofstream snapfile;
+  snapfile.open("sphere_snap.dat");
+
+  for (Spin<double, D, double>* s : sphere.s) {
+    Spin<double, D, double> rs = sphere.s0.inverse().act(*s);
+    snapfile << rs.s << " " << rs.x.transpose() << "\n";
+    delete s;
+  }
+
+  snapfile.close();
+  */
+
+  return 0;
+}
+
diff --git a/euclidean.hpp b/euclidean.hpp
index b1ce2c6..fb2dbcd 100644
--- a/euclidean.hpp
+++ b/euclidean.hpp
@@ -5,6 +5,15 @@
 #include "spin.hpp"
 #include "vector.hpp"
 
+typedef double Radius;
+typedef signed IsingSpin;
+
+template <class T, int D> class Dimer {
+public:
+  double radius;
+  Vector<T, D> relativePosition;
+};
+
 template <class T, int D> class Euclidean {
 public:
   Vector<T, D> t;
@@ -25,36 +34,27 @@ public:
     r = r0;
   }
 
-  Vector<T, D> act(const Vector<T, D>& x) const {
-    return t + r * x;
-  }
-
-  template <class S> Spin<T, D, S> act(const Spin<T, D, S>& s) const {
-    Spin<T, D, S> s_new;
-
-    s_new.x = this->act(s.x);
-    s_new.s = s.s;
+  Vector<T, D> act(const Vector<T, D>& x) const { return t + r * x; }
 
-    return s_new;
+  template <typename U> Vector<U, D> act(const Vector<U, D>& x) const {
+    return t.template cast<U>() + r.template cast<U>() * x;
   }
 
-  Euclidean act(const Euclidean& x) const {
-    Vector<T, D> tnew = r * x.t + t;
-    Matrix<T, D> rnew = r * x.r;
+  Radius act(Radius r) const { return r; }
 
-    Euclidean pnew(tnew, rnew);
+  IsingSpin act(IsingSpin s) const { return s; }
 
-    return pnew;
+  Dimer<T, D> act(const Dimer<T, D>& d) const {
+    return {.radius = d.radius, .relativePosition = r * d.relativePosition};
   }
 
-  Euclidean inverse() const {
-    Vector<T, D> tnew = -r.transpose() * t;
-    Matrix<T, D> rnew = r.transpose();
+  template <class S> Spin<T, D, S> act(const Spin<T, D, S>& s) const {
+    return {.x = act(s.x), .s = act(s.s)};
+  }
 
-    Euclidean pnew(tnew, rnew);
+  Euclidean act(const Euclidean& x) const { return Euclidean(r * x.t + t, r * x.r); }
 
-    return pnew;
-  }
+  Euclidean inverse() const { return Euclidean(-r.transpose() * t, r.transpose()); }
 };
 
 template <class T, int D> class TorusGroup {
@@ -105,6 +105,16 @@ public:
     return s_new;
   }
 
+  template <typename U> Vector<U, D> act(const Vector<U, D>& s) const {
+    Vector<U, D> s_new = t.template cast<U>() + r.template cast<U>() * s;
+
+    for (unsigned i = 0; i < D; i++) {
+      s_new(i) = fmod(L + s_new(i), L);
+    }
+
+    return s_new;
+  }
+
   TorusGroup act(const TorusGroup& x) const {
     Vector<T, D> tnew = r * x.t + t;
     Matrix<T, D> rnew = r * x.r;
diff --git a/ising.hpp b/ising.hpp
index 6f140be..a97da26 100644
--- a/ising.hpp
+++ b/ising.hpp
@@ -91,15 +91,15 @@ void isingPopulate(isingModel& m, unsigned L, double pop, double mag) {
   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);
+        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);
       }
     }
   }
diff --git a/ising_animate.cpp b/ising_animate.cpp
index cd65b6b..f3ab314 100644
--- a/ising_animate.cpp
+++ b/ising_animate.cpp
@@ -1,45 +1,5 @@
+#include "animation.hpp"
 #include "ising.hpp"
-#include <GL/glut.h>
-
-class Animation : public measurement<signed, D, TorusGroup<signed, D>, signed> {
-private:
-  bool color;
-
-public:
-  Animation(double L, unsigned w, bool tcolor, int argc, char* argv[]) {
-    color = tcolor;
-
-    glutInit(&argc, argv);
-    glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);
-    glutInitWindowSize(w, w);
-    glutCreateWindow("wolffWindow");
-    glClearColor(0.0, 0.0, 0.0, 0.0);
-    glMatrixMode(GL_PROJECTION);
-    glLoadIdentity();
-    gluOrtho2D(0, L, 0, L);
-  }
-
-  void post_cluster(const isingModel& m) override {
-    glClearColor(1.0, 1.0, 1.0, 1.0);
-    glClear(GL_COLOR_BUFFER_BIT);
-    for (const Spin<signed, 2, signed>* s : m.s) {
-      if (s->s == 1) {
-        if (color)
-          glColor3f(1.0, 0.0, 0.0);
-        else
-          glColor3f(1.0, 1.0, 1.0);
-      } else if (s->s == -1) {
-        if (color)
-          glColor3f(0.0, 0.0, 1.0);
-        else
-          glColor3f(0.0, 0.0, 0.0);
-      }
-      Vector<signed, 2> xx = m.s0.inverse().act(*s).x;
-      glRecti(xx(0), xx(1), xx(0) + 1, xx(1) + 1);
-    }
-    glFlush();
-  }
-};
 
 int main(int argc, char* argv[]) {
   unsigned L = 32;
@@ -50,10 +10,11 @@ int main(int argc, char* argv[]) {
   double T = 2.0 / log(1.0 + sqrt(2.0));
   double H = 1.0;
   bool color = false;
+  unsigned wait = N;
 
   int opt;
 
-  while ((opt = getopt(argc, argv, "N:L:T:H:m:r:p:c")) != -1) {
+  while ((opt = getopt(argc, argv, "N:L:T:H:m:r:p:cw:")) != -1) {
     switch (opt) {
     case 'N':
       N = (unsigned)atof(optarg);
@@ -79,6 +40,9 @@ int main(int argc, char* argv[]) {
     case 'c':
       color = true;
       break;
+    case 'w':
+      wait = atoi(optarg);
+      break;
     default:
       exit(1);
     }
@@ -97,7 +61,7 @@ int main(int argc, char* argv[]) {
 
   auto g = isingGen(L);
 
-  Animation A(L, 750, color, argc, argv);
+  Animation<signed, D, TorusGroup<signed, D>, IsingSpin> A(L, 750, argc, argv, wait, true);
 
   ising.wolff(T, {g}, A, N);
 
diff --git a/spheres.hpp b/spheres.hpp
new file mode 100644
index 0000000..040313c
--- /dev/null
+++ b/spheres.hpp
@@ -0,0 +1,112 @@
+
+#include "space_wolff.hpp"
+
+template <int D>
+std::function<double(const Spin<double, D, Radius>&, const Spin<double, D, Radius>&)>
+zSpheres(double a, double k) {
+  return [a, k](const Spin<double, D, Radius>& s1, const Spin<double, D, Radius>& s2) -> double {
+    Vector<double, D> d = s1.x - s2.x;
+
+    double σ = s1.s + s2.s;
+    double δ = σ - sqrt(d.transpose() * d);
+
+    if (δ > -a * σ) {
+      return 0.5 * k * (2 * pow(a * σ, 2) - pow(δ, 2));
+    } else if (δ > -2 * a * σ) {
+      return 0.5 * k * pow(δ + 2 * a * σ, 2);
+    } else {
+      return 0;
+    }
+  };
+}
+
+template <int D, class S> std::function<double(Spin<double, D, S>)> bCenter(double H) {
+  return [H](Spin<double, D, S> s) -> double { return H * s.x.norm(); };
+}
+
+template <int D> Euclidean<double, D> flipAround(double θ, Vector<double, D>& t₀) {
+  Matrix<double, D> m;
+
+  m(0, 0) = -cos(2 * θ);
+  m(1, 1) = cos(2 * θ);
+  m(0, 1) = -2 * cos(θ) * sin(θ);
+  m(1, 0) = -2 * cos(θ) * sin(θ);
+
+  return Euclidean<double, D>(t₀ - m * t₀, m);
+}
+
+template <int D, class S> Gen<double, D, Euclidean<double, D>, S> nudgeGen(double ε) {
+  return [ε](Model<double, D, Euclidean<double, D>, S>& M,
+             Rng& rng) -> Transformation<double, D, Euclidean<double, D>, S>* {
+    double θ = rng.uniform((double)0.0, 2 * M_PI);
+
+    Spin<double, D, S>* s = rng.pick(M.s);
+    Vector<double, D> t = s->x;
+    for (unsigned j = 0; j < D; j++) {
+      t(j) += rng.variate<double, std::normal_distribution>(0.0, ε);
+    }
+
+    return new SpinFlip<double, D, Euclidean<double, D>, S>(M, flipAround(θ, t), s);
+  };
+}
+
+template <int D, class S> Gen<double, D, Euclidean<double, D>, S> swapGen(double ε) {
+  return [ε](Model<double, D, Euclidean<double, D>, S>& M,
+             Rng& rng) -> Transformation<double, D, Euclidean<double, D>, S>* {
+    Spin<double, D, S>* s1 = rng.pick(M.s);
+    Spin<double, D, S>* s2 = rng.pick(M.s);
+
+    while (s1 == s2) {
+      s2 = rng.pick(M.s);
+    }
+
+    Vector<double, D> t1 = s1->x;
+    Vector<double, D> t2 = s2->x;
+    Vector<double, D> t12 = t1 - t2;
+    Vector<double, D> t = (t1 + t2) / 2;
+
+    double θ =
+        atan2(t12(1), t12(0)) + rng.variate<double, std::normal_distribution>(0.0, ε) / t12.norm();
+
+    return new PairFlip<double, D, Euclidean<double, D>, S>(M, flipAround(θ, t), s1, s2);
+  };
+}
+
+template <int D, class S> Gen<double, D, Euclidean<double, D>, S> accrossGen(double ε) {
+  return [ε](Model<double, D, Euclidean<double, D>, S>& M,
+             Rng& rng) -> Transformation<double, D, Euclidean<double, D>, S>* {
+    Spin<double, D, S>* s1 = rng.pick(M.s);
+    Spin<double, D, S>* s2 = rng.pick(M.s);
+
+    while (s1 == s2) {
+      s2 = rng.pick(M.s);
+    }
+
+    Vector<double, D> t1 = s1->x;
+    Vector<double, D> t2 = s2->x;
+    Vector<double, D> t12 = t1 - t2;
+    Vector<double, D> t = t2;
+
+    double θ =
+        atan2(t12(1), t12(0)) + rng.variate<double, std::normal_distribution>(0.0, ε) / t12.norm();
+
+    return new SpinFlip<double, D, Euclidean<double, D>, S>(M, flipAround(θ, t), s1);
+  };
+}
+
+template <int D, class S> Gen<double, D, Euclidean<double, D>, S> centerGen(double ε) {
+  return [ε](Model<double, D, Euclidean<double, D>, S>& M,
+             Rng& rng) -> Transformation<double, D, Euclidean<double, D>, S>* {
+    double θ = rng.uniform((double)0.0, 2 * M_PI);
+
+    Vector<double, D> t = M.s0.t;
+    for (unsigned j = 0; j < D; j++) {
+      t(j) += rng.variate<double, std::normal_distribution>(0.0, ε);
+    }
+
+    Spin<double, D, S>* s = rng.pick(M.s);
+
+    return new SpinFlip<double, D, Euclidean<double, D>, S>(M, flipAround(θ, t), s);
+  };
+}
+
diff --git a/spheres_infinite.cpp b/spheres_infinite.cpp
index 56abe77..af91427 100644
--- a/spheres_infinite.cpp
+++ b/spheres_infinite.cpp
@@ -4,226 +4,12 @@
 #include <chrono>
 
 #include "space_wolff.hpp"
-#include <GL/glut.h>
+#include "spheres.hpp"
+#include "animation.hpp"
 
 const unsigned D = 2;
 typedef Model<double, D, Euclidean<double, D>, double> model;
 
-class animation : public measurement<double, 2, Euclidean<double, D>, double> {
-private:
-  uint64_t t1;
-  uint64_t t2;
-  unsigned n;
-  unsigned tmp;
-  unsigned wait;
-  double L;
-  Euclidean<double, D> s0_tmp;
-  std::ofstream& file;
-
-public:
-  animation(double L, unsigned w, int argc, char* argv[], unsigned wait, std::ofstream& file) : s0_tmp(0), wait(wait), L(50 * L), file(file) {
-    t1 = 0;
-    t2 = 0;
-    n = 0;
-
-    glutInit(&argc, argv);
-    glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);
-    glutInitWindowSize(w, w);
-    glutCreateWindow("wolffWindow");
-    glClearColor(0.0, 0.0, 0.0, 0.0);
-    glMatrixMode(GL_PROJECTION);
-    glLoadIdentity();
-    gluOrtho2D(-L, L, -L, L);
-  }
-
-  void pre_cluster(const model& m, unsigned, const Transformation<double, D, Euclidean<double, D>, double>* t) override {
-    s0_tmp = m.s0.inverse();
-    tmp = 0;
-    glClearColor(1.0f, 1.0f, 1.0f, 1.0f);
-    glClear(GL_COLOR_BUFFER_BIT);
-    glColor3d(1.0, 0.0, 0.0);
-    glBegin(GL_LINES);
-      Vector<double, 2> r = t->r.t / 2;
-      double θ = atan2(r(1), r(0));
-      Vector<double, 2> v1 = s0_tmp.act({r(0) + L * sin(θ), r(1) - L * cos(θ)});
-      Vector<double, 2> v2 = s0_tmp.act({r(0) - L * sin(θ), r(1) + L * cos(θ)});
-      glVertex2d(v1(0), v1(1));
-      glVertex2d(v2(0), v2(1));
-      if (n > wait)
-        file << v1.transpose() << " " << v2.transpose() << std::endl;
-    glEnd();
-    for (const Spin<double, 2, double>* s : m.s) {
-      glBegin(GL_POLYGON);
-      unsigned n_points = 50;
-      glColor3f(0.0f, 0.0f, 0.0f);
-      if (n > wait)
-        file << s << " " << s->s << " " << s0_tmp.act(*s).x.transpose() << " ";
-      for (unsigned i = 0; i < n_points; i++) {
-        glVertex2d(s0_tmp.act(*s).x(0) + s->s * cos(2 * i * M_PI / n_points),
-                   s0_tmp.act(*s).x(1) + s->s * sin(2 * i * M_PI / n_points));
-      }
-      glEnd();
-    }
-    if (n > wait)
-      file << std::endl;
-    glLineWidth(3);
-    glFlush();
-  }
-
-  void plain_site_transformed(const model& m, const Transformation<double, D, Euclidean<double, D>, double>& t) override {
-    if (n > wait) {
-      if (t.current().size() > 1) {
-        glColor3f(0.0f, 0.0f, 1.0f);
-      } else {
-        glColor3f(0.0f, 1.0f, 0.0f);
-      }
-      for (const Spin<double, 2, double>* s : t.current()) {
-        if (s != NULL) {
-          file << s << " " << s0_tmp.act(t.r.act(*s)).x.transpose() << " ";
-          glBegin(GL_POLYGON);
-          unsigned n_points = 50;
-          for (unsigned i = 0; i < n_points; i++) {
-            glVertex2d(s0_tmp.act(*s).x(0) + s->s * cos(2 * i * M_PI / n_points),
-                       s0_tmp.act(*s).x(1) + s->s * sin(2 * i * M_PI / n_points));
-          }
-          glEnd();
-        } else {
-          file << s << " " << s0_tmp.act({0, 0}).transpose() << " ";
-        }
-      }
-    }
-    tmp++;
-  }
-
-  void post_cluster(const model& m) override {
-    t1 += tmp;
-    t2 += tmp * tmp;
-    if (n > wait) {
-      file << std::endl;
-      glFlush();
-      sleep(2);
-    }
-    n++;
-  }
-
-  void clear() {
-    t1 = 0;
-    t2 = 0;
-    n = 0;
-  }
-
-  double var() { return (t2 - t1 * t1 / (double)n) / (double)n; }
-};
-
-Gen<double, D, Euclidean<double, D>, double> eGen(double ε) {
-  return [ε](model& M, Rng& rng) -> Transformation<double, D, Euclidean<double, D>, double>* {
-    Vector<double, D> t;
-    Matrix<double, D> m;
-
-    double θ = rng.uniform((double)0.0, 2 * M_PI);
-    m(0, 0) = -cos(2 * θ);
-    m(1, 1) = cos(2 * θ);
-    m(0, 1) = -2 * cos(θ) * sin(θ);
-    m(1, 0) = -2 * cos(θ) * sin(θ);
-
-    Spin<double, D, double>* s = rng.pick(M.s);
-    t = s->x;
-    for (unsigned j = 0; j < D; j++) {
-      t(j) += rng.variate<double, std::normal_distribution>(0.0, ε);
-    }
-
-    Euclidean<double, D> g(t - m * t, m);
-    return new SpinFlip<double, D, Euclidean<double, D>, double>(M, g, s);
-  };
-}
-
-Gen<double, D, Euclidean<double, D>, double> mGen(double ε) {
-  return [ε](model& M, Rng& rng) -> Transformation<double, D, Euclidean<double, D>, double>* {
-    Matrix<double, D> m;
-
-    Spin<double, D, double>* s1 = rng.pick(M.s);
-    Spin<double, D, double>* s2 = rng.pick(M.s);
-
-    while (s1 == s2) {
-      s2 = rng.pick(M.s);
-    }
-
-    Vector<double, D> t1 = s1->x;
-    Vector<double, D> t2 = s2->x;
-    Vector<double, D> t12 = t1 - t2;
-    Vector<double, D> t = (t1 + t2) / 2;
-
-    double θ = atan2(t12(1), t12(0)) + rng.variate<double, std::normal_distribution>(0.0, ε) / t12.norm();
-
-    m(0, 0) = -cos(2 * θ);
-    m(1, 1) = cos(2 * θ);
-    m(0, 1) = -2 * cos(θ) * sin(θ);
-    m(1, 0) = -2 * cos(θ) * sin(θ);
-
-    Vector<double, D> t3 = t - m * t;
-
-    if (t3(0) != t3(0)) {
-      std::cout << t3 << "\n" << t << "\n" << m << "\n" << t12 << "\n" ;
-      getchar();
-    }
-
-    Euclidean<double, D> g(t3, m);
-    return new PairFlip<double, D, Euclidean<double, D>, double>(M, g, s1, s2);
-  };
-}
-
-Gen<double, D, Euclidean<double, D>, double> tGen(double ε) {
-  return [ε](model& M, Rng& rng) -> Transformation<double, D, Euclidean<double, D>, double>* {
-    Matrix<double, D> m;
-
-    Spin<double, D, double>* s1 = rng.pick(M.s);
-    Spin<double, D, double>* s2 = rng.pick(M.s);
-
-    while (s1 == s2) {
-      s2 = rng.pick(M.s);
-    }
-
-    Vector<double, D> t1 = s1->x;
-    Vector<double, D> t2 = s2->x;
-    Vector<double, D> t12 = t1 - t2;
-    Vector<double, D> t = t2;
-
-    double θ = atan2(t12(1), t12(0)) + rng.variate<double, std::normal_distribution>(0.0, ε) / t12.norm();
-
-    m(0, 0) = -cos(2 * θ);
-    m(1, 1) = cos(2 * θ);
-    m(0, 1) = -2 * cos(θ) * sin(θ);
-    m(1, 0) = -2 * cos(θ) * sin(θ);
-
-    Vector<double, D> t3 = t - m * t;
-
-    Euclidean<double, D> g(t3, m);
-    return new SpinFlip<double, D, Euclidean<double, D>, double>(M, g, s1);
-  };
-}
-
-Gen<double, D, Euclidean<double, D>, double> rGen(double ε) {
-  return [ε](model& M, Rng& rng) -> Transformation<double, D, Euclidean<double, D>, double>* {
-    Vector<double, D> t;
-    Matrix<double, D> m;
-
-    double θ = rng.uniform((double)0.0, 2 * M_PI);
-    m(0, 0) = -cos(2 * θ);
-    m(1, 1) = cos(2 * θ);
-    m(0, 1) = -2 * cos(θ) * sin(θ);
-    m(1, 0) = -2 * cos(θ) * sin(θ);
-
-    Spin<double, D, double>* s = rng.pick(M.s);
-    t = M.s0.t;
-    for (unsigned j = 0; j < D; j++) {
-      t(j) += rng.variate<double, std::normal_distribution>(0.0, ε);
-    }
-
-    Euclidean<double, D> g(t - m * t, m);
-    return new SpinFlip<double, D, Euclidean<double, D>, double>(M, g, s);
-  };
-}
-
 int main(int argc, char* argv[]) {
   const unsigned D = 2;
 
@@ -270,26 +56,6 @@ int main(int argc, char* argv[]) {
     }
   }
 
-  std::function<double(const Spin<double, D, double>&, const Spin<double, D, double>&)> Z =
-      [L, a, k](const Spin<double, D, double>& s1, const Spin<double, D, double>& s2) -> double {
-    Vector<double, D> d = s1.x - s2.x;
-
-    double σ = s1.s + s2.s;
-    double δ = σ - sqrt(d.transpose() * d);
-
-    if (δ > -a * σ) {
-      return 0.5 * k * (2 * pow(a * σ, 2) - pow(δ, 2));
-    } else if (δ > -2 * a * σ) {
-      return 0.5 * k * pow(δ + 2 * a * σ, 2);
-    } else {
-      return 0;
-    }
-  };
-
-  std::function<double(Spin<double, D, double>)> B = [L, H](Spin<double, D, double> s) -> double {
-    return H * s.x.norm();
-  };
-
   std::function<double(Spin<double, D, double>)> B_hard = [L, H](Spin<double, D, double> s) -> double {
     if (fabs(s.x(0)) < 3 * L / 4 && fabs(s.x(1)) < 3 * L / 4) {
       return 0;
@@ -298,10 +64,10 @@ int main(int argc, char* argv[]) {
     }
   };
 
-  auto g1 = eGen(1);
-  auto g2 = mGen(0.1);
-  auto g3 = tGen(0.1);
-  auto g4 = rGen(0);
+  auto g1 = nudgeGen<D, Radius>(1);
+  auto g2 = swapGen<D, Radius>(0.1);
+  auto g3 = accrossGen<D, Radius>(0.1);
+  auto g4 = centerGen<D, Radius>(0);
 
   auto tag = std::chrono::high_resolution_clock::now();
 
@@ -311,8 +77,8 @@ int main(int argc, char* argv[]) {
   std::ofstream file;
   file.open(filename);
 
-  animation A(L, 750, argc, argv, wait, file);
-  model sphere(1.0, Z, B);
+  Animation<double, D, Euclidean<double, D>, Radius> A(L, 750, argc, argv, wait);
+  model sphere(1.0, zSpheres<D>(a, k), bCenter<D, Radius>(H));
 
   Rng rng;
 
diff --git a/spin.hpp b/spin.hpp
index bb8369b..d4d39cc 100644
--- a/spin.hpp
+++ b/spin.hpp
@@ -7,7 +7,5 @@ template <class T, int D, class S> class Spin {
 public:
   Vector<T, D> x;
   S s;
-
-  Spin() { x = Vector<T, D>::Zero(); }
 };