began transition to more general transformations

1 files changed, 34 insertions, 6 deletions

diff --git a/spheres_infinite.cpp b/spheres_infinite.cpp
index 9bc86d7..a9cd7d0 100644
--- a/spheres_infinite.cpp
+++ b/spheres_infinite.cpp
@@ -67,7 +67,7 @@ public:
   double var() { return (t2 - t1 * t1 / (double)n) / (double)n; }
 };
 
-std::function<Euclidean<double, D>(const model&, randutils::mt19937_rng&)> eGen(double ε) {
+Gen<double, D, Euclidean<double, D>, double> eGen(double ε) {
   return [ε](const model& M, randutils::mt19937_rng& rng) -> Euclidean<double, 2> {
     Vector<double, D> t;
     Matrix<double, D> m;
@@ -89,6 +89,31 @@ std::function<Euclidean<double, D>(const model&, randutils::mt19937_rng&)> eGen(
   };
 }
 
+Gen<double, D, Euclidean<double, D>, double> mGen(double ε) {
+  return [ε](const model& M, randutils::mt19937_rng& rng) -> Euclidean<double, 2> {
+    Matrix<double, D> m;
+
+    unsigned f_ind1 = rng.uniform((unsigned)0, (unsigned)M.s.size());
+    unsigned f_ind2 = rng.uniform((unsigned)0, (unsigned)M.s.size() - 1);
+    if (f_ind2 >= f_ind1) f_ind2++;
+
+    Vector<double, D> t1 = M.s[f_ind1].x;
+    Vector<double, D> t2 = M.s[f_ind2].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, ε);
+
+    m(0, 0) = -cos(2 * θ);
+    m(1, 1) = cos(2 * θ);
+    m(0, 1) = -2 * cos(θ) * sin(θ);
+    m(1, 0) = -2 * cos(θ) * sin(θ);
+
+    Euclidean<double, D> g(t - m * t, m);
+    return g;
+  };
+}
+
 int main(int argc, char* argv[]) {
   const unsigned D = 2;
 
@@ -122,8 +147,8 @@ int main(int argc, char* argv[]) {
     }
   }
 
-  double k = 1e8;
-  double a = 0.0;
+  double k = 1000;
+  double a = 0.05;
 
   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 {
@@ -145,7 +170,8 @@ int main(int argc, char* argv[]) {
     return H * s.x.norm();
   };
 
-  auto g = eGen(0.05);
+  auto g1 = eGen(0.25);
+  auto g2 = mGen(0.1);
   animation A(L, 750, argc, argv);
   model sphere(1, Z, B);
 
@@ -156,14 +182,15 @@ int main(int argc, char* argv[]) {
   unsigned nx = floor(sqrt(n));
   for (unsigned i = 0; i < n; i++) {
     Vector<double, D> pos = {(i / nx) * L / nx, (i % nx) * L / nx};
-    sphere.s.push_back({pos, 0.5});
+    sphere.s.push_back({pos, rng.uniform<double>(0.45, 0.45)});
     sphere.dict.insert(&sphere.s.back());
   }
 
-  sphere.wolff(T, g, A, N);
+  sphere.wolff(T, {g1, g2}, A, N);
   std::ofstream outfile;
   outfile.open("test.dat");
 
+  /*
   for (signed i = -10; i <= 10; i++) {
     A.clear();
     double ε = pow(2, -4 + i / 2.0);
@@ -172,6 +199,7 @@ int main(int argc, char* argv[]) {
     outfile << ε << " " << A.var() / sphere.s.size() << std::endl;
     std::cout << ε << " " << A.var() / sphere.s.size() << std::endl;
   }
+  */
 
   outfile.close();