Some refactoring.

1 files changed, 30 insertions, 37 deletions

diff --git a/soft_spheres.cpp b/soft_spheres.cpp
index cdd6f31..7944aa0 100644
--- a/soft_spheres.cpp
+++ b/soft_spheres.cpp
@@ -1,55 +1,48 @@
 #include "animation.hpp"
 
+template <int D, int n> class SoftSphere : public Sphere<D> {
+private:
+  const double c0 = -pow(2, 2*n-3) * pow(5, -n) * (8 + 6*n + pow(n, 2));
+  const double c2 =  pow(2, 2+2*n) * pow(5, -n-2) * n * (4+n);
+  const double c4 = -pow(2, 5+2*n) * pow(5, -n-4) * n * (2+n);
+public:
+  using Sphere<D>::Sphere;
+
+  double U(const Position<D>& y, double ry) const {
+    double r2 = Sphere<D>::regularizedDistanceSquared(y, ry);
+    if (r2 < pow(1.25, 2)) {
+      return pow(1 / r2, n / 2) + c0 + c2 * r2 + c4 * pow(r2, 2);
+    } else {
+      return 0;
+    }
+  }
+};
+
 int main(int argc, char* argv[]) {
-  const unsigned n = 12;
+  const unsigned n = 24;
 
   unsigned N = 1200;
-  double R = 0.023;
-  double ε = 0.01;
+  double R = 0.025;
+  double ε = 0.005;
   unsigned steps = 1e6;
-  double β = 1;
+  double β = 0.5;
 
   initializeAnimation(argc, argv);
 
   Rng r;
 
-  Model<2, SoftSphere<2, n>> m(N, 2 * R * 1.25);
+  Model<2, SoftSphere<2, n>> m(N, 2 * R * 1.25, 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(SoftSphere<2, n>(x, rad));
-  }
-
-  for (unsigned i = 0; i < steps * N; i++) {
-    SoftSphere<2, n>& s = r.pick(m.particles);
-    Vector<2> Δx;
-    for (double& Δxi : Δx) {
-      Δxi = r.variate<double, std::normal_distribution>(0, ε);
+  for (unsigned i = 0; i < steps; i++) {
+    for (unsigned j = 0; j < N; j++) {
+      m.randomMove(β, ε, r);
+      m.randomSwap(β, r);
     }
 
-    for (unsigned j = 0; j < N; j++)
-      m.metropolis(β, s, Δx, r);
-
-    SoftSphere<2, n>& s1 = r.pick(m.particles);
-    SoftSphere<2, n>& 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;
+//    std::cout << m.randomClusterSwap(β, r) << std::endl;
     
-    draw(m);
+    if (i % 3 == 0)
+      draw(m);
   }
 
   return 0;