InterferenceHardDisk.cpp

Go to the documentation of this file.

//  ************************************************************************************************
//
//  BornAgain: simulate and fit reflection and scattering
//
//! @file      Sample/Aggregate/InterferenceHardDisk.cpp
//! @brief     Implements class InterferenceHardDisk.
//!
//! @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/InterferenceHardDisk.h"
#include "Base/Math/Bessel.h"
#include "Base/Math/IntegratorGK.h"
#include "Fit/Param/RealLimits.h"
#include <cmath>
 
namespace {
 
const double p = 7.0 / 3.0 - 4.0 * std::sqrt(3.0) / M_PI;
 
double Czero(double packing)
{
    double numerator = 1.0 + packing + 3.0 * p * packing * packing - p * std::pow(packing, 3);
    double denominator = std::pow(1.0 - packing, 3);
    return -numerator / denominator;
}
 
double S2(double packing)
{
    double factor = 3.0 * packing * packing / 8.0;
    double numerator = 8.0 * (1.0 - 2.0 * p) + (25.0 - 9.0 * p) * p * packing
                       - (7.0 - 3.0 * p) * p * packing * packing;
    double denominator = 1.0 + packing + 3.0 * p * packing * packing - p * std::pow(packing, 3);
    return factor * numerator / denominator;
}
 
double W2(double x)
{
    return 2.0 * (std::acos(x) - x * std::sqrt(1.0 - x * x)) / M_PI;
}
 
} // namespace
 
 
InterferenceHardDisk::InterferenceHardDisk(double radius, double density, double position_var)
    : IInterference(position_var)
    , m_radius(radius)
    , m_density(density)
{
    if (m_radius < 0.0 || m_density < 0.0 || packingRatio() > 0.65)
        throw std::runtime_error("InterferenceHardDisk::validateParameters: "
                                 "radius and density must be positive and packing ratio between "
                                 "0 and 0.65");
    RealLimits::nonnegative().check("Radius", m_radius);
    RealLimits::nonnegative().check("TotalParticleDensity", m_density);
}
 
InterferenceHardDisk* InterferenceHardDisk::clone() const
{
    auto* result = new InterferenceHardDisk(m_radius, m_density, m_position_var);
    return result;
}
 
double InterferenceHardDisk::particleDensity() const
{
    return m_density;
}
 
double InterferenceHardDisk::radius() const
{
    return m_radius;
}
 
double InterferenceHardDisk::density() const
{
    return m_density;
}
 
double InterferenceHardDisk::iff_without_dw(const R3 q) const
{
    const double qx = q.x();
    const double qy = q.y();
    const double q2r = 2.0 * std::sqrt(qx * qx + qy * qy) * m_radius;
    const double packing = packingRatio();
    const double c_zero = Czero(packing);
    const double s2 = S2(packing);
    const double c_q = 2.0 * M_PI
                       * RealIntegrator().integrate(
                           [=](double x) -> double {
                               double cx =
                                   c_zero * (1.0 + 4.0 * packing * (W2(x / 2.0) - 1.0) + s2 * x);
                               return x * cx * Math::Bessel::J0(q2r * x);
                           },
                           0.0, 1.0);
    const double rho = 4.0 * packing / M_PI;
    return 1.0 / (1.0 - rho * c_q);
}
 
double InterferenceHardDisk::packingRatio() const
{
    return M_PI * m_radius * m_radius * m_density;
}