AutocorrelationModels.cpp

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

#include "Sample/Interface/AutocorrelationModels.h"
#include "Base/Py/PyFmt.h"
#include "Base/Util/Assert.h"
#include <gsl/gsl_sf_bessel.h>
#include <numbers>

using std::numbers::pi;

//! @param sigma: rms of the roughness in nanometers
//! @param hurstParameter: hurst parameter which describes how jagged the interface,
//! dimensionless (0.0, 1.0], where 0.0 gives more spikes, 1.0 more smoothness
//! @param lateralCorrLength: lateral correlation length of the roughness in nanometers
K_CorrelationModel::K_CorrelationModel(double sigma, double hurst, double lateralCorrLength)
    : m_sigma(sigma)
    , m_hurst_parameter(hurst)
    , m_lateral_corr_length(lateralCorrLength)
{
    validateOrThrow();
}

K_CorrelationModel* K_CorrelationModel::clone() const
{
    return new K_CorrelationModel(m_sigma, m_hurst_parameter, m_lateral_corr_length);
}

std::string K_CorrelationModel::validate() const
{
    std::vector<std::string> errs;
    requestGe0(errs, m_sigma, "sigma");
    if (m_sigma > 0) {
        requestIn(errs, m_hurst_parameter, "hurst", 0, 1);
        requestGe0(errs, m_lateral_corr_length, "lateralCorrLength"); // may be zero if unused
    }
    if (!errs.empty())
        return jointError(errs);
    m_validated = true;
    return "";
}

std::string K_CorrelationModel::pythonArguments() const
{
    return Py::Fmt::printArguments({{m_sigma, parDefs()[0].unit},
                                    {m_hurst_parameter, parDefs()[1].unit},
                                    {m_lateral_corr_length, parDefs()[2].unit}});
}

//! Power spectral density of the surface roughness is a result of two-dimensional
//! Fourier transform of the correlation function of the roughness profile.
//!
//! Based on Palasantzas, Phys Rev B, 48, 14472 (1993)
double K_CorrelationModel::spectralFunction(const R3& k) const
{
    ASSERT(m_validated);
    double H = m_hurst_parameter;
    double clength2 = m_lateral_corr_length * m_lateral_corr_length;
    double Qpar2 = k.magxy2();
    return 4.0 * pi * H * m_sigma * m_sigma * clength2 * std::pow(1 + Qpar2 * clength2, -1 - H);
}

//! Correlation function of the roughness profile
double K_CorrelationModel::corrFunction(const R3& k) const
{
    ASSERT(m_validated);
    double H = m_hurst_parameter;
    double clength = m_lateral_corr_length;
    double R = k.magxy();
    return m_sigma * m_sigma * std::pow(2., 1 - H) / tgamma(H) * std::pow(R / clength, H)
           * gsl_sf_bessel_Knu(H, R / clength);
}