//  ************************************************************************************************
//
//  libformfactor: efficient and accurate computation of scattering form factors
//
//! @file      ff/Face.cpp
//! @brief     Implements class Face
//!
//! @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)
//
//  ************************************************************************************************

#include "ff/Face.h"
#include "ff/Edge.h"
#include "ff/Factorial.h"
#include <algorithm>
#include <functional>
#include <stdexcept>
#define _USE_MATH_DEFINES
#include <math.h>

namespace {

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

constexpr auto ReciprocalFactorialArray = ff_aux::generateReciprocalFactorialArray<171>();

complex_t sinc(const complex_t z) // cardinal sine function, sin(x)/x
{
    // This is an exception from the rule that we must not test floating-point numbers for equality.
    // For small non-zero arguments, sin(z) returns quite accurately z or z-z^3/6.
    // There is no loss of precision in computing sin(z)/z.
    // Therefore there is no need for an expensive test like abs(z)<eps.
    if (z == complex_t(0., 0.))
        return 1.0;
    return std::sin(z) / z;
}

std::vector<ff::Edge> precompute_edges_raw(const std::vector<R3>& V, double radius_2d)
{
    size_t NV = V.size();
    if (!NV)
        throw std::runtime_error("Invalid polyhedral face: no edges given");
    if (NV < 3)
        throw std::runtime_error("Invalid polyhedral face: less than three edges");
    // Initialize list of 'edges'.
    // Do not create an edge if two vertices are too close to each other.
    // TODO This is implemented in a somewhat sloppy way: we just skip an edge if it would
    //      be too short. This leaves tiny open edges. In a clean implementation, we
    //      rather should merge adjacent vertices before generating edges.
    std::vector<ff::Edge> result;
    for (size_t j = 0; j < NV; ++j) {
        size_t jj = (j + 1) % NV;
        if ((V[j] - V[jj]).mag() < 1e-14 * radius_2d)
            continue; // distance too short -> skip this edge
        result.emplace_back(V[j], V[jj]);
    }
    return result;
}

std::vector<ff::Edge> precompute_edges_with_symmetry(const std::vector<R3>& V, bool sym_S2,
                                                     const R3& face_center, double radius_2d)
{
    // only now deal with inversion symmetry
    std::vector<ff::Edge> result = precompute_edges_raw(V, radius_2d);
    size_t NE = result.size();
    if (sym_S2) {
        if (NE & 1)
            throw std::runtime_error("Invalid polyhedral face: odd #edges violates symmetry S2");
        NE /= 2;
        for (size_t j = 0; j < NE; ++j) {
            if ((result[j].R() + result[j + NE].R() - 2 * face_center).mag2()
                > pow(2e-12 * radius_2d, 2))
                throw std::runtime_error(
                    "Invalid polyhedral face: edge centers violate symmetry S2");
            if ((result[j].E() + result[j + NE].E()).mag() > 1e-12 * radius_2d)
                throw std::runtime_error(
                    "Invalid polyhedral face: edge vectors violate symmetry S2");
        }
        // keep only half of the egdes
        result.erase(result.begin() + NE, result.end());
    }
    return result;
}

R3 precompute_normal(const std::vector<ff::Edge>& edges)
{
    size_t NE = edges.size();
    if (NE < 3)
        throw std::invalid_argument("Face has less than three non-vanishing edges");
    R3 result;
    for (size_t j = 0; j < NE; ++j)
        result += edges[j].E().cross(edges[(j + 1) % NE].E());
    return result.unit_or_null();
}

//! Returns distance of figure plane from origin
double precompute_rperp(const std::vector<R3>& V, const R3& normal, double radius_2d)
{
    const size_t NV = V.size();
    const double estimate = V[0].dot(normal);
    double sum1 = 0;
    for (size_t j = 1; j < NV; ++j) {
        const double x = V[j].dot(normal);
        sum1 += (x - estimate);
    }
    const double rperp = estimate + sum1 / NV;
    return rperp;
}

double precompute_area(const std::vector<R3>& V, const R3& normal)
{
    double result = 0;
    size_t NV = V.size();
    // compute m_area
    for (size_t j = 0; j < NV; ++j) {
        size_t jj = (j + 1) % NV;
        result += normal.dot(V[j].cross(V[jj])) / 2;
    }
    return result;
}

//! Returns radius of circle that contains all vertices.
double precompute_radius(const std::vector<R3>& V)
{
    double diameter = 0;
    for (size_t j = 0; j < V.size(); ++j)
        for (size_t jj = j + 1; jj < V.size(); ++jj)
            diameter = std::max(diameter, (V[j] - V[jj]).mag());
    return diameter / 2;
}

} // namespace


#ifdef ALGORITHM_DIAGNOSTIC
void ff::PolyhedralDiagnosis::reset()
{
    order = 0;
    algo = 0;
    msg.clear();
};
std::string ff::PolyhedralDiagnosis::message() const
{
    std::string result = "algo=" + std::to_string(algo) + ", order=" + std::to_string(order);
    if (!msg.empty())
        result += ", msg:\n" + msg;
    return result;
}
bool ff::PolyhedralDiagnosis::operator==(const PolyhedralDiagnosis& other) const
{
    return order == other.order && algo == other.algo;
}
bool ff::PolyhedralDiagnosis::operator!=(const PolyhedralDiagnosis& other) const
{
    return !(*this == other);
}
#endif

//  ************************************************************************************************
//  PolyhedralFace implementation
//  ************************************************************************************************

//! Sets internal variables for given vertex chain.

//! @param V oriented vertex list
//! @param _sym_S2 true if face has a perpedicular two-fold symmetry axis

ff::Face::Face(const std::vector<R3>& V, bool sym_S2)
    : m_sym_S2(sym_S2)
    , m_radius_2d(precompute_radius(V))
    , m_normal(precompute_normal(precompute_edges_raw(V, m_radius_2d)))
    , m_rperp(precompute_rperp(V, m_normal, m_radius_2d))
    , m_area(precompute_area(V, m_normal))
    , m_edges(precompute_edges_with_symmetry(V, m_sym_S2, m_rperp * m_normal, m_radius_2d))
{
}

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

const std::vector<ff::Edge>& ff::Face::edges() const
{
    return m_edges;
}

R3 ff::Face::center_of_polygon() const
{
    if (m_sym_S2)
        return m_normal * m_rperp;
    const R3 scaled_normal = (2 / (3. * m_area)) * m_normal;
    R3 meanR;
    for (const Edge& e : m_edges)
        meanR += e.R();
    meanR /= m_edges.size();
    R3 result = meanR;
    for (const Edge& e : m_edges)
        result += (e.R() - meanR).cross(e.E()).dot(scaled_normal) * (e.R() - meanR);
    return result;
}

//! Sets qperp and qpa according to argument q and to this polygon's normal.

void ff::Face::decompose_q(C3 q, complex_t& qperp, C3& qpa) const
{
    qperp = m_normal.dot(q);
    qpa = q - qperp * m_normal;
    // improve numeric accuracy:
    qpa -= m_normal.dot(qpa) * m_normal;
    if (qpa.mag() < eps * std::abs(qperp))
        qpa = C3(0., 0., 0.);
}

//! Returns core contribution to f_n

complex_t ff::Face::ff_n_core(int m, C3 qpa, complex_t qperp) const
{
    const C3 prevec = 2. * m_normal.cross(qpa); // complex conjugation not here but in .dot
    complex_t result = 0;
    const complex_t qrperp = qperp * m_rperp;
    for (size_t i = 0; i < m_edges.size(); ++i) {
        const Edge& e = m_edges[i];
        const complex_t vfac = prevec.dot(e.E());
        const complex_t tmp = e.contrib(m + 1, qpa, qrperp);
        result += vfac * tmp;
        //     std::cout << std::scientific << std::showpos << std::setprecision(16)
        //               << "DBX ff_n_core " << m << " " << vfac << " " << tmp
        //               << " term=" << vfac * tmp << " sum=" << result << "\n";
    }
    return result;
}

//! Returns contribution qn*f_n [of order q^(n+1)] from this face to the polyhedral form factor.

complex_t ff::Face::ff_n(int n, C3 q) const
{
    complex_t qn = q.dot(m_normal); // conj(q)*normal (dot is antilinear in 'this' argument)
    if (std::abs(qn) < eps * q.mag())
        return 0.;
    complex_t qperp;
    C3 qpa;
    decompose_q(q, qperp, qpa);
    double qpa_mag2 = qpa.mag2();
    if (qpa_mag2 == 0.)
        return qn * pow(qperp * m_rperp, n) * m_area * ReciprocalFactorialArray[n];
    if (m_sym_S2)
        return qn * (ff_n_core(n, qpa, qperp) + ff_n_core(n, -qpa, qperp)) / qpa_mag2;
    complex_t tmp = ff_n_core(n, qpa, qperp);
    // std::cout << "DBX ff_n " << n << " " << qn << " " << tmp << " " << qpa_mag2 << "\n";
    return qn * tmp / qpa_mag2;
}

//! Returns sum of n>=1 terms of qpa expansion of 2d form factor

complex_t ff::Face::expansion(complex_t fac_even, complex_t fac_odd, C3 qpa, double abslevel) const
{
#ifdef ALGORITHM_DIAGNOSTIC
    polyhedralDiagnosis.algo += 1;
#endif
    complex_t sum = 0;
    complex_t n_fac = I;
    int count_return_condition = 0;
    for (int n = 1; n < n_limit_series; ++n) {
#ifdef ALGORITHM_DIAGNOSTIC
        polyhedralDiagnosis.order = std::max(polyhedralDiagnosis.order, n);
#endif
        complex_t term = n_fac * (n & 1 ? fac_odd : fac_even) * ff_n_core(n, qpa, 0) / qpa.mag2();
        // std::cout << std::setprecision(16) << "    sum=" << sum << " +term=" << term << "\n";
        sum += term;
        if (std::abs(term) <= eps * std::abs(sum) || std::abs(sum) < eps * abslevel)
            ++count_return_condition;
        else
            count_return_condition = 0;
        if (count_return_condition > 2)
            return sum; // regular exit
        n_fac = mul_I(n_fac);
    }
    throw std::runtime_error("Numeric error in polyhedral face: series f(q_pa) not converged");
}

//! Returns core contribution to analytic 2d form factor.

complex_t ff::Face::edge_sum_ff(C3 q, C3 qpa, bool sym_Ci) const
{
    C3 prevec = m_normal.cross(qpa); // complex conjugation will take place in .dot
    complex_t sum = 0;
    complex_t vfacsum = 0;
    for (size_t i = 0; i < m_edges.size(); ++i) {
        const ff::Edge& e = m_edges[i];
        complex_t qE = e.qE(qpa);
        complex_t qR = e.qR(qpa);
        complex_t Rfac = m_sym_S2 ? sin(qR) : (sym_Ci ? cos(e.qR(q)) : exp_I(qR));
        complex_t vfac;
        if (m_sym_S2 || i < m_edges.size() - 1) {
            vfac = prevec.dot(e.E());
            vfacsum += vfac;
        } else {
            vfac = -vfacsum; // to improve numeric accuracy: qcE_J = - sum_{j=0}^{J-1} qcE_j
        }
        complex_t term = vfac * sinc(qE) * Rfac;
        sum += term;
        //    std::cout << std::scientific << std::showpos << std::setprecision(16)
        //              << "    sum=" << sum << " term=" << term << " vf=" << vfac << " qE=" << qE
        //              << " qR=" << qR << " sinc=" << sinc(qE) << " Rfac=" << Rfac << "\n";
    }
    return sum;
}

//! Returns the contribution ff(q) of this face to the polyhedral form factor.

complex_t ff::Face::ff(C3 q, bool sym_Ci) const
{
    complex_t qperp;
    C3 qpa;
    decompose_q(q, qperp, qpa);
    double qpa_red = m_radius_2d * qpa.mag();
    complex_t qr_perp = qperp * m_rperp;
    complex_t ff0 = (sym_Ci ? 2. * I * sin(qr_perp) : exp_I(qr_perp)) * m_area;
    if (qpa_red == 0)
        return ff0;
    if (qpa_red < qpa_limit_series && !m_sym_S2) {
        // summation of power series
        complex_t fac_even;
        complex_t fac_odd;
        if (sym_Ci) {
            fac_even = 2. * mul_I(sin(qr_perp));
            fac_odd = 2. * cos(qr_perp);
        } else {
            fac_even = exp_I(qr_perp);
            fac_odd = fac_even;
        }
        return ff0 + expansion(fac_even, fac_odd, qpa, std::abs(ff0));
    }
    // direct evaluation of analytic formula
    complex_t prefac;
    if (m_sym_S2)
        prefac = sym_Ci ? -8. * sin(qr_perp) : 4. * mul_I(exp_I(qr_perp));
    else
        prefac = sym_Ci ? 4. : 2. * exp_I(qr_perp);
    // std::cout << "       qrperp=" << qr_perp << " => prefac=" << prefac << "\n";
    return prefac * edge_sum_ff(q, qpa, sym_Ci) / mul_I(qpa.mag2());
}

//! Two-dimensional form factor, for use in prism, from power series.

complex_t ff::Face::ff_2D_expanded(C3 qpa) const
{
    return m_area + expansion(1., 1., qpa, std::abs(m_area));
}

//! Two-dimensional form factor, for use in prism, from sum over edge form factors.

complex_t ff::Face::ff_2D_direct(C3 qpa) const
{
    return (m_sym_S2 ? 4. : 2. / I) * edge_sum_ff(qpa, qpa, false) / qpa.mag2();
}

//! Returns the two-dimensional form factor of this face, for use in a prism.

complex_t ff::Face::ff_2D(C3 qpa) const
{
    if (std::abs(qpa.dot(m_normal)) > eps * qpa.mag())
        throw std::runtime_error(
            "Numeric error in polyhedral formfactor: ff_2D called with perpendicular q component");
    double qpa_red = m_radius_2d * qpa.mag();
    if (qpa_red == 0)
        return m_area;
    if (qpa_red < qpa_limit_series && !m_sym_S2)
        return ff_2D_expanded(qpa);
    return ff_2D_direct(qpa);
}

//! Throws if deviation from inversion symmetry is detected. Does not check vertices.

void ff::Face::assert_Ci(const Face& other) const
{
    if (std::abs(m_rperp - other.m_rperp) > std::abs(1e-15 * (m_rperp + other.m_rperp)))
        throw std::runtime_error(
            "Invalid polyhedron: faces with different distance from origin violate symmetry Ci");
    if (std::abs(m_area - other.m_area) > 1e-15 * (m_area + other.m_area))
        throw std::runtime_error(
            "Invalid polyhedron: faces with different areas violate symmetry Ci");
    if ((m_normal + other.m_normal).mag() > 1e-14)
        throw std::runtime_error(
            "Invalid polyhedron: faces do not have opposite orientation, violating symmetry Ci");
}

//! Returns true if point v, is located inside the infinite extruded volume defined by this polygon.

bool ff::Face::is_inside(const R3& v) const
{
    int n_rightOf = 0; // number of occasions v is located right of the edge
    for (const Edge& edge : m_edges) {
        R3 p = v - edge.R();
        double o = (m_normal.cross(edge.E())).dot(p);
        if (o >= 0)
            n_rightOf++;
        else if (o < 0)
            n_rightOf--;
        if (m_sym_S2) {
            p = v + edge.R();
            double o = (m_normal.cross(-edge.E())).dot(p);
            if (o >= 0)
                n_rightOf++;
            else if (o < 0)
                n_rightOf--;
        }
    }
    int edge_size = m_sym_S2 ? m_edges.size() * 2 : m_edges.size();
    return abs(n_rightOf) == edge_size;
}