TruncatedSpheroid.cpp

Go to the documentation of this file.

//  ************************************************************************************************
//
//  BornAgain: simulate and fit reflection and scattering
//
//! @file      Sample/HardParticle/TruncatedSpheroid.cpp
//! @brief     Implements class TruncatedSpheroid.
//!
//! @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/HardParticle/TruncatedSpheroid.h"
#include "Base/Math/Bessel.h"
#include "Base/Math/Constants.h"
#include "Base/Math/IntegratorGK.h"
#include "Sample/Shapes/TruncatedEllipsoidNet.h"
#include <limits>
 
TruncatedSpheroid::TruncatedSpheroid(const std::vector<double> P)
    : IFormFactor(P)
    , m_radius(m_P[0])
    , m_height(m_P[1])
    , m_height_flattening(m_P[2])
    , m_dh(m_P[3])
{
    checkNodeArgs();
    if (m_height > 2. * m_radius * m_height_flattening || m_dh > m_height) {
        std::ostringstream ostr;
        ostr << "::TruncatedSpheroid() -> Error in class initialization with parameters ";
        ostr << " radius:" << m_radius;
        ostr << " height:" << m_height;
        ostr << " height_flattening:" << m_height_flattening << "\n\n";
        ostr << "Check for 'height <= 2.*radius*height_flattening' failed.";
        throw std::runtime_error(ostr.str());
    }
    m_shape3D = std::make_unique<TruncatedEllipsoidNet>(
        m_radius, m_radius, m_height_flattening * m_radius, m_height, m_dh);
}
 
TruncatedSpheroid::TruncatedSpheroid(double radius, double height, double height_flattening,
                                     double dh)
    : TruncatedSpheroid(std::vector<double>{radius, height, height_flattening, dh})
{
}
 
complex_t TruncatedSpheroid::formfactor_at_bottom(C3 q) const
{
    double H = m_height;
    double R = m_radius;
    double fp = m_height_flattening;
 
    if (std::abs(q.mag()) <= std::numeric_limits<double>::epsilon())
        return M_PI / 3. / fp * (H * H * (3. * R - H / fp) - m_dh * m_dh * (3. * R - m_dh / fp));
 
    complex_t z_part = exp_I(q.z() * (H - fp * R));
    return M_TWOPI * z_part
           * ComplexIntegrator().integrate(
               [=](double Z) {
                   double R = m_radius;
                   double fp = m_height_flattening;
 
                   double Rz = std::sqrt(R * R - Z * Z / (fp * fp));
                   complex_t qrRz = std::sqrt(q.x() * q.x() + q.y() * q.y()) * Rz;
                   complex_t J1_qrRz_div_qrRz = Math::Bessel::J1c(qrRz);
 
                   return Rz * Rz * J1_qrRz_div_qrRz * exp_I(q.z() * Z);
               },
               fp * R - H, fp * R - m_dh);
}