//  ************************************************************************************************
//
//  BornAgain: simulate and fit reflection and scattering
//
//! @file      Sample/HardParticle/Polyhedra.cpp
//! @brief     Implements class Box.
//!
//! @homepage  http://www.bornagainproject.org
//! @license   GNU General Public License v3 or higher (see COPYING)
//! @copyright Forschungszentrum Jülich GmbH 2018
//! @authors   Scientific Computing Group at MLZ (see CITATION, AUTHORS)
//
//  ************************************************************************************************

#include "Sample/HardParticle/Polyhedra.h"
#include "Base/Const/Units.h"
#include <ff/Make.h>
#include <numbers>

using std::numbers::pi;

//  ************************************************************************************************
//  Prisms
//  ************************************************************************************************

Box::Box(const std::vector<double> P)
    : IFormfactorPrism(P)
    , m_length(m_P[0])
    , m_width(m_P[1])
    , m_height(m_P[2])
{
    pimpl.reset(ff::make::Box(m_length, m_width, m_height));
    m_validated = true;
}

bool Box::contains(const R3& position) const
{
    if (std::abs(position.x()) <= length() / 2 && std::abs(position.y()) <= width() / 2
        && (position.z() >= 0 && position.z() <= height()))
        return true;

    return false;
}

Box::Box(double length, double width, double height)
    : Box(std::vector<double>{length, width, height})
{
}


Prism3::Prism3(const std::vector<double> P)
    : IFormfactorPrism(P)
    , m_base_edge(m_P[0])
    , m_height(m_P[1])
{
    pimpl.reset(ff::make::Prism3(m_base_edge, m_height));
    m_validated = true;
}

bool Prism3::contains(const R3& position) const
{
    double B = baseEdge();
    double H = height();

    double l = B * std::sin(pi / 3);
    double x_shift = B / 2 * std::tan(pi / 6);

    if (position.x() + x_shift < 0 || position.x() + x_shift > l || std::abs(position.y()) > B / 2
        || position.z() < 0 || position.z() > H)
        return false;

    double theta = 0; // angle between position & x-axis in position.z() plane
    if (position.x() + x_shift != 0 || position.y() != 0)
        theta =
            std::asin(std::abs(position.y())
                      / std::sqrt(std::pow(position.x() + x_shift, 2) + std::pow(position.y(), 2)));

    double k = l / (std::sin(theta) / std::tan(pi / 6) + std::cos(theta));

    if (std::pow(position.x() + x_shift, 2) + std::pow(position.y(), 2) <= std::pow(k, 2))
        return true;

    return false;
}

Prism3::Prism3(double base_edge, double height)
    : Prism3(std::vector<double>{base_edge, height})
{
}


Prism6::Prism6(const std::vector<double> P)
    : IFormfactorPrism(P)
    , m_base_edge(m_P[0])
    , m_height(m_P[1])
{
    pimpl.reset(ff::make::Prism6(m_base_edge, m_height));
    m_validated = true;
}

bool Prism6::contains(const R3& position) const
{
    double B = baseEdge();
    double H = height();

    if (std::abs(position.x()) > B || std::abs(position.y()) > B || position.z() < 0
        || position.z() > H)
        return false;

    double theta_prime = 0; // angle between position & x-axis in position.z() plane
    if (position.x() != 0 || position.y() != 0)
        theta_prime = Units::rad2deg(
            std::asin(std::abs(position.y())
                      / std::sqrt(std::pow(position.x(), 2) + std::pow(position.y(), 2))));
    int c = static_cast<int>(theta_prime / 60); // multiplicative constant
    double theta = Units::deg2rad(theta_prime - c * 60);
    double k_z = B / (std::cos(theta) + std::sin(theta) / std::tan(pi / 3));

    if (std::pow(position.x(), 2) + std::pow(position.y(), 2) <= std::pow(k_z, 2))
        return true;

    return false;
}

Prism6::Prism6(double base_edge, double height)
    : Prism6(std::vector<double>{base_edge, height})
{
}

//  ************************************************************************************************
//  Platonic solids
//  ************************************************************************************************

PlatonicTetrahedron::PlatonicTetrahedron(const std::vector<double> P)
    : IFormfactorPolyhedron(P)
    , m_edge(m_P[0])
{
    pimpl.reset(ff::make::Tetrahedron(m_edge));
    m_validated = true;
}

bool PlatonicTetrahedron::contains(const R3&) const
{
    // TODO: Implement Platonic tetrahedron
    std::ostringstream ostr;
    ostr << "Outer shape Platonic tetrahedron not yet implemented for Mesocrystal";
    ostr << "\n\nStay tuned!";
    throw std::runtime_error(ostr.str());
}

PlatonicTetrahedron::PlatonicTetrahedron(double edge)
    : PlatonicTetrahedron(std::vector<double>{edge})
{
}


PlatonicOctahedron::PlatonicOctahedron(const std::vector<double> P)
    : IFormfactorPolyhedron(P)
    , m_edge(m_P[0])
{
    pimpl.reset(ff::make::Octahedron(m_edge));
    m_validated = true;
}

bool PlatonicOctahedron::contains(const R3&) const
{
    // TODO: Implement Platonic octahedron
    std::ostringstream ostr;
    ostr << "Outer shape Platonic octahedron not yet implemented for Mesocrystal";
    ostr << "\n\nStay tuned!";
    throw std::runtime_error(ostr.str());
}

PlatonicOctahedron::PlatonicOctahedron(double edge)
    : PlatonicOctahedron(std::vector<double>{edge})
{
}


Dodecahedron::Dodecahedron(const std::vector<double> P)
    : IFormfactorPolyhedron(P)
    , m_edge(m_P[0])
{
    pimpl.reset(ff::make::Dodecahedron(m_edge));
    m_validated = true;
}

bool Dodecahedron::contains(const R3&) const
{
    // TODO: Implement Dodecahedron
    std::ostringstream ostr;
    ostr << "Outer shape Dodecahedron not yet implemented for Mesocrystal";
    ostr << "\n\nStay tuned!";
    throw std::runtime_error(ostr.str());
}

Dodecahedron::Dodecahedron(double edge)
    : Dodecahedron(std::vector<double>{edge})
{
}


Icosahedron::Icosahedron(const std::vector<double> P)
    : IFormfactorPolyhedron(P)
    , m_edge(m_P[0])
{
    pimpl.reset(ff::make::Icosahedron(m_edge));
    m_validated = true;
}

bool Icosahedron::contains(const R3&) const
{
    // TODO: Implement Icosahedron
    std::ostringstream ostr;
    ostr << "Outer shape Icosahedron not yet implemented for Mesocrystal";
    ostr << "\n\nStay tuned!";
    throw std::runtime_error(ostr.str());
}

Icosahedron::Icosahedron(double edge)
    : Icosahedron(std::vector<double>{edge})
{
}

//  ************************************************************************************************
//  Pyramids
//  ************************************************************************************************

Pyramid2::Pyramid2(const std::vector<double> P)
    : IFormfactorPolyhedron(P)
    , m_length(m_P[0])
    , m_width(m_P[1])
    , m_height(m_P[2])
    , m_alpha(m_P[3])
{
    pimpl.reset(ff::make::Pyramid2(m_length, m_width, m_height, m_alpha));
    m_validated = true;
}

bool Pyramid2::contains(const R3& position) const
{
    double l_z = length() / 2
                 - position.z() / std::tan(alpha()); // half-length of rectangle at a given height
    double w_z =
        width() / 2 - position.z() / std::tan(alpha()); // half-width of rectangle at a given height
    if (std::abs(position.x()) <= l_z && std::abs(position.y()) <= w_z
        && (position.z() >= 0 && position.z() <= height()))
        return true;

    return false;
}

Pyramid2::Pyramid2(double length, double width, double height, double alpha)
    : Pyramid2(std::vector<double>{length, width, height, alpha})
{
}


Pyramid3::Pyramid3(const std::vector<double> P)
    : IFormfactorPolyhedron(P)
    , m_base_edge(m_P[0])
    , m_height(m_P[1])
    , m_alpha(m_P[2])
{
    pimpl.reset(ff::make::Pyramid3(m_base_edge, m_height, m_alpha));
    m_validated = true;
}

bool Pyramid3::contains(const R3& position) const
{
    double B = baseEdge();
    double H = height();

    double B_z = B - position.z() * 2 / std::tan(alpha()); // edge of triangle at a given height

    double l = B_z * std::sin(pi / 3);
    double x_shift = B_z / 2 * std::tan(pi / 6);

    if (position.x() + x_shift < 0 || position.x() + x_shift > l || std::abs(position.y()) > B_z / 2
        || position.z() < 0 || position.z() > H)
        return false;

    double theta = 0; // angle between position & x-axis in position.z() plane
    if (position.x() + x_shift != 0 || position.y() != 0)
        theta =
            std::asin(std::abs(position.y())
                      / std::sqrt(std::pow(position.x() + x_shift, 2) + std::pow(position.y(), 2)));

    double k = l / (std::sin(theta) / std::tan(pi / 6) + std::cos(theta));

    if (std::pow(position.x() + x_shift, 2) + std::pow(position.y(), 2) <= std::pow(k, 2))
        return true;

    return false;
}

Pyramid3::Pyramid3(double base_edge, double height, double alpha)
    : Pyramid3(std::vector<double>{base_edge, height, alpha})
{
}


Pyramid4::Pyramid4(const std::vector<double> P)
    : IFormfactorPolyhedron(P)
    , m_base_edge(m_P[0])
    , m_height(m_P[1])
    , m_alpha(m_P[2])
{
    pimpl.reset(ff::make::Pyramid4(m_base_edge, m_height, m_alpha));
    m_validated = true;
}

bool Pyramid4::contains(const R3& position) const
{
    double B = baseEdge();
    double H = height();

    double l_z =
        B / 2 - position.z() / std::tan(alpha()); // half-length of square at a given height
    if (std::abs(position.x()) <= l_z && std::abs(position.y()) <= l_z
        && (position.z() >= 0 && position.z() <= H))
        return true;

    return false;
}

Pyramid4::Pyramid4(double base_edge, double height, double alpha)
    : Pyramid4(std::vector<double>{base_edge, height, alpha})
{
}


Pyramid6::Pyramid6(const std::vector<double> P)
    : IFormfactorPolyhedron(P)
    , m_base_edge(m_P[0])
    , m_height(m_P[1])
    , m_alpha(m_P[2])
{
    pimpl.reset(ff::make::Pyramid6(m_base_edge, m_height, m_alpha));
    m_validated = true;
}

bool Pyramid6::contains(const R3& position) const
{
    double B = baseEdge();

    if (std::abs(position.x()) > B || std::abs(position.y()) > B || position.z() < 0
        || position.z() > height())
        return false;

    double l_z = B - position.z() / std::tan(alpha()); // edge of hexagon at a given height
    double theta_prime = 0; // angle between position & x-axis in position.z() plane
    if (position.x() != 0 || position.y() != 0)
        theta_prime = Units::rad2deg(
            std::asin(std::abs(position.y())
                      / std::sqrt(std::pow(position.x(), 2) + std::pow(position.y(), 2))));
    int c = static_cast<int>(theta_prime / 60); // multiplication constant
    double theta = Units::deg2rad(theta_prime - c * 60);
    double k_z = l_z / (std::cos(theta) + std::sin(theta) / std::tan(pi / 3));

    if (std::pow(position.x(), 2) + std::pow(position.y(), 2) <= std::pow(k_z, 2))
        return true;

    return false;
}

Pyramid6::Pyramid6(double base_edge, double height, double alpha)
    : Pyramid6(std::vector<double>{base_edge, height, alpha})
{
}


Bipyramid4::Bipyramid4(const std::vector<double> P)
    : IFormfactorPolyhedron(P)
    , m_length(m_P[0])
    , m_base_height(m_P[1])
    , m_height_ratio(m_P[2])
    , m_alpha(m_P[3])
{
    pimpl.reset(ff::make::Bipyramid4(m_length, m_base_height, m_height_ratio, m_alpha));
    m_validated = true;
}

bool Bipyramid4::contains(const R3& position) const
{
    double L = length();
    double H = base_height();
    double rH = heightRatio();

    double total_Height = H + rH * H;

    if (std::abs(position.x()) > L / 2 || std::abs(position.y()) > L / 2 || position.z() < 0
        || position.z() > total_Height)
        return false;

    // half-length of square (i.e. horizontal cross-section of Bipyramid4) at a given height
    double l_z = L / 2 - std::abs(H - position.z()) / std::tan(alpha());
    if (std::abs(position.x()) <= l_z && std::abs(position.y()) <= l_z)
        return true;

    return false;
}

Bipyramid4::Bipyramid4(double length, double base_height, double height_ratio, double alpha)
    : Bipyramid4(std::vector<double>{length, base_height, height_ratio, alpha})
{
}

//  ************************************************************************************************
//  Others
//  ************************************************************************************************

CantellatedCube::CantellatedCube(const std::vector<double> P)
    : IFormfactorPolyhedron(P)
    , m_length(m_P[0])
    , m_removed_length(m_P[1])
{
    pimpl.reset(ff::make::CantellatedCube(m_length, m_removed_length));
    m_validated = true;
}

bool CantellatedCube::contains(const R3&) const
{
    // TODO: Implement Cantellated cube
    std::ostringstream ostr;
    ostr << "Outer shape Cantellated cube not yet implemented for Mesocrystal";
    ostr << "\n\nStay tuned!";
    throw std::runtime_error(ostr.str());
}

CantellatedCube::CantellatedCube(double length, double removed_length)
    : CantellatedCube(std::vector<double>{length, removed_length})
{
}


TruncatedCube::TruncatedCube(const std::vector<double> P)
    : IFormfactorPolyhedron(P)
    , m_length(m_P[0])
    , m_removed_length(m_P[1])
{
    pimpl.reset(ff::make::TruncatedCube(m_length, m_removed_length));
    m_validated = true;
}

bool TruncatedCube::contains(const R3&) const
{
    // TODO: Implement Truncated cube
    std::ostringstream ostr;
    ostr << "Outer shape Truncated cube not yet implemented for Mesocrystal";
    ostr << "\n\nStay tuned!";
    throw std::runtime_error(ostr.str());
}

TruncatedCube::TruncatedCube(double length, double removed_length)
    : TruncatedCube(std::vector<double>{length, removed_length})
{
}