diff options
author | Jaron Kent-Dobias <jaron@kent-dobias.com> | 2020-02-26 00:00:16 -0500 |
---|---|---|
committer | Jaron Kent-Dobias <jaron@kent-dobias.com> | 2020-02-26 00:00:16 -0500 |
commit | 6461165f4daea01c5baa1d969ee5c726fb4c559e (patch) | |
tree | 9946c3a04f5dd62e0cc78fcd34089613d0dd0c48 /spheres.hpp | |
parent | 98350b34d1826e1d6687f18c45b0fbc6a3488742 (diff) | |
download | space_wolff-6461165f4daea01c5baa1d969ee5c726fb4c559e.tar.gz space_wolff-6461165f4daea01c5baa1d969ee5c726fb4c559e.tar.bz2 space_wolff-6461165f4daea01c5baa1d969ee5c726fb4c559e.zip |
Finished fixes of new features.
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); }; } |