//  ************************************************************************************************
//
//  BornAgain: simulate and fit reflection and scattering
//
//! @file      Resample/Processed/Slicer.cpp
//! @brief     Implements function Compute::Slicing::sliceFormfactor.
//!
//! @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 "Resample/Processed/Slicer.h"
#include "Base/Math/Functions.h"
#include "Base/Type/Span.h"
#include "Base/Util/Assert.h"
#include "Base/Vector/RotMatrix.h"
#include "Resample/Particle/ReCompound.h"
#include "Resample/Particle/ReMesocrystal.h"
#include "Resample/Particle/ReParticle.h"
#include "Resample/Slice/ZLimits.h"
#include "Sample/HardParticle/HardParticles.h"
#include "Sample/Material/MaterialUtil.h"
#include "Sample/Particle/Compound.h"
#include "Sample/Particle/CoreAndShell.h"
#include "Sample/Particle/Crystal.h"
#include "Sample/Particle/Mesocrystal.h"
#include "Sample/Particle/Particle.h"
#include "Sample/Scattering/Rotations.h"
#include <numbers>

using std::numbers::pi;

namespace {

IFormfactor* doSlice(const IFormfactor* ff, double dz_bottom, double dz_top)
{
    if (const auto* f = dynamic_cast<const Pyramid2*>(ff)) {
        double dbase_edge = 2 * dz_bottom * Math::cot(f->alpha());
        return new Pyramid2(f->length() - dbase_edge, f->width() - dbase_edge,
                            f->height() - dz_bottom - dz_top, f->alpha());
    }
    if (const auto* f = dynamic_cast<const Box*>(ff)) {
        return new Box(f->length(), f->width(), f->height() - dz_bottom - dz_top);
    }
    if (const auto* f = dynamic_cast<const Cone*>(ff)) {
        double dradius = dz_bottom * Math::cot(f->alpha());
        return new Cone(f->radius() - dradius, f->height() - dz_bottom - dz_top, f->alpha());
    }
    if (const auto* f = dynamic_cast<const Pyramid6*>(ff)) {
        double dbase_edge = 2 / sqrt(3) * dz_bottom * Math::cot(f->alpha());
        return new Pyramid6(f->baseEdge() - dbase_edge, f->height() - dz_bottom - dz_top,
                            f->alpha());
    }
    if (const auto* f = dynamic_cast<const Bipyramid4*>(ff)) {
        if (dz_bottom > f->base_height()) {
            double dbase_edge = 2 * (dz_bottom - f->base_height()) * Math::cot(f->alpha());
            return new Pyramid4(f->length() - dbase_edge, f->height() - dz_bottom - dz_top,
                                f->alpha());
        }
        if (dz_top > f->heightRatio() * f->base_height()) {
            double dbase_edge = 2 * (f->base_height() - dz_bottom) * Math::cot(f->alpha());
            return new Pyramid4(f->length() - dbase_edge, f->height() - dz_bottom - dz_top,
                                pi - f->alpha());
        }
        {
            return new Bipyramid4(f->length(), f->base_height() - dz_bottom,
                                  f->heightRatio() * f->base_height() - dz_top, f->alpha());
        }
    }
    if (const auto* f = dynamic_cast<const Cylinder*>(ff)) {
        return new Cylinder(f->radius(), f->height() - dz_bottom - dz_top);
    }
    if (const auto* f = dynamic_cast<const EllipsoidalCylinder*>(ff)) {
        return new EllipsoidalCylinder(f->radiusX(), f->radiusY(),
                                       f->height() - dz_bottom - dz_top);
    }
    if (const auto* f = dynamic_cast<const HorizontalCylinder*>(ff)) {
        return new HorizontalCylinder(f->radius(), f->length(), f->slice_bottom() + dz_bottom,
                                      f->slice_top() - dz_top);
    }
    if (const auto* f = dynamic_cast<const Sphere*>(ff)) {
        return new SphericalSegment(f->radius(), dz_top, dz_bottom);
    }
    if (const auto* f = dynamic_cast<const Spheroid*>(ff)) {
        return new SpheroidalSegment(f->radiusXY(), f->radiusZ(), dz_top, dz_bottom);
    }
    if (const auto* f = dynamic_cast<const LongBoxGauss*>(ff)) {
        return new LongBoxGauss(f->length(), f->width(), f->height() - dz_bottom - dz_top);
    }
    if (const auto* f = dynamic_cast<const LongBoxLorentz*>(ff)) {
        return new LongBoxLorentz(f->length(), f->width(), f->height() - dz_bottom - dz_top);
    }
    if (const auto* f = dynamic_cast<const Prism3*>(ff)) {
        return new Prism3(f->baseEdge(), f->height() - dz_bottom - dz_top);
    }
    if (const auto* f = dynamic_cast<const Prism6*>(ff)) {
        return new Prism6(f->baseEdge(), f->height() - dz_bottom - dz_top);
    }
    if (const auto* f = dynamic_cast<const Pyramid4*>(ff)) {
        double dbaseEdge = 2 * dz_bottom * Math::cot(f->alpha());
        return new Pyramid4(f->baseEdge() - dbaseEdge, f->height() - dz_bottom - dz_top,
                            f->alpha());
    }
    if (const auto* f = dynamic_cast<const Pyramid3*>(ff)) {
        double dbaseEdge = 2 * sqrt(3) * dz_bottom * Math::cot(f->alpha());
        return new Pyramid3(f->baseEdge() - dbaseEdge, f->height() - dz_bottom - dz_top,
                            f->alpha());
    }
    if (const auto* f = dynamic_cast<const SphericalSegment*>(ff)) {
        return new SphericalSegment(f->radius(), f->cutFromTop() + dz_top,
                                    f->cutFromBottom() + dz_bottom);
    }
    if (const auto* f = dynamic_cast<const SpheroidalSegment*>(ff)) {
        return new SpheroidalSegment(f->radiusXY(), f->radiusZ(), f->cutFromTop() + dz_top,
                                     f->cutFromBottom() + dz_bottom);
    }
    throw std::runtime_error("Slicing of " + ff->className() + " not supported");
}

ReParticle* createParticleSlice(const IFormfactor* ff, ZLimits limits, const R3& position,
                                const IRotation* rot)
{
    const Span span = ff->spanZ(rot) + position.z();
    if (span.hig() <= limits.low() || span.low() >= limits.hig())
        return nullptr;

    const RotMatrix* rotMatrix =
        rot && !rot->isIdentity() ? new RotMatrix(rot->rotMatrix()) : nullptr;

    if (limits.low() <= span.low() && span.hig() <= limits.hig())
        return new ReParticle(ff->clone(), new R3(position), rotMatrix);

    if (!ff->canSliceAnalytically(rot))
        throw std::runtime_error("Slicing of " + ff->className()
                                 + " not supported for the given rotation");

    const double height = span.hig() - span.low();
    R3 new_position(position);
    double z_bottom = position.z();
    double z_top = position.z() + height;
    double dz_top = std::isinf(limits.hig()) ? -1 : z_top - limits.hig();
    double dz_bottom = std::isinf(limits.low()) ? -1 : limits.low() - z_bottom;
    ASSERT(dz_top >= 0 || dz_bottom >= 0);
    ASSERT(dz_bottom <= height);
    ASSERT(dz_top <= height);
    if (dz_bottom < 0)
        dz_bottom = 0;
    if (dz_top < 0)
        dz_top = 0;
    if (dz_bottom > 0)
        new_position.setZ(limits.low());

    IFormfactor* slicedff = doSlice(ff, dz_bottom, dz_top);
    return new ReParticle(slicedff, new R3(new_position), rotMatrix);
}

//! Recursively processes the basis of a mesocrystal.
IReParticle* processBasis(const IParticle* basis, const Material& ambientMat)
{
    if (const auto* b = dynamic_cast<const Compound*>(basis)) {
        const auto& particles = b->decompose();
        ASSERT(!particles.empty());
        auto* result = new ReCompound;
        for (const auto* particle : particles) {
            std::unique_ptr<IReParticle> re(processBasis(particle, ambientMat));
            result->addFormfactor(*re);
        }
        return result;
    } else if (dynamic_cast<const CoreAndShell*>(basis)) {
        throw std::runtime_error("Mesocrystal with CoreAndShell basis not yet supported");
    }
    // Remaining case: the basis is a single particle.
    const auto* p = dynamic_cast<const Particle*>(basis);
    ASSERT(p);
    const IRotation* rot = p->rotation();
    auto particleSlice = std::make_unique<ReParticle>(
        p->pFormfactor()->clone(), new const R3(p->particlePosition()),
        rot && !rot->isIdentity() ? new const RotMatrix(rot->rotMatrix()) : nullptr);
    if (!particleSlice)
        return {};
    double volume = particleSlice->volume();
    Material transformed_material(p->rotation()
                                      ? p->material()->rotatedMaterial(p->rotation()->rotMatrix())
                                      : *p->material());
    particleSlice->setMaterial(transformed_material);
    particleSlice->setAdmixedFraction(volume);
    particleSlice->setAdmixedMaterial(transformed_material);
    particleSlice->setAmbientMaterial(ambientMat);
    return particleSlice.release();
}

} // namespace


OwningVector<IReParticle> Compute::Slicing::particlesInSlice(const IParticle* particle,
                                                             const ZLimits& limits,
                                                             const Material& ambientMat,
                                                             const MesoOptions& meso_options)
{
    if (const auto* p = dynamic_cast<const Particle*>(particle)) {
        ASSERT(p->pFormfactor());
        std::unique_ptr<ReParticle> particleSlice(
            createParticleSlice(p->pFormfactor(), limits, p->particlePosition(), p->rotation()));
        if (!particleSlice)
            return {};
        double volume = particleSlice->volume();
        Material transformed_material(
            p->rotation() ? p->material()->rotatedMaterial(p->rotation()->rotMatrix())
                          : *p->material());
        particleSlice->setMaterial(transformed_material);
        particleSlice->setAdmixedFraction(volume);
        particleSlice->setAdmixedMaterial(transformed_material);
        particleSlice->setAmbientMaterial(ambientMat);
        OwningVector<IReParticle> result;
        result.push_back(particleSlice.release());
        return result;
    }

    if (const auto* p = dynamic_cast<const CoreAndShell*>(particle)) {
        const Particle* core = p->coreParticle();
        const Particle* shell = p->shellParticle();
        ASSERT(core && shell);

        // shell
        std::unique_ptr<Particle> P_shell(shell->clone());
        if (p->rotation())
            P_shell->rotate(*p->rotation());
        P_shell->translate(p->particlePosition());
        OwningVector<IReParticle> shell_slices =
            particlesInSlice(P_shell.get(), limits, ambientMat, meso_options);
        if (shell_slices.empty())
            return {};
        ASSERT(shell_slices.size() == 1);
        IReParticle* sliced_shell = shell_slices.release_back();
        if (!sliced_shell)
            return {};
        const Material& shell_material = sliced_shell->admixed().material;

        // core
        std::unique_ptr<Particle> P_core(core->clone());
        if (p->rotation())
            P_core->rotate(*p->rotation());
        P_core->translate(p->particlePosition());
        OwningVector<IReParticle> core_slices =
            particlesInSlice(P_core.get(), limits, shell_material, meso_options);
        ASSERT(core_slices.size() == 1);
        IReParticle* sliced_core = core_slices.release_back();

        // if core out of limits, return sliced shell
        if (!sliced_core) {
            OwningVector<IReParticle> result;
            result.push_back(sliced_shell);
            return result;
        }

        sliced_shell->setAdmixedFraction(sliced_shell->admixedFraction()
                                         - sliced_core->admixedFraction());

        OwningVector<IReParticle> result;
        result.push_back(sliced_core);
        result.push_back(sliced_shell);
        return result;
    }

    if (dynamic_cast<const Compound*>(particle))
        throw std::runtime_error("Compound does not yet support slicing");

    if (const auto* p = dynamic_cast<const Mesocrystal*>(particle)) {
        const Crystal* crystal = &p->particleStructure();
        const IFormfactor* meso_formfactor = p->outerShape();
        ASSERT(crystal && meso_formfactor);

        double unit_cell_volume = crystal->lattice()->unitCellVolume();
        if (unit_cell_volume <= 0)
            return {};

        std::unique_ptr<ReParticle> new_shape(
            createParticleSlice(meso_formfactor, limits, p->particlePosition(), p->rotation()));

        const std::unique_ptr<Crystal> new_crystal(
            crystal->transformed(p->particlePosition(), p->rotation()));

        const std::unique_ptr<IReParticle> new_basis(
            processBasis(new_crystal->basis(), ambientMat));

        auto mc = std::make_unique<ReMesocrystal>(std::optional<size_t>{}, *new_crystal->lattice(),
                                                  *new_basis, *new_shape, meso_options,
                                                  new_crystal->position_variance());

        if (!meso_options.use_reciprocal_sum)
            mc->setBasisIndexes(p->calcBasisIndexes());

        std::vector<double> weights;
        std::vector<Material> materials;
        double volume = 0;
        for (const auto* pp : crystal->basis()->decompose()) {
            volume += pp->volume();
            weights.push_back(pp->volume());
            materials.push_back(pp->avgeMaterial());
        }
        mc->setAdmixedFraction(volume / unit_cell_volume);
        mc->setAdmixedMaterial(
            MaterialUtil::averagedMaterial("MesocrystalAvgeMat", weights, materials));

        OwningVector<IReParticle> result;
        result.push_back(mc.release());
        return result;
    }
    ASSERT_NEVER;
}