2 files changed, 53 insertions, 9 deletions

diff --git a/euclidean.hpp b/euclidean.hpp
index b1a5e80..49ec2d9 100644
--- a/euclidean.hpp
+++ b/euclidean.hpp
@@ -25,10 +25,14 @@ 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 = t + r * s.x;
+    s_new.x = this->act(s.x);
     s_new.s = s.s;
 
     return s_new;
diff --git a/spheres_infinite.cpp b/spheres_infinite.cpp
index e41ea17..7879f34 100644
--- a/spheres_infinite.cpp
+++ b/spheres_infinite.cpp
@@ -15,10 +15,11 @@ private:
   unsigned n;
   unsigned tmp;
   unsigned wait;
+  double L;
   Euclidean<double, D> s0_tmp;
 
 public:
-  animation(double L, unsigned w, int argc, char* argv[]) : s0_tmp(0), wait(1000) {
+  animation(double L, unsigned w, int argc, char* argv[]) : s0_tmp(0), wait(1000), L(L) {
     t1 = 0;
     t2 = 0;
     n = 0;
@@ -60,11 +61,6 @@ public:
     tmp = 0;
     glClearColor(1.0f, 1.0f, 1.0f, 1.0f);
     glClear(GL_COLOR_BUFFER_BIT);
-    /*
-    glBegin(GL_LINE);
-    Vector<double, 2> r_center = s0_tmp.act({t->r.t, 0});
-    glEnd();
-    */
     for (const Spin<double, 2, double>* s : m.s) {
       glBegin(GL_POLYGON);
       unsigned n_points = 50;
@@ -75,6 +71,19 @@ public:
       }
       glEnd();
     }
+    glColor3d(1.0, 0.0, 0.0);
+    glLineWidth(3);
+    glBegin(GL_LINES);
+      Euclidean<double, D> a_tmp = t->r;
+      a_tmp.t /= 2;
+      Euclidean<double, D> r_tmp = t->r;
+      Vector<double, 2> r_center = r_tmp.t / 2;
+      double θ = acos(r_tmp.r(1, 1)) / 2;
+      Vector<double, 2> v1 = s0_tmp.act({r_center(0) + 50*L * sin(θ), r_center(1) - 50*L * cos(θ)});
+      Vector<double, 2> v2 = s0_tmp.act({r_center(0) - 50*L * sin(θ), r_center(1) + 50*L * cos(θ)});
+      glVertex2d(v1(0), v1(1));
+      glVertex2d(v2(0), v2(1));
+    glEnd();
     glFlush();
   }
 
@@ -154,6 +163,36 @@ Gen<double, D, Euclidean<double, D>, double> mGen(double ε) {
   };
 }
 
+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;
@@ -240,7 +279,8 @@ int main(int argc, char* argv[]) {
 
   auto g1 = eGen(1);
   auto g2 = mGen(0.1);
-  auto g3 = rGen(1);
+  auto g3 = tGen(0.1);
+  auto g4 = rGen(1);
   animation A(L, 750, argc, argv);
   model sphere(1.0, Z, B);
 
@@ -257,7 +297,7 @@ int main(int argc, char* argv[]) {
     sphere.dict.insert(ss);
   }
 
-  sphere.wolff(T, {g1, g2}, A, N);
+  sphere.wolff(T, {g1, g2, g3, g4}, A, N);
   std::ofstream outfile;
   outfile.open("test.dat");