From 1586d6b2247c510c11d9a4847bde67cb0947243d Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Wed, 14 Aug 2019 11:45:12 -0400
Subject: fixed bugs in Ising, p for bonds pre-computed

---
 ising.cpp | 110 +++++++++++++++++++++++++++++++++++++++++++++-----------------
 1 file changed, 80 insertions(+), 30 deletions(-)

(limited to 'ising.cpp')

diff --git a/ising.cpp b/ising.cpp
index ca2e2e7..fffce19 100644
--- a/ising.cpp
+++ b/ising.cpp
@@ -3,7 +3,7 @@
 
 std::function<double(spin<signed, 2, signed>)> B_sin(unsigned L, unsigned n, double H) {
   return [n, H, L] (spin<signed, 2, signed> s) -> double {
-    return H * cos(2 * M_PI * n * s.x[0] / ((double)L));
+    return H * s.s * cos(2 * M_PI * n * s.x(0) / ((double)L));
   };
 }
 
@@ -14,13 +14,14 @@ int main(int argc, char* argv[]) {
   unsigned N = 1000;
   unsigned mod = 0;
   unsigned multi = 1e4;
+  double mag = 0.5;
   double T = 2.0 / log(1.0 + sqrt(2.0));
   double H = 1.0;
   double ε = 0.1;
 
   int opt;
 
-  while ((opt = getopt(argc, argv, "N:L:T:H:e:m:M")) != -1) {
+  while ((opt = getopt(argc, argv, "N:L:T:H:e:m:M:n:")) != -1) {
     switch (opt) {
       case 'N': 
         N = (unsigned)atof(optarg);
@@ -43,42 +44,76 @@ int main(int argc, char* argv[]) {
       case 'M':
         multi = atoi(optarg);
         break;
+      case 'n':
+        mag = atof(optarg);
+        break;
       default:
         exit(1);
     }
   }
 
-  std::function<double(spin<signed, D, signed>, spin<signed, D, signed>)> Z =
-    [] (spin<signed, D, signed> s1, spin<signed, D, signed> s2) -> double {
-      bool one_one = false;
-      bool many_ones = false;
-      bool any_two = false;
+  double pZ = 1.0 - exp(-1.0 / T);
+
+  std::function<double(spin<signed, D, signed>, spin<signed, D, signed>, spin<signed, D, signed>)> Z =
+    [L, pZ] (spin<signed, D, signed> s1, spin<signed, D, signed> s2, spin<signed, D, signed> s1_new) -> double {
+      bool old_one_one = false;
+      bool old_many_ones = false;
+      bool old_any_two = false;
+      bool new_one_one = false;
+      bool new_many_ones = false;
+      bool new_any_two = false;
+
+      vector<signed, D> old_diff = diff<signed, D>(L, s1.x, s2.x);
+      vector<signed, D> new_diff = diff<signed, D>(L, s1_new.x, s2.x);
       
       for (unsigned i = 0; i < D; i++) {
-        unsigned diff = abs(s1.x(i) - s2.x(i));
-        if (diff == 1 && !one_one) {
-          one_one = true;
-        } else if (diff == 1 && one_one) {
-          many_ones = true;
-          break;
-        } else if (diff > 1) {
-          any_two = true;
-          break;
+        if (old_diff(i) == 1 && !old_one_one) {
+          old_one_one = true;
+        } else if (old_diff(i) == 1 && old_one_one) {
+          old_many_ones = true;
+        } else if (old_diff(i) > 1) {
+          old_any_two = true;
+        }
+        if (new_diff(i) == 1 && !new_one_one) {
+          new_one_one = true;
+        } else if (new_diff(i) == 1 && new_one_one) {
+          new_many_ones = true;
+        } else if (new_diff(i) > 1) {
+          new_any_two = true;
         }
       }
 
-      if (!one_one && !any_two) {
-        return -std::numeric_limits<double>::infinity();
-      } else if (one_one && !many_ones && !any_two) {
-        return s1.s * s2.s;
+      bool were_on_someone = !old_one_one && !old_any_two;
+      bool are_on_someone = !new_one_one && !new_any_two;
+      bool were_nearest_neighbors = old_one_one && !old_many_ones && !old_any_two;
+      bool are_nearest_neighbors = new_one_one && !new_many_ones && !new_any_two;
+
+      if (were_on_someone) {
+        return 0.0;
+      } else if (are_on_someone) {
+        return 1.0;
+      } else if (were_nearest_neighbors && are_nearest_neighbors) {
+        return 0.0;
+      } else if (were_nearest_neighbors) {
+        if (s1.s * s2.s == 1) {
+          return pZ;
+        } else {
+          return 0.0;
+        }
+      } else if (are_nearest_neighbors) {
+        if (s1_new.s * s2.s == -1) {
+          return pZ;
+        } else {
+          return 0.0;
+        }
       } else {
-        return 0;
+        return 0.0;
       }
     };
 
   std::function<double(spin<signed, D, signed>)> B_face =
     [L, H] (spin<signed, D, signed> s) -> double {
-      return H * s.s * smiley[s.x(1) * 16 / L][s.x(0) * 16 / L];
+      return H * s.s * smiley[s.x(0) * 16 / L][s.x(1) * 16 / L];
     };
 
   std::function<double(spin<signed, D, signed>)> B;
@@ -118,7 +153,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 ((coin(rng) && up < pow(L, 2) / 2) || down >= pow(L, 2) / 2) {
+      if ((coin(rng) && up < pow(L, 2) * mag) || down >= pow(L, 2) * mag) {
         ising.s.push_back({{i, j}, 1});
         up++;
       } else {
@@ -154,13 +189,6 @@ int main(int argc, char* argv[]) {
     }
   }
 
-  std::vector<signed> output(pow(L, D));
-
-  for (spin<signed, D, signed> s : ising.s) {
-    spin<signed, D, signed> rs = ising.s0.inverse().act(s);
-    output[L * rs.x(1) + rs.x(0)] = s.s;
-  }
-
   std::ofstream outfile;
   outfile.open("out.dat", std::ios::app);
 
@@ -173,6 +201,28 @@ int main(int argc, char* argv[]) {
   }
   outfile << "\n";
 
+  std::ofstream snapfile;
+  snapfile.open("snap.dat");
+
+  std::vector<std::vector<signed>> snap(L);
+  for (std::vector<signed>& line : snap) {
+    line.resize(L);
+  }
+
+  for (spin<signed, D, signed> s : ising.s) {
+    spin<signed, D, signed> snew = ising.s0.inverse().act(s);
+    snap[snew.x(0)][snew.x(1)] = snew.s;
+  }
+
+  for (std::vector<signed> row : snap) {
+    for (signed s : row) {
+      snapfile << s << " ";
+    }
+    snapfile << "\n";
+  }
+
+  snapfile.close();
+
   return 0;
 }
 
-- 
cgit v1.2.3-54-g00ecf