1 files changed, 64 insertions, 35 deletions
diff --git a/hard_spheres.cpp b/hard_spheres.cpp
index 6a14953..b8be1cf 100644
--- a/hard_spheres.cpp
+++ b/hard_spheres.cpp
@@ -1,53 +1,82 @@
 #include "animation.hpp"
 
+template <int D> class HardSphere : public Sphere<D> {
+public:
+  using Sphere<D>::Sphere;
+
+  double U(const Position<D>& y, double ry) const {
+    if (Sphere<D>::regularizedDistanceSquared(y, ry) < 1) {
+      return 1e6;
+    } else {
+      return 0;
+    }
+  }
+};
+
 int main(int argc, char* argv[]) {
   unsigned N = 1200;
-  double R = 0.021;
-  double ε = 0.01;
+  double R = 0.024;
+  double ε = 0.005;
   unsigned steps = 1e6;
 
   initializeAnimation(argc, argv);
 
   Rng r;
 
-  Model<2, HardSphere<2>> m(N, 2 * R);
+  Model<2, HardSphere<2>> m(N, 2 * R, gContinuousPolydisperse(R), r);
 
-  for (unsigned i = 0; i < N; i++) {
-    Position<2> x = {r.uniform(0.0, 1.0), r.uniform(0.0, 1.0)};
-    double rad = pow(r.uniform(pow(2.219, -2), 1.0), -0.5) * R / 2.219;
-    m.insert(HardSphere<2>(x, rad));
-  }
+  for (unsigned i = 0; i < steps; i++) {
+    for (unsigned j = 0; j < N; j++) {
+      HardSphere<2>& p = r.pick(m.particles);
+      Position<2> xNew = p.x;
+      Vector<2> Δx;
+      for (double& xi : Δx) {
+        xi = r.variate<double, std::normal_distribution>(0, ε);
+      }
+      xNew += Δx;
+      bool reject = false;
+      for (const HardSphere<2>* neighbor : m.nl.neighborsOf(xNew)) {
+        if (neighbor != &p) {
+          if (neighbor->regularizedDistanceSquared(xNew, p.r) < 1 && neighbor->regularizedDistanceSquared(p.x, p.r) >= 1) {
+            reject = true;
+            break;
+          }
+        }
+      }
+      if (!reject) {
+        m.move(p, xNew);
+      }
+      
+      HardSphere<2>& p1 = r.pick(m.particles);
+      HardSphere<2>& p2 = r.pick(m.particles);
 
-  for (unsigned i = 0; i < steps * N; i++) {
-    HardSphere<2>& s = r.pick(m.particles);
-    Vector<2> Δx;
-    for (double& Δxi : Δx) {
-      Δxi = r.variate<double, std::normal_distribution>(0, ε);
+      double rat = p1.r / p2.r;
+      if (rat < 1.2 && rat > 0.83) {
+        reject = false;
+        for (const HardSphere<2>* neighbor : m.nl.neighborsOf(p2.x)) {
+          if (neighbor != &p1 && neighbor != &p2) {
+            if (neighbor->regularizedDistanceSquared(p2.x, p1.r) < 1 && neighbor->regularizedDistanceSquared(p1.x, p1.r) >= 1) {
+              reject = true;
+              break;
+            }
+          }
+        }
+        for (const HardSphere<2>* neighbor : m.nl.neighborsOf(p1.x)) {
+          if (neighbor != &p1 && neighbor != &p2) {
+            if (neighbor->regularizedDistanceSquared(p1.x, p2.r) < 1 && neighbor->regularizedDistanceSquared(p2.x, p2.r) >= 1) {
+              reject = true;
+              break;
+            }
+          }
+        }
+        if (!reject) {
+          std::swap(p1.r, p2.r);
+        }
+      }
     }
 
-    m.metropolis(1, s, Δx, r);
-
-    HardSphere<2>& s1 = r.pick(m.particles);
-    HardSphere<2>& s2 = r.pick(m.particles);
-
-    Vector<2> t1 = s1.x;
-    Vector<2> t2 = s2.x;
-    Vector<2> t = (t1 + t2) / 2;
-
-    Matrix<2> mat;
-
-    mat(0, 0) = -1;
-    mat(1, 1) = -1;
-    mat(0, 1) = 0;
-    mat(1, 0) = 0;
-
-    Euclidean<2> g(t - mat * t, mat);
-
-    std::cout << m.clusterFlip(1, g, s1, r) << std::endl;
-
-    if (i % (N) == 0) {
+    if (i % 10 == 0)
       draw(m);
-    }
   }
 
   return 0;