1 files changed, 19 insertions, 0 deletions
diff --git a/spheres.hpp b/spheres.hpp
index 040313c..3dc4a5e 100644
--- a/spheres.hpp
+++ b/spheres.hpp
@@ -20,6 +20,25 @@ zSpheres(double a, double k) {
   };
 }
 
+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;
+    }
+  };
+}
+
 template <int D, class S> std::function<double(Spin<double, D, S>)> bCenter(double H) {
   return [H](Spin<double, D, S> s) -> double { return H * s.x.norm(); };
 }