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