spherocylinder.h

#pragma once

#include "potentials_base.h"
#include "units.h"

namespace Faunus::SpheroCylinder {

inline Point vec_perpproject(const Point& a, const Point& b)
{
    return a - b * a.dot(b);
}

Point mindist_segment2segment(const Point& dir1, double halfl1, const Point& dir2, double halfl2,
                              const Point& r_cm);

[[maybe_unused]] inline Point mindist_segment2point(const Point& segment_direction,
                                                    const double half_length,
                                                    const Point& separation)
{
    const auto c = segment_direction.dot(separation);
    double d;
    if (c > half_length) {
        d = half_length;
    }
    else {
        d = (c > -half_length) ? c : -half_length;
    }
    return -separation + (segment_direction * d);
}

int find_intersect_plane(const Cigar& part1, const Cigar& part2, const Point& r_cm,
                         const Point& w_vec, double cutoff_squared, double cospatch,
                         std::array<double, 5>& intersections);

int test_intrpatch(const Cigar& part1, Point& vec, double cospatch, double ti,
                   std::array<double, 5>& intersections);

int find_intersect_planec(const Cigar& part1, const Cigar& part2, const Point& r_cm,
                          const Point& w_vec, double rcut2, double cospatch,
                          std::array<double, 5>& intersections);

int psc_intersect(const Cigar& particle1, const Cigar& particle2, const Point& r_cm,
                  std::array<double, 5>& intersections, double cutoff_squared);

int cpsc_intersect(const Cigar& part1, const Cigar& part2, const Point& r_cm,
                   std::array<double, 5>& intersections, double rcut2);

inline double fanglscale(const double a, const Cigar& cigar)
{
    // a = r_ij * n_i
    if (a <= cigar.pcanglsw) {
        return 0.0;
    }
    if (a >= cigar.pcangl) {
        return 1.0;
    }
    return 0.5 - ((cigar.pcanglsw + cigar.pcangl) * 0.5 - a) / (cigar.pcangl - cigar.pcanglsw);
}

} // namespace Faunus::SpheroCylinder

// --------------------

namespace Faunus::pairpotential {

class HardSpheroCylinder : public PairPotential
{
  public:
    double operator()(const Particle& particle1, const Particle& particle2,
                      [[maybe_unused]] double d,
                      const Point& center_to_center_distance) const override
    {
        assert(particle1.hasExtension() && particle2.hasExtension());
        auto minimum_distance_squared =
            SpheroCylinder::mindist_segment2segment(
                particle1.ext->scdir, particle1.ext->half_length, particle2.ext->scdir,
                particle1.ext->half_length, center_to_center_distance)
                .squaredNorm();
        const auto contact_distance = 0.5 * (particle1.traits().sigma + particle2.traits().sigma);
        if (minimum_distance_squared < contact_distance * contact_distance) {
            return pc::infty;
        }
        return 0.0;
    } 

    HardSpheroCylinder();
    void to_json(json& j) const override;
    void from_json(const json& j) override;
};

template <pairpotential::RequirePairPotential PatchPotential,
          pairpotential::RequirePairPotential CylinderPotential>
class CigarWithSphere : public PairPotential
{
  private:
    PatchPotential patch_potential;       
    CylinderPotential cylinder_potential; 

  public:
    CigarWithSphere();
    void to_json(json& j) const override;
    void from_json(const json& j) override;

    inline double operator()(const Particle& cigar, const Particle& sphere,
                             [[maybe_unused]] double distance_squared,
                             const Point& center_separation) const override
    {
        assert(cigar.hasExtension());
        const auto c = cigar.getExt().scdir.dot(center_separation);
        double contt;
        if (c > cigar.ext->half_length) {
            contt = cigar.ext->half_length;
        }
        else {
            if (c > -cigar.ext->half_length) {
                contt = c;
            }
            else {
                contt = -cigar.ext->half_length;
            }
        }
        const Point distvec = -center_separation + (cigar.ext->scdir * contt);
        if (cigar.traits().sphero_cylinder.type == SpheroCylinderData::PatchType::None &&
            sphere.traits().sphero_cylinder.type == SpheroCylinderData::PatchType::None) {
            return cylinder_potential(cigar, sphere, distvec.squaredNorm(), Point::Zero());
        }

        // patchy interaction
        const auto cutoff_squared = patch_potential.cutOffSquared(cigar.id, sphere.id);
        // scaling function: angular dependence of patch1
        const Point vec1 = SpheroCylinder::vec_perpproject(distvec, cigar.ext->scdir).normalized();
        const auto s = vec1.dot(cigar.ext->patchdir);
        const auto f1 = SpheroCylinder::fanglscale(s, cigar.getExt());

        // scaling function for the length of spherocylinder within cutoff
        const auto ndist_squared = distvec.dot(distvec);
        const auto t = sqrt(cutoff_squared - ndist_squared); // TODO cutoff
        double f0;
        if (contt + t > cigar.ext->half_length) {
            f0 = cigar.ext->half_length;
        }
        else {
            f0 = contt + t;
        }
        if (contt - t < -cigar.ext->half_length) {
            f0 -= -cigar.ext->half_length;
        }
        else {
            f0 -= contt - t;
        }
        return f1 * (f0 + 1.0) * patch_potential(cigar, sphere, ndist_squared, Point::Zero()) +
               cylinder_potential(cigar, sphere, ndist_squared, Point::Zero());
    }
};

template <pairpotential::RequirePairPotential PatchPotential,
          pairpotential::RequirePairPotential CylinderPotential>
class CigarWithCigar : public PairPotential
{
  private:
    PatchPotential patch_potential;       
    CylinderPotential cylinder_potential; 

    [[nodiscard]] double patchyPatchyEnergy(const Particle& particle1, const Particle& particle2,
                                            const Point& center_separation) const;

    [[nodiscard]] double isotropicIsotropicEnergy(const Particle& particle1,
                                                  const Particle& particle2,
                                                  const Point& center_separation) const;

  public:
    CigarWithCigar();
    void to_json(json& j) const override;
    void from_json(const json& j) override;

    inline double operator()(const Particle& particle1, const Particle& particle2,
                             [[maybe_unused]] double center_separation_squared,
                             const Point& center_separation) const override
    {
        assert(particle1.hasExtension() && particle2.hasExtension());
        if (particle1.traits().sphero_cylinder.type != SpheroCylinderData::PatchType::None &&
            particle2.traits().sphero_cylinder.type != SpheroCylinderData::PatchType::None) {
            return patchyPatchyEnergy(particle1, particle2, center_separation);
        }
        return isotropicIsotropicEnergy(particle1, particle2, center_separation);
    }
};

template <pairpotential::RequirePairPotential PatchPotential,
          pairpotential::RequirePairPotential CylinderPotential,
          pairpotential::RequirePairPotential SphereWithSphere = CylinderPotential>
class CompleteCigarPotential : public PairPotential
{
  private:
    SphereWithSphere sphere_sphere;                                // pair potential between spheres
    CigarWithCigar<PatchPotential, CylinderPotential> cigar_cigar; // pair potential between cigars
    CigarWithSphere<PatchPotential, CylinderPotential>
        cigar_sphere; // pair potential cigar <-> sphere

  public:
    inline double operator()(const Particle& a, const Particle& b,
                             [[maybe_unused]] double distance_squared,
                             const Point& distance) const override
    {
        const double small_number = 1e-6;
        const auto a_is_sphere = !a.hasExtension() || a.ext->half_length < small_number;
        const auto b_is_sphere = !a.hasExtension() || b.ext->half_length < small_number;

        if (!a_is_sphere && !b_is_sphere) {
            return cigar_cigar(a, b, distance_squared, distance);
        }
        if (a_is_sphere && b_is_sphere) {
            return sphere_sphere(a, b, distance_squared, distance);
        }
        if (b_is_sphere) {
            return cigar_sphere(a, b, distance_squared, distance);
        }
        return cigar_sphere(b, a, distance_squared, distance);
    }

    void to_json([[maybe_unused]] json& j) const override;
    void from_json(const json& j) override;
    CompleteCigarPotential();
};

} // namespace Faunus::pairpotential