InterferenceRadialParaCrystal.cpp

Go to the documentation of this file.

//  ************************************************************************************************
//
//  BornAgain: simulate and fit reflection and scattering
//
//! @file      Sample/Aggregate/InterferenceRadialParaCrystal.cpp
//! @brief     Implements class InterferenceRadialParaCrystal.
//!
//! @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/Aggregate/InterferenceRadialParaCrystal.h"
#include "Fit/Param/RealLimits.h"
#include <limits>
 
//! Constructor of interference function of radial paracrystal.
//! @param peak_distance: average distance to the next neighbor in nanometers
//! @param damping_length: the damping (coherence) length of the paracrystal in nanometers
InterferenceRadialParaCrystal::InterferenceRadialParaCrystal(double peak_distance,
                                                             double damping_length)
    : IInterference(0)
    , m_peak_distance(peak_distance)
    , m_damping_length(damping_length)
    , m_use_damping_length(true)
    , m_kappa(0.0)
    , m_domain_size(0.0)
{
    if (m_damping_length == 0.0)
        m_use_damping_length = false;
 
    RealLimits::nonnegative().check("PeakDistance", m_peak_distance);
    RealLimits::nonnegative().check("DampingLength", m_damping_length);
    RealLimits::nonnegative().check("SizeSpaceCoupling", m_kappa);
    RealLimits::nonnegative().check("DomainSize", m_domain_size);
}
 
InterferenceRadialParaCrystal* InterferenceRadialParaCrystal::clone() const
{
    auto* result = new InterferenceRadialParaCrystal(m_peak_distance, m_damping_length);
    result->setPositionVariance(m_position_var);
    if (m_pdf)
        result->setProbabilityDistribution(*m_pdf);
    result->setKappa(m_kappa);
    result->setDomainSize(m_domain_size);
    return result;
}
 
//! Sets size spacing coupling parameter of the Size Spacing Correlation Approximation.
void InterferenceRadialParaCrystal::setKappa(double kappa)
{
    m_kappa = kappa;
}
 
double InterferenceRadialParaCrystal::kappa() const
{
    return m_kappa;
}
 
//! Sets domain size (finite size corrections).
//! @param size: size of coherence domain along the lattice main axis in nanometers
 
void InterferenceRadialParaCrystal::setDomainSize(double size)
{
    m_domain_size = size;
}
 
complex_t InterferenceRadialParaCrystal::FTPDF(double qpar) const
{
    complex_t phase = exp_I(qpar * m_peak_distance);
    double amplitude = m_pdf->standardizedFT(qpar);
    complex_t result = phase * amplitude;
    if (m_use_damping_length)
        result *= std::exp(-m_peak_distance / m_damping_length);
    return result;
}
 
//! Sets one-dimensional probability distribution.
//! @param pdf: probability distribution (Fourier transform of probability density)
 
void InterferenceRadialParaCrystal::setProbabilityDistribution(const IProfile1D& pdf)
{
    m_pdf.reset(pdf.clone());
}
 
std::vector<const INode*> InterferenceRadialParaCrystal::nodeChildren() const
{
    return std::vector<const INode*>() << m_pdf;
}
 
double InterferenceRadialParaCrystal::iff_without_dw(const R3 q) const
{
    if (!m_pdf)
        throw std::runtime_error("InterferenceRadialParaCrystal::"
                                 "evaluate() -> Error! Probability distribution for "
                                 "interference function not properly initialized");
    double result = 0.0;
    double qxr = q.x();
    double qyr = q.y();
    double qpar = std::sqrt(qxr * qxr + qyr * qyr);
    int n = static_cast<int>(std::abs(m_domain_size / m_peak_distance));
    auto nd = static_cast<double>(n);
    complex_t fp = FTPDF(qpar);
    if (n < 1) {
        if (std::abs(1.0 - fp) < 10. * std::numeric_limits<double>::epsilon())
            result = m_pdf->qSecondDerivative() / m_peak_distance / m_peak_distance;
        else
            result = ((1.0 + fp) / (1.0 - fp)).real();
    } else {
        if (std::norm(1.0 - fp) < 10. * std::numeric_limits<double>::epsilon())
            result = nd;
        // for (1-fp)*nd small, take the series expansion to second order in nd*(1-fp)
        else if (std::abs(1.0 - fp) * nd < 2e-4) {
            complex_t intermediate =
                (nd - 1.0) / 2.0 + (nd * nd - 1.0) * (fp - 1.0) / 6.0
                + (nd * nd * nd - 2.0 * nd * nd - nd + 2.0) * (fp - 1.0) * (fp - 1.0) / 24.0;
            result = 1.0 + 2.0 * intermediate.real();
        } else {
            complex_t tmp;
            if (std::abs(fp) == 0.0
                || std::log(std::abs(fp)) * nd < std::log(std::numeric_limits<double>::min())) {
                tmp = 0.0;
            } else {
                tmp = std::pow(fp, n);
            }
            complex_t intermediate =
                fp / (1.0 - fp) - fp * (1.0 - tmp) / nd / (1.0 - fp) / (1.0 - fp);
            result = 1.0 + 2.0 * intermediate.real();
        }
    }
    return result;
}