5 files changed, 24 insertions, 23 deletions

diff --git a/dimers.cpp b/dimers.cpp
index 4416333..d5cfdf4 100644
--- a/dimers.cpp
+++ b/dimers.cpp
@@ -85,7 +85,7 @@ int main(int argc, char* argv[]) {
   file.open(filename);
 
   Animation<double, D, Euclidean<double, D>, Dimer<double, D>> A(L, 750, argc, argv, wait);
-  model sphere(1.0, Z, bCenter<D, Dimer<double, D>>(H));
+  model sphere(1, Z, bCenter<D, Dimer<double, D>>(H));
 
   Rng rng;
 
diff --git a/dimers_torus.cpp b/dimers_torus.cpp
index 74cf47e..ec7c73a 100644
--- a/dimers_torus.cpp
+++ b/dimers_torus.cpp
@@ -86,7 +86,7 @@ int main(int argc, char* argv[]) {
 
   sphere.s.resize(n);
 
-  unsigned nx = floor(sqrt(n));
+  unsigned nx = ceil(sqrt(n));
   for (unsigned i = 0; i < sphere.s.size(); i++) {
     Spin<double, 2, Dimer<double, D>>* ss = new Spin<double, 2, Dimer<double, D>>();
     ss->x = {(i / nx) * L / nx, (i % nx) * L / nx};
diff --git a/spheres.cpp b/spheres.cpp
index 9ca2036..4b2237d 100644
--- a/spheres.cpp
+++ b/spheres.cpp
@@ -75,7 +75,7 @@ int main(int argc, char* argv[]) {
 
   sphere.s.resize(n);
 
-  unsigned nx = floor(sqrt(n));
+  unsigned nx = ceil(sqrt(n));
   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};
diff --git a/spheres.hpp b/spheres.hpp
index 3dc4a5e..60154f6 100644
--- a/spheres.hpp
+++ b/spheres.hpp
@@ -2,13 +2,8 @@
 #include "space_wolff.hpp"
 
 template <int D>
-std::function<double(const Spin<double, D, Radius>&, const Spin<double, D, Radius>&)>
-zSpheres(double a, double k) {
-  return [a, k](const Spin<double, D, Radius>& s1, const Spin<double, D, Radius>& s2) -> double {
-    Vector<double, D> d = s1.x - s2.x;
-
-    double σ = s1.s + s2.s;
-    double δ = σ - sqrt(d.transpose() * d);
+double softPotential(double a, double k, const Vector<double, D>& diff, double σ) {
+    double δ = σ - sqrt(diff.transpose() * diff);
 
     if (δ > -a * σ) {
       return 0.5 * k * (2 * pow(a * σ, 2) - pow(δ, 2));
@@ -17,6 +12,23 @@ zSpheres(double a, double k) {
     } else {
       return 0;
     }
+}
+
+template <int D>
+double hardPotential(const Vector<double, D>& diff, double σ) {
+    if (pow(σ, 2) < diff.transpose() * diff) {
+      return 0;
+    } else {
+      return -std::numeric_limits<double>::infinity();
+    }
+}
+
+template <int D>
+std::function<double(const Spin<double, D, Radius>&, const Spin<double, D, Radius>&)>
+zSpheres(double a, double k) {
+  return [a, k](const Spin<double, D, Radius>& s1, const Spin<double, D, Radius>& s2) -> double {
+    Vector<double, D> d = s1.x - s2.x;
+    return softPotential(a, k, d, s1.s + s2.s);
   };
 }
 
@@ -24,18 +36,7 @@ template <int D>
 std::function<double(const Spin<double, D, Radius>&, const Spin<double, D, Radius>&)>
 zSpheresTorus(double L, double a, double k) {
       return [L, a, k](const Spin<double, D, double>& s1, const Spin<double, D, double>& s2) -> double {
-    Vector<double, D> d = diff(L, s1.x, s2.x);
-
-    double σ = s1.s + s2.s;
-    double δ = σ - sqrt(d.transpose() * d);
-
-    if (δ > -a * σ) {
-      return 0.5 * k * (2 * pow(a * σ, 2) - pow(δ, 2));
-    } else if (δ > -2 * a * σ) {
-      return 0.5 * k * pow(δ + 2 * a * σ, 2);
-    } else {
-      return 0;
-    }
+      return softPotential(a, k, diff(L, s1.x, s2.x), s1.s + s2.s);
   };
 }
 
diff --git a/vector.hpp b/vector.hpp
index f642816..04a1a97 100644
--- a/vector.hpp
+++ b/vector.hpp
@@ -10,7 +10,7 @@ template <class T, int D> Vector<T, D> diff(T L, Vector<T, D> v1, Vector<T, D> v
   Vector<T, D> v;
 
   for (unsigned i = 0; i < D; i++) {
-    v(i) = std::abs(mod(v1(i), L) - mod(v2(i), L));
+    v(i) = mod(std::abs(v1(i) - v2(i)), L);
     if (v(i) > L / 2) {
       v(i) = L - v(i);
     }