From 1cd8e2c66fc6c5fb9a5bdaee8e59e0a172889343 Mon Sep 17 00:00:00 2001
From: "Joachim Wuttke (h)" <j.wuttke@fz-juelich.de>
Date: Fri, 14 Aug 2020 10:22:55 +0200
Subject: [PATCH] Ripples: move profile ff computation to namespace ripples
---
Core/Correlations/FTDistributions1D.cpp | 2 +-
Core/HardParticle/ProfileRipple1.cpp | 25 +----
Core/HardParticle/ProfileRipple1.h | 1 -
Core/HardParticle/ProfileRipple2.cpp | 28 +----
Core/HardParticle/ProfileRipple2.h | 2 +-
Core/HardParticle/Ripples.cpp | 61 +++++++++++
Core/HardParticle/Ripples.h | 4 +
Core/StandardSamples/TwoDimLatticeBuilder.cpp | 3 +-
auto/Wrap/libBornAgainCore.py | 8 ++
auto/Wrap/libBornAgainCore_wrap.cpp | 102 ++++++++++++++++++
10 files changed, 182 insertions(+), 54 deletions(-)
diff --git a/Core/Correlations/FTDistributions1D.cpp b/Core/Correlations/FTDistributions1D.cpp
index d76afd283b3..bfc4d918a58 100644
--- a/Core/Correlations/FTDistributions1D.cpp
+++ b/Core/Correlations/FTDistributions1D.cpp
@@ -29,7 +29,7 @@ const double CosineDistributionFactor = 1.0 / 3.0 - 2.0 / M_PI / M_PI;
IFTDistribution1D::IFTDistribution1D(const NodeMeta& meta, const std::vector<double>& PValues)
: INode(nodeMetaUnion({{"Omega", "nm", "Half-width", 0, INF, 1.}}, meta), PValues),
- m_omega(m_P[0])
+ m_omega(m_P[0])
{
}
diff --git a/Core/HardParticle/ProfileRipple1.cpp b/Core/HardParticle/ProfileRipple1.cpp
index c2bf9699a23..6255ac8d7ad 100644
--- a/Core/HardParticle/ProfileRipple1.cpp
+++ b/Core/HardParticle/ProfileRipple1.cpp
@@ -14,8 +14,8 @@
#include "Core/HardParticle/ProfileRipple1.h"
#include "Core/Basics/Exceptions.h"
-#include "Core/Basics/MathConstants.h"
#include "Core/Parametrization/RealParameter.h"
+#include "Core/HardParticle/Ripples.h"
#include "Core/Shapes/RippleCosine.h" // from Shapes/
#include "Core/Tools/MathFunctions.h"
#include "Fit/Tools/RealLimits.h"
@@ -41,28 +41,7 @@ double ProfileRipple1::radialExtension() const
//! Complex form factor.
complex_t ProfileRipple1::factor_yz(complex_t qy, complex_t qz) const
{
- complex_t factor = m_width / M_PI;
-
- // analytical expressions for some particular cases
- if (qz == 0.) {
- if (qy == 0.)
- return factor * M_PI_2 * m_height;
- complex_t aaa = qy * m_width / (M_TWOPI);
- complex_t aaa2 = aaa * aaa;
- if (aaa2 == 1.)
- return factor * M_PI_4 * m_height;
- return factor * M_PI_2 * m_height * MathFunctions::sinc(qy * m_width * 0.5) / (1.0 - aaa2);
- }
-
- // numerical integration otherwise
- const complex_t ay = qy * m_width / M_TWOPI;
- const complex_t az = complex_t(0, 1) * qz * (m_height / 2);
-
- const auto integrand = [&](double u) -> complex_t {
- return sin(u) * exp(az * std::cos(u)) * (ay == 0. ? u : sin(ay * u) / ay);
- };
- complex_t integral = ComplexIntegrator().integrate(integrand, 0, M_PI);
- return factor * integral * exp(az) * (m_height / 2);
+ return ripples::profile_yz_cosine(qy, qz, m_width, m_height);
}
void ProfileRipple1::onChange()
diff --git a/Core/HardParticle/ProfileRipple1.h b/Core/HardParticle/ProfileRipple1.h
index c0155a08a0d..e4b00dbdc19 100644
--- a/Core/HardParticle/ProfileRipple1.h
+++ b/Core/HardParticle/ProfileRipple1.h
@@ -16,7 +16,6 @@
#define BORNAGAIN_CORE_HARDPARTICLE_PROFILERIPPLE1_H
#include "Core/Scattering/IFormFactorBorn.h"
-#include "Core/Tools/Integrator.h"
//! Base class for form factors with a cosine ripple profile in the yz plane.
diff --git a/Core/HardParticle/ProfileRipple2.cpp b/Core/HardParticle/ProfileRipple2.cpp
index 83c05077f8e..123dc6afa88 100644
--- a/Core/HardParticle/ProfileRipple2.cpp
+++ b/Core/HardParticle/ProfileRipple2.cpp
@@ -15,6 +15,7 @@
#include "Core/HardParticle/ProfileRipple2.h"
#include "Core/Basics/Exceptions.h"
#include "Core/Basics/MathConstants.h"
+#include "Core/HardParticle/Ripples.h"
#include "Core/Parametrization/RealParameter.h"
#include "Core/Shapes/RippleSawtooth.h" // from Shapes/
#include "Core/Tools/MathFunctions.h"
@@ -42,32 +43,7 @@ double ProfileRipple2::radialExtension() const
//! Complex form factor.
complex_t ProfileRipple2::factor_yz(complex_t qy, complex_t qz) const
{
- complex_t result;
- const complex_t factor = m_height * m_width;
- const complex_t qyW2 = qy * m_width * 0.5;
- const complex_t qyd = qy * m_asymmetry;
- const complex_t qzH = qz * m_height;
- const complex_t a = qzH + qyd;
- // dimensionless scale factors
- const double a_scale = std::abs(a);
- const double w_scale = std::abs(qyW2);
-
- if (w_scale < 1.e-5) { // |q_y*W| << 1
- if (a_scale < 1e-5) { // |q_y*W| << 1 && |q_z*H + q_y*d| << 1
- // relative error is O((q_y*W)^2) and O((q_z*H + q_y*d)^2)
- result = exp_I(-qyd) * (0.5 + mul_I(a) / 6.);
- } else {
- // relative error is O((q_y*W)^2)
- result = exp_I(-qyd) * (1.0 + mul_I(a) - exp_I(a)) / (a * a);
- }
- } else {
- const complex_t gamma_p = (a + qyW2) * 0.5;
- const complex_t gamma_m = (a - qyW2) * 0.5;
- result = exp_I(gamma_m) * MathFunctions::sinc(gamma_p)
- - exp_I(gamma_p) * MathFunctions::sinc(gamma_m);
- result = mul_I(exp_I(-qyd) * result / (qyW2 * 2.));
- }
- return factor * result;
+ return ripples::profile_yz_triangular(qy, qz, m_width, m_height, m_asymmetry);
}
void ProfileRipple2::onChange()
diff --git a/Core/HardParticle/ProfileRipple2.h b/Core/HardParticle/ProfileRipple2.h
index 718c78dc849..069a9954b0a 100644
--- a/Core/HardParticle/ProfileRipple2.h
+++ b/Core/HardParticle/ProfileRipple2.h
@@ -18,7 +18,7 @@
#include "Core/Scattering/IFormFactorBorn.h"
#include "Core/Tools/Integrator.h"
-//! Base class for form factors with a cosine ripple profile in the yz plane.
+//! Base class for form factors with a triangular ripple profile in the yz plane.
class BA_CORE_API_ ProfileRipple2 : public IFormFactorBorn
{
diff --git a/Core/HardParticle/Ripples.cpp b/Core/HardParticle/Ripples.cpp
index a2d88df4678..cd7dee8f5ab 100644
--- a/Core/HardParticle/Ripples.cpp
+++ b/Core/HardParticle/Ripples.cpp
@@ -13,6 +13,8 @@
// ************************************************************************** //
#include "Core/HardParticle/Ripples.h"
+#include "Core/Tools/Integrator.h"
+#include "Core/Basics/MathConstants.h"
#include "Core/Tools/MathFunctions.h"
complex_t ripples::factor_x_box(complex_t q, double r)
@@ -29,3 +31,62 @@ complex_t ripples::factor_x_Lorentz(complex_t q, double r)
{
return r / (1.0 + (q * r) * (q * r));
}
+
+//! Complex form factor of triangular ripple.
+complex_t ripples::profile_yz_cosine(complex_t qy, complex_t qz, double width, double height)
+{
+ complex_t factor = width / M_PI;
+
+ // analytical expressions for some particular cases
+ if (qz == 0.) {
+ if (qy == 0.)
+ return factor * M_PI_2 * height;
+ complex_t aaa = qy * width / (M_TWOPI);
+ complex_t aaa2 = aaa * aaa;
+ if (aaa2 == 1.)
+ return factor * M_PI_4 * height;
+ return factor * M_PI_2 * height * MathFunctions::sinc(qy * width * 0.5) / (1.0 - aaa2);
+ }
+
+ // numerical integration otherwise
+ const complex_t ay = qy * width / M_TWOPI;
+ const complex_t az = complex_t(0, 1) * qz * (height / 2);
+
+ const auto integrand = [&](double u) -> complex_t {
+ return sin(u) * exp(az * std::cos(u)) * (ay == 0. ? u : sin(ay * u) / ay);
+ };
+ complex_t integral = ComplexIntegrator().integrate(integrand, 0, M_PI);
+ return factor * integral * exp(az) * (height / 2);
+}
+
+//! Complex form factor of triangular ripple.
+complex_t ripples::profile_yz_triangular(complex_t qy, complex_t qz, double width, double height,
+ double asymmetry)
+{
+ complex_t result;
+ const complex_t factor = height * width;
+ const complex_t qyW2 = qy * width * 0.5;
+ const complex_t qyd = qy * asymmetry;
+ const complex_t qzH = qz * height;
+ const complex_t a = qzH + qyd;
+ // dimensionless scale factors
+ const double a_scale = std::abs(a);
+ const double w_scale = std::abs(qyW2);
+
+ if (w_scale < 1.e-5) { // |q_y*W| << 1
+ if (a_scale < 1e-5) { // |q_y*W| << 1 && |q_z*H + q_y*d| << 1
+ // relative error is O((q_y*W)^2) and O((q_z*H + q_y*d)^2)
+ result = exp_I(-qyd) * (0.5 + mul_I(a) / 6.);
+ } else {
+ // relative error is O((q_y*W)^2)
+ result = exp_I(-qyd) * (1.0 + mul_I(a) - exp_I(a)) / (a * a);
+ }
+ } else {
+ const complex_t gamma_p = (a + qyW2) * 0.5;
+ const complex_t gamma_m = (a - qyW2) * 0.5;
+ result = exp_I(gamma_m) * MathFunctions::sinc(gamma_p)
+ - exp_I(gamma_p) * MathFunctions::sinc(gamma_m);
+ result = mul_I(exp_I(-qyd) * result / (qyW2 * 2.));
+ }
+ return factor * result;
+}
diff --git a/Core/HardParticle/Ripples.h b/Core/HardParticle/Ripples.h
index 92f37c36fbc..e700b70800c 100644
--- a/Core/HardParticle/Ripples.h
+++ b/Core/HardParticle/Ripples.h
@@ -25,6 +25,10 @@ complex_t factor_x_box(complex_t q, double l);
complex_t factor_x_Gauss(complex_t q, double l);
complex_t factor_x_Lorentz(complex_t q, double l);
+complex_t profile_yz_cosine(complex_t qy, complex_t qz, double width, double height);
+complex_t profile_yz_triangular(complex_t qy, complex_t qz, double width, double height,
+ double asymmetry);
+
} // namespace ripples
#endif // BORNAGAIN_CORE_HARDPARTICLE_RIPPLES_H
diff --git a/Core/StandardSamples/TwoDimLatticeBuilder.cpp b/Core/StandardSamples/TwoDimLatticeBuilder.cpp
index 06e645d858e..ec77b444569 100644
--- a/Core/StandardSamples/TwoDimLatticeBuilder.cpp
+++ b/Core/StandardSamples/TwoDimLatticeBuilder.cpp
@@ -155,8 +155,7 @@ MultiLayer* RotatedSquareLatticeBuilder::buildSample() const
std::unique_ptr<InterferenceFunction2DLattice> P_interference_function{
InterferenceFunction2DLattice::createSquare(10.0 * Units::nanometer, 30.0 * Units::degree)};
FTDecayFunction2DCauchy pdf(300.0 * Units::nanometer / 2.0 / M_PI,
- 100.0 * Units::nanometer / 2.0 / M_PI,
- 30.0 * Units::degree);
+ 100.0 * Units::nanometer / 2.0 / M_PI, 30.0 * Units::degree);
P_interference_function->setDecayFunction(pdf);
ParticleLayout particle_layout;
diff --git a/auto/Wrap/libBornAgainCore.py b/auto/Wrap/libBornAgainCore.py
index 46912fc9f6e..761dac428f7 100644
--- a/auto/Wrap/libBornAgainCore.py
+++ b/auto/Wrap/libBornAgainCore.py
@@ -17715,6 +17715,14 @@ def factor_x_Lorentz(q, l):
"""
return _libBornAgainCore.factor_x_Lorentz(q, l)
+
+def profile_yz_cosine(qy, qz, width, height):
+ r"""profile_yz_cosine(complex_t qy, complex_t qz, double width, double height) -> complex_t"""
+ return _libBornAgainCore.profile_yz_cosine(qy, qz, width, height)
+
+def profile_yz_triangular(qy, qz, width, height, asymmetry):
+ r"""profile_yz_triangular(complex_t qy, complex_t qz, double width, double height, double asymmetry) -> complex_t"""
+ return _libBornAgainCore.profile_yz_triangular(qy, qz, width, height, asymmetry)
class FormFactorGaussSphere(IFormFactorBorn):
r"""Proxy of C++ FormFactorGaussSphere class."""
diff --git a/auto/Wrap/libBornAgainCore_wrap.cpp b/auto/Wrap/libBornAgainCore_wrap.cpp
index 9a12df0aca5..d2ae11d1204 100644
--- a/auto/Wrap/libBornAgainCore_wrap.cpp
+++ b/auto/Wrap/libBornAgainCore_wrap.cpp
@@ -99998,6 +99998,106 @@ fail:
}
+SWIGINTERN PyObject *_wrap_profile_yz_cosine(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+ PyObject *resultobj = 0;
+ complex_t arg1 ;
+ complex_t arg2 ;
+ double arg3 ;
+ double arg4 ;
+ std::complex< double > val1 ;
+ int ecode1 = 0 ;
+ std::complex< double > val2 ;
+ int ecode2 = 0 ;
+ double val3 ;
+ int ecode3 = 0 ;
+ double val4 ;
+ int ecode4 = 0 ;
+ PyObject *swig_obj[4] ;
+ complex_t result;
+
+ if (!SWIG_Python_UnpackTuple(args, "profile_yz_cosine", 4, 4, swig_obj)) SWIG_fail;
+ ecode1 = SWIG_AsVal_std_complex_Sl_double_Sg_(swig_obj[0], &val1);
+ if (!SWIG_IsOK(ecode1)) {
+ SWIG_exception_fail(SWIG_ArgError(ecode1), "in method '" "profile_yz_cosine" "', argument " "1"" of type '" "complex_t""'");
+ }
+ arg1 = static_cast< complex_t >(val1);
+ ecode2 = SWIG_AsVal_std_complex_Sl_double_Sg_(swig_obj[1], &val2);
+ if (!SWIG_IsOK(ecode2)) {
+ SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "profile_yz_cosine" "', argument " "2"" of type '" "complex_t""'");
+ }
+ arg2 = static_cast< complex_t >(val2);
+ ecode3 = SWIG_AsVal_double(swig_obj[2], &val3);
+ if (!SWIG_IsOK(ecode3)) {
+ SWIG_exception_fail(SWIG_ArgError(ecode3), "in method '" "profile_yz_cosine" "', argument " "3"" of type '" "double""'");
+ }
+ arg3 = static_cast< double >(val3);
+ ecode4 = SWIG_AsVal_double(swig_obj[3], &val4);
+ if (!SWIG_IsOK(ecode4)) {
+ SWIG_exception_fail(SWIG_ArgError(ecode4), "in method '" "profile_yz_cosine" "', argument " "4"" of type '" "double""'");
+ }
+ arg4 = static_cast< double >(val4);
+ result = ripples::profile_yz_cosine(arg1,arg2,arg3,arg4);
+ resultobj = SWIG_From_std_complex_Sl_double_Sg_(static_cast< std::complex<double> >(result));
+ return resultobj;
+fail:
+ return NULL;
+}
+
+
+SWIGINTERN PyObject *_wrap_profile_yz_triangular(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+ PyObject *resultobj = 0;
+ complex_t arg1 ;
+ complex_t arg2 ;
+ double arg3 ;
+ double arg4 ;
+ double arg5 ;
+ std::complex< double > val1 ;
+ int ecode1 = 0 ;
+ std::complex< double > val2 ;
+ int ecode2 = 0 ;
+ double val3 ;
+ int ecode3 = 0 ;
+ double val4 ;
+ int ecode4 = 0 ;
+ double val5 ;
+ int ecode5 = 0 ;
+ PyObject *swig_obj[5] ;
+ complex_t result;
+
+ if (!SWIG_Python_UnpackTuple(args, "profile_yz_triangular", 5, 5, swig_obj)) SWIG_fail;
+ ecode1 = SWIG_AsVal_std_complex_Sl_double_Sg_(swig_obj[0], &val1);
+ if (!SWIG_IsOK(ecode1)) {
+ SWIG_exception_fail(SWIG_ArgError(ecode1), "in method '" "profile_yz_triangular" "', argument " "1"" of type '" "complex_t""'");
+ }
+ arg1 = static_cast< complex_t >(val1);
+ ecode2 = SWIG_AsVal_std_complex_Sl_double_Sg_(swig_obj[1], &val2);
+ if (!SWIG_IsOK(ecode2)) {
+ SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "profile_yz_triangular" "', argument " "2"" of type '" "complex_t""'");
+ }
+ arg2 = static_cast< complex_t >(val2);
+ ecode3 = SWIG_AsVal_double(swig_obj[2], &val3);
+ if (!SWIG_IsOK(ecode3)) {
+ SWIG_exception_fail(SWIG_ArgError(ecode3), "in method '" "profile_yz_triangular" "', argument " "3"" of type '" "double""'");
+ }
+ arg3 = static_cast< double >(val3);
+ ecode4 = SWIG_AsVal_double(swig_obj[3], &val4);
+ if (!SWIG_IsOK(ecode4)) {
+ SWIG_exception_fail(SWIG_ArgError(ecode4), "in method '" "profile_yz_triangular" "', argument " "4"" of type '" "double""'");
+ }
+ arg4 = static_cast< double >(val4);
+ ecode5 = SWIG_AsVal_double(swig_obj[4], &val5);
+ if (!SWIG_IsOK(ecode5)) {
+ SWIG_exception_fail(SWIG_ArgError(ecode5), "in method '" "profile_yz_triangular" "', argument " "5"" of type '" "double""'");
+ }
+ arg5 = static_cast< double >(val5);
+ result = ripples::profile_yz_triangular(arg1,arg2,arg3,arg4,arg5);
+ resultobj = SWIG_From_std_complex_Sl_double_Sg_(static_cast< std::complex<double> >(result));
+ return resultobj;
+fail:
+ return NULL;
+}
+
+
SWIGINTERN PyObject *_wrap_new_FormFactorGaussSphere__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
PyObject *resultobj = 0;
std::vector< double,std::allocator< double > > arg1 ;
@@ -131691,6 +131791,8 @@ static PyMethodDef SwigMethods[] = {
"complex_t ripples::factor_x_Lorentz(complex_t q, double l)\n"
"\n"
""},
+ { "profile_yz_cosine", _wrap_profile_yz_cosine, METH_VARARGS, "profile_yz_cosine(complex_t qy, complex_t qz, double width, double height) -> complex_t"},
+ { "profile_yz_triangular", _wrap_profile_yz_triangular, METH_VARARGS, "profile_yz_triangular(complex_t qy, complex_t qz, double width, double height, double asymmetry) -> complex_t"},
{ "new_FormFactorGaussSphere", _wrap_new_FormFactorGaussSphere, METH_VARARGS, "\n"
"FormFactorGaussSphere(vdouble1d_t P)\n"
"new_FormFactorGaussSphere(double mean_radius) -> FormFactorGaussSphere\n"
--
GitLab