From 468c25eabd1aaac7d02988fef97b66bb378988b3 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Wed, 26 Feb 2020 09:14:04 -0500
Subject: More cleanup and temperary saving routines for presentation

---
 spheres.hpp | 55 +++++++++++++++++++++++++++++++++++--------------------
 1 file changed, 35 insertions(+), 20 deletions(-)

(limited to 'spheres.hpp')

diff --git a/spheres.hpp b/spheres.hpp
index 60154f6..2ca0851 100644
--- a/spheres.hpp
+++ b/spheres.hpp
@@ -1,26 +1,26 @@
 
 #include "space_wolff.hpp"
 
-template <int 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));
-    } else if (δ > -2 * a * σ) {
-      return 0.5 * k * pow(δ + 2 * a * σ, 2);
-    } else {
-      return 0;
-    }
+template <int D> using Sphere = Spin<double, D, Radius>;
+
+template <int 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));
+  } else if (δ > -2 * a * σ) {
+    return 0.5 * k * pow(δ + 2 * a * σ, 2);
+  } 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> 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>
@@ -35,8 +35,23 @@ 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 {
-      return softPotential(a, k, diff(L, s1.x, s2.x), s1.s + s2.s);
+  return [L, a, k](const Spin<double, D, double>& s1, const Spin<double, D, double>& s2) -> double {
+    return softPotential(a, k, diff(L, s1.x, s2.x), s1.s + s2.s);
+  };
+}
+
+template <int D>
+std::function<double(const Spin<double, D, Dimer<double, D>>&,
+                     const Spin<double, D, Dimer<double, D>>&)>
+zDimers(std::function<double(const Sphere<D>&, const Sphere<D>&)> zSingle) {
+  return [zSingle](const Spin<double, D, Dimer<double, D>>& s1,
+            const Spin<double, D, Dimer<double, D>>& s2) -> double {
+    Spin<double, D, Radius> s11 = {.x = s1.x + s1.s.relativePosition, .s = s1.s.radius};
+    Spin<double, D, Radius> s12 = {.x = s1.x - s1.s.relativePosition, .s = s1.s.radius};
+    Spin<double, D, Radius> s21 = {.x = s2.x + s2.s.relativePosition, .s = s2.s.radius};
+    Spin<double, D, Radius> s22 = {.x = s2.x - s2.s.relativePosition, .s = s2.s.radius};
+
+    return zSingle(s11, s21) + zSingle(s12, s21) + zSingle(s11, s22) + zSingle(s12, s22);
   };
 }
 
-- 
cgit v1.2.3-54-g00ecf