diff options
Diffstat (limited to 'spheres.hpp')
-rw-r--r-- | spheres.hpp | 39 |
1 files changed, 20 insertions, 19 deletions
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); }; } |