3 files changed, 44 insertions, 19 deletions

diff --git a/ising.cpp b/ising.cpp
index 38dce3b..838ec9e 100644
--- a/ising.cpp
+++ b/ising.cpp
@@ -13,9 +13,11 @@ private:
   uint64_t t2;
   unsigned n;
   unsigned tmp;
+  bool color;
 
 public:
-  animation(double L, unsigned w, int argc, char* argv[]) {
+  animation(double L, unsigned w, bool tcolor, int argc, char* argv[]) {
+    color = tcolor;
     t1 = 0;
     t2 = 0;
     n = 0;
@@ -31,14 +33,19 @@ public:
   }
 
   void post_cluster(const model& m) override {
-    glClearColor(1.0f, 1.0f, 1.0f, 1.0f);
+    glClearColor(1.0, 1.0, 1.0, 1.0);
     glClear(GL_COLOR_BUFFER_BIT);
-    glClearColor(0.0, 1.0, 0.0, 0.0);
     for (const Spin<signed, 2, signed>* s : m.s) {
       if (s->s == 1) {
-        glColor3f(0.0, 0.0, 0.0);
+        if (color)
+          glColor3f(1.0, 0.0, 0.0);
+        else
+          glColor3f(1.0, 1.0, 1.0);
       } else if (s->s == -1) {
-        glColor3f(1.0, 1.0, 1.0);
+        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);
@@ -54,13 +61,15 @@ int main(int argc, char* argv[]) {
   unsigned mod = 0;
   unsigned multi = 1e4;
   double mag = 0.5;
+  double pop = 1.0;
   double T = 2.0 / log(1.0 + sqrt(2.0));
   double H = 1.0;
   double ε = 0.1;
+  bool color = false;
 
   int opt;
 
-  while ((opt = getopt(argc, argv, "N:L:T:H:e:m:M:n:")) != -1) {
+  while ((opt = getopt(argc, argv, "N:L:T:H:e:m:M:r:p:c")) != -1) {
     switch (opt) {
     case 'N':
       N = (unsigned)atof(optarg);
@@ -83,9 +92,15 @@ int main(int argc, char* argv[]) {
     case 'M':
       multi = atoi(optarg);
       break;
-    case 'n':
+    case 'r':
       mag = atof(optarg);
       break;
+    case 'p':
+      pop = atof(optarg);
+      break;
+    case 'c':
+      color = true;
+      break;
     default:
       exit(1);
     }
@@ -124,12 +139,18 @@ int main(int argc, char* argv[]) {
     }
   };
 
-  std::function<double(const Spin<signed, D, signed>&)> B_face =
-      [L, H](const Spin<signed, D, signed>& s) -> double {
-    return H * s.s * smiley[15 - s.x(1) * 16 / L][s.x(0) * 16 / L];
-  };
+  std::function<double(Spin<signed, D, signed>)> B;
+
+  if (mod == 0) {
+    B = [L, H](const Spin<signed, D, signed>& s) -> double {
+      return H * s.s * smiley[15 - s.x(1) * 16 / L][s.x(0) * 16 / L];
+    };
+  } else {
+    B = [mod, L, H](const Spin<signed, D, signed>& s) -> double {
+      return H * s.s * sin(s.x(0) * mod * 2 * M_PI / L);
+    };
+  }
 
-  std::function<double(Spin<signed, D, signed>)> B = B_face;
 
   Model<signed, D, TorusGroup<signed, D>, signed> ising(L, Z, B);
 
@@ -140,6 +161,7 @@ int main(int argc, char* argv[]) {
   unsigned down = 0;
   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) {
@@ -152,6 +174,7 @@ int main(int argc, char* argv[]) {
       ising.s.push_back(ss);
       ising.dict.insert(ss);
       n++;
+      }
     }
   }
 
@@ -182,10 +205,14 @@ int main(int argc, char* argv[]) {
 
     std::set<Spin<signed, 2, signed>*> s2s = M.dict.at(s2.x);
 
-    return new PairFlip<signed, 2, TorusGroup<signed, 2>, signed>(M, g, ss, *s2s.begin());
+    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());
+    }
   };
 
-  animation A(L, 750, argc, argv);
+  animation A(L, 750, color, argc, argv);
 
   ising.wolff(T, {g}, A, N);
 
diff --git a/space_wolff.hpp b/space_wolff.hpp
index c4a7a21..9d529a3 100644
--- a/space_wolff.hpp
+++ b/space_wolff.hpp
@@ -31,8 +31,6 @@ public:
   virtual void post_cluster(const Model<U, D, R, S>&){};
 };
 
-
-
 template <class U, int D, class R, class S> class Model {
 public:
   U L;
diff --git a/spheres_infinite.cpp b/spheres_infinite.cpp
index 026486e..5287a84 100644
--- a/spheres_infinite.cpp
+++ b/spheres_infinite.cpp
@@ -179,8 +179,8 @@ int main(int argc, char* argv[]) {
     return H * s.x.norm();
   };
 
-  auto g1 = eGen(0.5);
-  auto g2 = mGen(0.2);
+  auto g1 = eGen(1);
+  auto g2 = mGen(0.5);
   animation A(L, 750, argc, argv);
   model sphere(1.0, Z, B);
 
@@ -192,7 +192,7 @@ int main(int argc, char* argv[]) {
   for (unsigned i = 0; i < sphere.s.size(); i++) {
     Spin<double, 2, double>* ss = new Spin<double, 2, double>();
     ss->x = {(i / nx) * L / nx, (i % nx) * L / nx};
-    ss->s = rng.uniform<double>(0.25, 0.45);
+    ss->s = rng.uniform<double>(0.45, 0.45);
     sphere.s[i] = ss;
     sphere.dict.insert(ss);
   }