From 6461165f4daea01c5baa1d969ee5c726fb4c559e Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Wed, 26 Feb 2020 00:00:16 -0500
Subject: Finished fixes of new features.

---
 spheres.hpp | 39 ++++++++++++++++++++-------------------
 1 file changed, 20 insertions(+), 19 deletions(-)

(limited to 'spheres.hpp')

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);
   };
 }
 
-- 
cgit v1.2.3-54-g00ecf