1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
|
// ************************************************************************************************
//
// BornAgain: simulate and fit reflection and scattering
//
//! @file Sample/HardParticle/Ripples.cpp
//! @brief Implements computations in namespace ripples.
//!
//! @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/Ripples.h"
#include "Base/Math/Functions.h"
#include "Base/Math/IntegratorGK.h"
#include <numbers>
using std::numbers::pi;
complex_t ripples::factor_x_box(complex_t q, double l)
{
return l * Math::sinc(q * l / 2.0);
}
complex_t ripples::factor_x_Gauss(complex_t q, double l)
{
return l * exp(-pow(q * l, 2) / 2.);
}
complex_t ripples::factor_x_Lorentz(complex_t q, double l)
{
return l / (1.0 + pow(q * l, 2));
}
//! Complex form factor of rectangular ripple (bar).
complex_t ripples::profile_yz_bar(complex_t qy, complex_t qz, double width, double height)
{
const complex_t qyWdiv2 = width * qy / 2.0;
const complex_t qzHdiv2 = height * qz / 2.0;
return height * width * exp_I(qzHdiv2) * Math::sinc(qyWdiv2) * Math::sinc(qzHdiv2);
}
//! Complex form factor of cosine ripple.
complex_t ripples::profile_yz_cosine(complex_t qy, complex_t qz, double width, double height)
{
complex_t factor = width / pi;
// analytical expressions for some particular cases
if (qz == 0.) {
if (qy == 0.)
return factor * (pi / 2) * height;
complex_t aaa = qy * width / ((2 * pi));
complex_t aaa2 = aaa * aaa;
if (aaa2 == 1.)
return factor * (pi / 4) * height;
return factor * (pi / 2) * height * Math::sinc(qy * width * 0.5) / (1.0 - aaa2);
}
// numerical integration otherwise
const complex_t ay = qy * width / (2 * pi);
const complex_t az = qz * (height / 2);
const auto integrand = [&](double u) -> complex_t {
return sin(u) * exp_I(az * std::cos(u)) * (ay == 0. ? u : sin(ay * u) / ay);
};
complex_t integral = ComplexIntegrator().integrate(integrand, 0, pi);
return factor * integral * exp_I(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) * Math::sinc(gamma_p) - exp_I(gamma_p) * Math::sinc(gamma_m);
result = mul_I(exp_I(-qyd) * result / (qyW2 * 2.));
}
return factor * result;
}
|
