1 files changed, 41 insertions, 14 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);