#pragma once
// Copyright (c) 2024-2026 Tarik Moussa.
// SPDX-License-Identifier: MIT


// Hyperbolic tetrahedron volume formulas.
// Ported from de.varylab.discreteconformal.functional.HyperIdealUtility (Java).

#include "clausen.hpp"
#include "constants.hpp"

#include <Eigen/Dense>
#include <cmath>
#include <complex>

namespace conformallab {

/// Volume of a generalized hyperbolic tetrahedron with dihedral
/// angles `A,…,F` via the Ushijima 2006 formula (DOI 10.1007/0-387-29555-0_13,
/// arxiv math/0309216). Note: sole author is Ushijima; "Meyerhoff" is not an author.
/// Same as Java `HyperIdealUtility.calculateTetrahedronVolume()`.
inline double calculateTetrahedronVolume(double A, double B, double C,
                                         double D, double E, double F) {
    // PI from constants.hpp (conformallab::PI)

    // Degenerate if any angle equals pi.
    if (A == PI || B == PI || C == PI || D == PI || E == PI || F == PI)
        return 0.0;

    const double sA = std::sin(A), sB = std::sin(B), sC = std::sin(C);
    const double sD = std::sin(D), sE = std::sin(E), sF = std::sin(F);
    const double cA = std::cos(A), cB = std::cos(B), cC = std::cos(C);
    const double cD = std::cos(D), cE = std::cos(E), cF = std::cos(F);

    // Unit complex numbers e^(i*angle).
    using Cx = std::complex<double>;
    auto polar = [](double angle) { return std::polar(1.0, angle); };

    Cx ad = polar(A + D), be = polar(B + E), cf = polar(C + F);
    Cx abc = polar(A + B + C), abf = polar(A + B + F);
    Cx ace = polar(A + C + E), aef = polar(A + E + F);
    Cx bcd = polar(B + C + D), bdf = polar(B + D + F);
    Cx def = polar(D + E + F), cde = polar(C + D + E);
    Cx abde = ad * be, acdf = ad * cf, bcef = be * cf;
    Cx abcdef = abc * def;

    Cx z = ad + be + cf + abf + ace + bcd + def + abcdef;

    // Gram matrix of the tetrahedron.
    Eigen::Matrix4d G;
    G <<  1.0, -cA, -cB, -cF,
         -cA,  1.0, -cC, -cE,
         -cB, -cC,  1.0, -cD,
         -cF, -cE,  -cD,  1.0;

    Cx sqrtG = std::sqrt(Cx(G.determinant(), 0.0));

    Cx f   = Cx(sA*sD + sB*sE + sC*sF, 0.0);
    Cx f1  = f - sqrtG;
    Cx f2  = f + sqrtG;
    Cx z1  = -2.0 * f1 / z;
    Cx z2  = -2.0 * f2 / z;

    auto U = [&](Cx zi) {
        return 0.5 * (
            + ImLi2(zi)
            + ImLi2(abde  * zi)
            + ImLi2(acdf  * zi)
            + ImLi2(bcef  * zi)
            - ImLi2(-abc  * zi)
            - ImLi2(-aef  * zi)
            - ImLi2(-bdf  * zi)
            - ImLi2(-cde  * zi)
        );
    };

    return (U(z1) - U(z2)) / 2.0;
}

/// Volume of a hyperideal tetrahedron with one ideal vertex at γ via
/// the Springborn 2008 formula (arxiv math/0603097).  Same as Java
/// `HyperIdealUtility.calculateTetrahedronVolumeWithIdealVertexAtGamma()`.
inline double calculateTetrahedronVolumeWithIdealVertexAtGamma(
    double gamma1, double gamma2, double gamma3,
    double alpha23, double alpha31, double alpha12)
{
    // PI from constants.hpp (conformallab::PI)
    auto L = [](double x) { return Lobachevsky(x); };

    double result = L(gamma1) + L(gamma2) + L(gamma3);

    result += L((PI + alpha31 - alpha12 - gamma1) / 2.0);
    result += L((PI + alpha12 - alpha23 - gamma2) / 2.0);
    result += L((PI + alpha23 - alpha31 - gamma3) / 2.0);

    result += L((PI - alpha31 + alpha12 - gamma1) / 2.0);
    result += L((PI - alpha12 + alpha23 - gamma2) / 2.0);
    result += L((PI - alpha23 + alpha31 - gamma3) / 2.0);

    result += L((PI + alpha31 + alpha12 - gamma1) / 2.0);
    result += L((PI + alpha12 + alpha23 - gamma2) / 2.0);
    result += L((PI + alpha23 + alpha31 - gamma3) / 2.0);

    result += L((PI - alpha31 - alpha12 - gamma1) / 2.0);
    result += L((PI - alpha12 - alpha23 - gamma2) / 2.0);
    result += L((PI - alpha23 - alpha31 - gamma3) / 2.0);

    return result / 2.0;
}

} // namespace conformallab