//  ************************************************************************************************
//
//  libformfactor: efficient and accurate computation of scattering form factors
//
//! @file      ff/Polyhedron.cpp
//! @brief     Implements class Polyhedron.
//!
//! @homepage  https://jugit.fz-juelich.de/mlz/libformfactor
//! @license   GNU General Public License v3 or higher (see LICENSE)
//! @copyright Forschungszentrum Jülich GmbH 2022
//! @author    Joachim Wuttke, Scientific Computing Group at MLZ (see CITATION)
//
//  ************************************************************************************************

//! The mathematics implemented here is described in full detail in a paper
//! by Joachim Wuttke, entitled
//! "Numerically stable form factor of any polygon and polyhedron"

#include "ff/Polyhedron.h"
#include "ff/Edge.h"
#include "ff/EulerOperations.h"
#include "ff/Face.h"
#include "ff/Topology.h"
#include <algorithm>
#include <fstream>
#include <stdexcept>
#define _USE_MATH_DEFINES
#include <math.h>
#ifdef ALGORITHM_DIAGNOSTIC_LEVEL2
#include <boost/format.hpp>
#endif

namespace {

const double eps = 2e-16;
const double q_limit_series = 1e-2;
const int n_limit_series = 20;

std::vector<R3> shift_vertices(const std::vector<R3>& vertices, const R3& location)
{
    std::vector<R3> result;
    result.reserve(vertices.size());
    for (const R3 v : vertices)
        result.emplace_back(v + location);
    return result;
}

std::vector<ff::Face> precompute_faces(const ff::Topology& topology,
                                       const std::vector<R3>& vertices)
{
    if (topology.faces.size() < 4)
        throw std::runtime_error("Invalid polyhedron: less than four non-vanishing faces");
    if (topology.symmetry_Ci) {
        std::vector<ff::Face> faces;
        size_t N = topology.faces.size();
        if (N & 1)
            throw std::runtime_error("Invalid polyhedron: odd #faces violates symmetry Ci");
        for (size_t k = 0; k < N / 2; ++k) {
            const ff::FacialTopology& tf = topology.faces[k];
            std::vector<R3> corners; // of one face
            corners.reserve(N / 2);
            for (int i : topology.faces[k].vertexIndices)
                corners.push_back(vertices[i]);
            faces.emplace_back(corners, tf.symmetry_S2);
            // check against counterface
            const ff::FacialTopology& tfc = topology.faces[N - 1 - k];
            corners.clear();
            corners.reserve(N / 2);
            for (int i : tfc.vertexIndices)
                corners.push_back(vertices[i]);
            ff::Face counterface(corners, tf.symmetry_S2);
            faces[k].assert_Ci(counterface);
        }
        return faces;
    }
    std::vector<ff::Face> faces;
    for (const ff::FacialTopology& tf : topology.faces) {
        std::vector<R3> corners; // of one face
        for (int i : tf.vertexIndices)
            corners.push_back(vertices[i]);
        faces.emplace_back(corners, tf.symmetry_S2);
    }
    return faces;
}

double precompute_radius(const std::vector<R3>& vertices)
{
    return std::max_element(vertices.begin(), vertices.end(),
                            [](const R3& a, const R3& b) { return a.mag() < b.mag(); })
        ->mag();
}

double precompute_volume(const std::vector<ff::Face>& faces, bool symmetry_Ci)
{
    double volume = 0;
    for (const ff::Face& Gk : faces)
        volume += Gk.pyramidalVolume();
    if (symmetry_Ci)
        volume = 2 * volume;
    return volume;
}

//! Returns true if a ray defined by r0 -> r1 passes through a plane defined by its
//! origin p0 and normal vector n.

bool ray_intersects_plane(R3 r0, R3 r1, R3 p0, R3 n)
{
    return (r1 - r0).dot(n) != 0;
}

//! Returns the intersection point of a ray with a plane.

R3 ray_plane_intersection(R3 r0, R3 r1, R3 p0, R3 n)
{
    if ((r1 - r0).dot(n) == 0)
        throw std::runtime_error("Invalid call to libformfactor, function ray_plane_intersection: "
                                 "ray is parallel to plane");
    double d = (p0 - r0).dot(n) / ((r1 - r0).dot(n)); // projection onto ray
    return r0 + d * (r1 - r0);
}

//! Returns weather i lies in the positive eay extension defined by r0 -> r1

bool intersects_in_positive_halfspace(const R3& r0, const R3& r1, const R3& i)
{
    return (r1 - r0).dot(i - r0) > 0;
}

//! Returns the number of faces crossed by the ray r0 -> r1. Set sym to true to mirror the faces

int ray_crossings(const R3& r0, const R3& r1, const std::vector<ff::Face>& faces, bool sym)
{
    int crossings = 0;
    int sign = sym ? -1 : 1;
    for (const ff::Face& polygon : faces) {
        const R3 origin = sym ? -polygon.edges()[0].R() : polygon.edges()[0].R();
        if (!ray_intersects_plane(r0, r1, origin, (sign * polygon.normal())))
            continue;
        const R3 i = ray_plane_intersection(r0, r1, origin, (sign * polygon.normal()));
        if (!intersects_in_positive_halfspace(r0, r1, i))
            continue;
        if (polygon.is_inside(i))
            crossings++;
    }
    return crossings;
}

//! Uses the fibonacci sphere algorithm to evenly generate points on as sphere
//! See: https://stackoverflow.com/questions/9600801/evenly-distributing-n-points-on-a-sphere
//! Or for more details: González, https://doi.org/10.1007/s11004-009-9257-x (2010)
std::vector<R3> fibonacci_sphere(const R3& origin, int nsamples)
{
    static const double phi = M_PI * (3.0 - sqrt(5.0));

    std::vector<R3> result(nsamples);
    for (int i = 0; i < nsamples; i++) {
        const double y = 1 - 1 * ((double)i / (nsamples - 1)) * 2;
        const double radius = sqrt(1 - y * y);
        const double theta = phi * i;
        const double x = cos(theta) * radius;
        const double z = sin(theta) * radius;
        result[i] = R3(x, y, z) + origin;
    }
    return result;
}

} // namespace


ff::Polyhedron::Polyhedron(const Topology& topology, const std::vector<R3>& vertices,
                           const R3& location)
    : IBody(location)
    , m_topology(std::make_unique<const Topology>(topology))
    , m_vertices(vertices)
    , m_absolute_vertices(shift_vertices(vertices, location))
    , m_faces(precompute_faces(topology, vertices))
    , m_radius(precompute_radius(vertices))
    , m_volume(precompute_volume(m_faces, m_topology->symmetry_Ci))
{
}

ff::Polyhedron::~Polyhedron() = default;

const ff::Topology& ff::Polyhedron::topology() const
{
    return *m_topology;
}

const std::vector<ff::Face>& ff::Polyhedron::faces() const
{
    return m_faces;
}

//! Returns the form factor F(q) of this polyhedron, with origin at z=0.

complex_t ff::Polyhedron::formfactor_at_center(const C3& q) const
{
    const bool sym_Ci = m_topology->symmetry_Ci;
    const double q_red = m_radius * q.mag();
#ifdef ALGORITHM_DIAGNOSTIC
    polyhedralDiagnosis.reset();
#endif
    if (q_red == 0)
        return m_volume;
    if (q_red < q_limit_series) {
        // summation of power series
#ifdef ALGORITHM_DIAGNOSTIC
        polyhedralDiagnosis.algo = 100;
#endif
        complex_t sum = 0;
        complex_t n_fac = (sym_Ci ? -2 : -1) / q.mag2();
        int count_return_condition = 0;
        for (int n = 2; n < n_limit_series; ++n) {
            if (sym_Ci && n & 1)
                continue;
#ifdef ALGORITHM_DIAGNOSTIC
            polyhedralDiagnosis.order = std::max(polyhedralDiagnosis.order, n);
#endif
            complex_t term = 0;
            for (const Face& Gk : m_faces)
                term += Gk.ff_n(n + 1, q);
            term *= n_fac;
#ifdef ALGORITHM_DIAGNOSTIC_LEVEL2
            polyhedralDiagnosis.msg +=
                boost::str(boost::format("  + term(n=%2i) = %23.17e+i*%23.17e\n") % n % term.real()
                           % term.imag());
#endif
            sum += term;
            if (std::abs(term) <= eps * std::abs(sum) || std::abs(sum) < eps * m_volume)
                ++count_return_condition;
            else
                count_return_condition = 0;
            if (count_return_condition > 2)
                return m_volume + sum; // regular exit
            n_fac = sym_Ci ? -n_fac : mul_I(n_fac);
        }
        throw std::runtime_error("Numeric failure in polyhedron: series F(q) not converged");
    }

    // direct evaluation of analytic formula (coefficients may involve series)
#ifdef ALGORITHM_DIAGNOSTIC
    polyhedralDiagnosis.algo = 200;
#endif
    complex_t sum = 0;
    for (const Face& Gk : m_faces) {
        complex_t qn = Gk.normalProjectionConj(q); // conj(q)*normal
        if (std::abs(qn) < eps * q.mag())
            continue;
        complex_t term = qn * Gk.ff(q, sym_Ci);
#ifdef ALGORITHM_DIAGNOSTIC //_LEVEL2
        polyhedralDiagnosis.msg += boost::str(boost::format("  + face_ff = %23.17e+i*%23.17e\n")
                                              % term.real() % term.imag());
#endif
        sum += term;
    }
#ifdef ALGORITHM_DIAGNOSTIC //_LEVEL2
    polyhedralDiagnosis.msg +=
        boost::str(boost::format(" -> raw sum = %23.17e+i*%23.17e; divisor = %23.17e\n")
                   % sum.real() % sum.imag() % q.mag2());
#endif
    return sum / I / q.mag2();
}

//! Returns true if a vertex v is located inside the polygon by counting odd vs even ray cast
//! intersections. Uses fibonacci_sphere as sampling method and determines the result by majority
//! vote.

bool ff::Polyhedron::is_inside(const R3& v, int nsamples) const
{
    const R3 v_rel = v - location();
    if ((v_rel).mag() > m_radius)
        return false;
    if (nsamples < 2)
        throw std::runtime_error(
            "Invalid call to libformfactor, function is_inside: nsamples < 2 is "
            "out of range.");
    int n_inside = 0; // number of rays that cross an odd number of faces
    std::vector<R3> samples = fibonacci_sphere(v_rel, nsamples);
    for (const R3& sample : samples) {
        int crossings = ray_crossings(v_rel, sample, m_faces, false);
        if (m_topology->symmetry_Ci)
            crossings += ray_crossings(v_rel, sample, m_faces, true);
        if (crossings % 2 != 0)
            n_inside++;
    }
    return n_inside > nsamples / 2;
}

R3 ff::Polyhedron::center_of_mass() const
{
    if (m_topology->symmetry_Ci)
        return {};
    R3 result;
    for (const ff::Face& face : m_faces)
        result += face.center_of_polygon() * face.pyramidalVolume();
    return result * (3. / 4) / m_volume;
}

double ff::Polyhedron::z_bottom() const
{
    return std::min_element(m_vertices.begin(), m_vertices.end(),
                            [](const R3& a, const R3& b) { return a.z() < b.z(); })
        ->z();
}

//! Returns a clipped copy of this Polyhedron,
//! zMin is the lower clipping planes offset,
//! zMax is the upper clipping planes offset,
//! The new Polyhedron consists of the volume between zMin and zMax
ff::Polyhedron* ff::Polyhedron::clipped(double zMin, double zMax) const
{
    if (zMin > zMax)
        throw std::runtime_error("Invalid call to libformfactor, function clipped: zMin > zMax");
    if (zMin == zMax)
        throw std::runtime_error("Invalid call to libformfactor, function clipped: zMin == zMax");
    std::unique_ptr<Polyhedron> result(ff::atomic::z_clip(*this, zMin,false));
    result.reset(ff::atomic::z_clip(*result, zMax,true));
    return result.release();
}