//  ************************************************************************************************
//
//  BornAgain: simulate and fit reflection and scattering
//
//! @file      Sample/Correlation/Profiles1D.cpp
//! @brief     Implements class interface IProfile2D and children thereof.
//!
//! @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/Correlation/Profiles1D.h"
#include "Base/Math/Functions.h"
#include "Base/Py/PyFmt.h"
#include "Base/Util/Assert.h"
#include <limits>
#include <numbers>

using std::numbers::pi;

namespace {

const double CosineDistributionFactor = 1.0 / 3.0 - 2.0 / pi / pi;

} // namespace


//  ************************************************************************************************
//  interface IProfile1D
//  ************************************************************************************************

IProfile1D::IProfile1D(const std::vector<double>& PValues)
    : INode(PValues)
    , m_omega(m_P[0])
{
}

std::string IProfile1D::pythonConstructor() const
{
    ASSERT(m_P.size() == 1);
    return Py::Fmt::printFunction(className(), m_omega, "nm");
}

std::string IProfile1D::validate() const
{
    std::vector<std::string> errs;
    requestGt0(errs, m_omega, "omega");
    if (!errs.empty())
        return jointError(errs);
    m_validated = true;
    return "";
}

//  ************************************************************************************************
//  class Profile1DCauchy
//  ************************************************************************************************

Profile1DCauchy::Profile1DCauchy(const std::vector<double> P)
    : IProfile1D(P)
{
    validateOrThrow();
}

Profile1DCauchy::Profile1DCauchy(double omega)
    : Profile1DCauchy(std::vector<double>{omega})
{
}

Profile1DCauchy* Profile1DCauchy::clone() const
{
    return new Profile1DCauchy(m_omega);
}

double Profile1DCauchy::standardizedFT(double q) const
{
    double sum_sq = q * q * m_omega * m_omega;
    return 1.0 / (1.0 + sum_sq);
}

double Profile1DCauchy::decayFT(double q) const
{
    double sum_sq = q * q * m_omega * m_omega;
    return m_omega * 2.0 / (1.0 + sum_sq);
}

double Profile1DCauchy::qSecondDerivative() const
{
    return 2.0 * m_omega * m_omega;
}

std::unique_ptr<IDistribution1DSampler> Profile1DCauchy::createSampler() const
{
    return std::make_unique<Distribution1DCauchySampler>(1 / m_omega);
}

//  ************************************************************************************************
//  class Profile1DGauss
//  ************************************************************************************************

Profile1DGauss::Profile1DGauss(const std::vector<double> P)
    : IProfile1D(P)
{
    validateOrThrow();
}

Profile1DGauss::Profile1DGauss(double omega)
    : Profile1DGauss(std::vector<double>{omega})
{
}

Profile1DGauss* Profile1DGauss::clone() const
{
    return new Profile1DGauss(m_omega);
}

double Profile1DGauss::standardizedFT(double q) const
{
    double sum_sq = q * q * m_omega * m_omega;
    return std::exp(-sum_sq / 2.0);
}

double Profile1DGauss::decayFT(double q) const
{
    double sum_sq = q * q * m_omega * m_omega;
    return m_omega * std::sqrt((2 * pi)) * std::exp(-sum_sq / 2.0);
}

double Profile1DGauss::qSecondDerivative() const
{
    return m_omega * m_omega;
}

std::unique_ptr<IDistribution1DSampler> Profile1DGauss::createSampler() const
{
    return std::make_unique<Distribution1DGaussSampler>(0.0, m_omega);
}

//  ************************************************************************************************
//  class Profile1DGate
//  ************************************************************************************************

Profile1DGate::Profile1DGate(const std::vector<double> P)
    : IProfile1D(P)
{
    validateOrThrow();
}

Profile1DGate::Profile1DGate(double omega)
    : Profile1DGate(std::vector<double>{omega})
{
}

Profile1DGate* Profile1DGate::clone() const
{
    return new Profile1DGate(m_omega);
}

double Profile1DGate::standardizedFT(double q) const
{
    return Math::sinc(q * m_omega);
}

double Profile1DGate::decayFT(double) const
{
    ASSERT_NEVER; // not yet needed, not yet implemented
}

double Profile1DGate::qSecondDerivative() const
{
    return m_omega * m_omega / 3.0;
}

std::unique_ptr<IDistribution1DSampler> Profile1DGate::createSampler() const
{
    return std::make_unique<Distribution1DGateSampler>(-m_omega, m_omega);
}

//  ************************************************************************************************
//  class Profile1DTriangle
//  ************************************************************************************************

Profile1DTriangle::Profile1DTriangle(const std::vector<double> P)
    : IProfile1D(P)
{
    validateOrThrow();
}

Profile1DTriangle::Profile1DTriangle(double omega)
    : Profile1DTriangle(std::vector<double>{omega})
{
}

Profile1DTriangle* Profile1DTriangle::clone() const
{
    return new Profile1DTriangle(m_omega);
}

double Profile1DTriangle::standardizedFT(double q) const
{
    double sincqw2 = Math::sinc(q * m_omega / 2.0);
    return sincqw2 * sincqw2;
}

double Profile1DTriangle::decayFT(double q) const
{
    double sincqw2 = Math::sinc(q * m_omega / 2.0);
    return m_omega * sincqw2 * sincqw2;
}

double Profile1DTriangle::qSecondDerivative() const
{
    return m_omega * m_omega / 6.0;
}

std::unique_ptr<IDistribution1DSampler> Profile1DTriangle::createSampler() const
{
    return std::make_unique<Distribution1DTriangleSampler>(m_omega);
}

//  ************************************************************************************************
//  class Profile1DCosine
//  ************************************************************************************************

Profile1DCosine::Profile1DCosine(const std::vector<double> P)
    : IProfile1D(P)
{
    validateOrThrow();
}

Profile1DCosine::Profile1DCosine(double omega)
    : Profile1DCosine(std::vector<double>{omega})
{
}

Profile1DCosine* Profile1DCosine::clone() const
{
    return new Profile1DCosine(m_omega);
}

double Profile1DCosine::standardizedFT(double q) const
{
    double qw = std::abs(q * m_omega);
    if (std::abs(1.0 - qw * qw / pi / pi) < std::numeric_limits<double>::epsilon())
        return 0.5;
    return Math::sinc(qw) / (1.0 - qw * qw / pi / pi);
}

double Profile1DCosine::decayFT(double) const
{
    ASSERT_NEVER; // not yet needed, not yet implemented
}

double Profile1DCosine::qSecondDerivative() const
{
    return CosineDistributionFactor * m_omega * m_omega;
}

std::unique_ptr<IDistribution1DSampler> Profile1DCosine::createSampler() const
{
    return std::make_unique<Distribution1DCosineSampler>(m_omega);
}

//  ************************************************************************************************
//  class Profile1DVoigt
//  ************************************************************************************************

Profile1DVoigt::Profile1DVoigt(const std::vector<double> P)
    : IProfile1D(P)
    , m_eta(m_P[1])
{
    validateOrThrow();
}

Profile1DVoigt::Profile1DVoigt(double omega, double eta)
    : Profile1DVoigt(std::vector<double>{omega, eta})
{
}

Profile1DVoigt* Profile1DVoigt::clone() const
{
    return new Profile1DVoigt(m_omega, m_eta);
}

double Profile1DVoigt::standardizedFT(double q) const
{
    double sum_sq = q * q * m_omega * m_omega;
    return m_eta * std::exp(-sum_sq / 2.0) + (1.0 - m_eta) * 1.0 / (1.0 + sum_sq);
}

double Profile1DVoigt::decayFT(double q) const
{
    double sum_sq = q * q * m_omega * m_omega;
    return m_eta * m_omega * std::sqrt((2 * pi)) * std::exp(-sum_sq / 2.0)
           + (1.0 - m_eta) * m_omega * 2.0 / (1.0 + sum_sq);
}

double Profile1DVoigt::qSecondDerivative() const
{
    return (2.0 - m_eta) * m_omega * m_omega;
}

std::unique_ptr<IDistribution1DSampler> Profile1DVoigt::createSampler() const
{
    throw std::runtime_error("Profile1DVoigt::createSampler not yet implemented."
                             " Please contact the maintainers");
}

std::string Profile1DVoigt::pythonConstructor() const
{
    ASSERT(m_P.size() == 2);
    return Py::Fmt::printFunction(className(), m_omega, "nm", m_eta, "");
}

std::string Profile1DVoigt::validate() const
{
    std::vector<std::string> errs;
    requestGt0(errs, m_omega, "omega");
    requestIn(errs, m_eta, "eta", 0, 1);
    if (!errs.empty())
        return jointError(errs);
    m_validated = true;
    return "";
}