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
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
|
// ************************************************************************************************
//
// BornAgain: simulate and fit reflection and scattering
//
//! @file Sample/Aggregate/Interference2DParacrystal.cpp
//! @brief Implements class Interference2DParacrystal.
//!
//! @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/Interference2DParacrystal.h"
#include "Base/Math/IntegratorGK.h"
#include "Base/Util/Assert.h"
#include <limits>
Interference2DParacrystal::Interference2DParacrystal(const Lattice2D& lattice,
double damping_length, double domain_size_1,
double domain_size_2)
: IInterference(0)
, m_integrate_xi(false)
, m_damping_length(damping_length)
, m_domain_sizes({domain_size_1, domain_size_2})
{
m_lattice.reset(lattice.clone());
validateOrThrow();
}
Interference2DParacrystal::~Interference2DParacrystal() = default;
Interference2DParacrystal* Interference2DParacrystal::clone() const
{
auto* result = new Interference2DParacrystal(*m_lattice, m_damping_length, m_domain_sizes[0],
m_domain_sizes[1]);
result->setPositionVariance(m_position_var);
if (m_pdf1 && m_pdf2)
result->setProbabilityDistributions(*m_pdf1, *m_pdf2);
result->setIntegrationOverXi(m_integrate_xi);
return result;
}
//! Sets the probability distributions (Fourier transformed) for the two lattice directions.
//! @param pdf_1: probability distribution in first lattice direction
//! @param pdf_2: probability distribution in second lattice direction
void Interference2DParacrystal::setProbabilityDistributions(const IProfile2D& pdf_1,
const IProfile2D& pdf_2)
{
m_pdf1.reset(pdf_1.clone());
m_pdf2.reset(pdf_2.clone());
}
//! Sets the damping length.
//! @param damping_length: the damping (coherence) length of the paracrystal in nanometers
void Interference2DParacrystal::setDampingLength(double damping_length)
{
m_damping_length = damping_length;
}
double Interference2DParacrystal::particleDensity() const
{
double area = m_lattice->unitCellArea();
return area == 0.0 ? 0.0 : 1.0 / area;
}
std::vector<const INode*> Interference2DParacrystal::nodeChildren() const
{
return std::vector<const INode*>() << m_pdf1 << m_pdf2 << m_lattice;
}
double Interference2DParacrystal::iff_without_dw(const R3& q) const
{
double range = pi;
if (!m_integrate_xi)
return interferenceForXi(m_lattice->rotationAngle(), q.x(), q.y());
return RealIntegrator().integrate(
[&](double xi) -> double { return interferenceForXi(xi, q.x(), q.y()); }, 0.0, range)
/ range;
}
//! Sets the sizes of coherence domains.
//! @param size_1: coherence domain size along the first basis vector in nanometers
//! @param size_2: coherence domain size along the second basis vector in nanometers
void Interference2DParacrystal::setDomainSizes(double size_1, double size_2)
{
m_domain_sizes = {size_1, size_2};
}
void Interference2DParacrystal::transformToPrincipalAxes(double qx, double qy, double gamma,
double delta, double& q_pa_1,
double& q_pa_2) const
{
q_pa_1 = qx * std::cos(gamma) + qy * std::sin(gamma);
q_pa_2 = qx * std::cos(gamma + delta) + qy * std::sin(gamma + delta);
}
//! Returns interference function for fixed angle xi.
double Interference2DParacrystal::interferenceForXi(double xi, double qx, double qy) const
{
// don't touch order of computation; problems under Windows
double rx = interference1D(qx, qy, xi, 0);
double ry = interference1D(qx, qy, xi + m_lattice->latticeAngle(), 1);
return rx * ry;
}
//! Returns interference function for fixed xi in the dimension determined by the given index.
double Interference2DParacrystal::interference1D(double qx, double qy, double xi,
size_t index) const
{
ASSERT(m_validated);
ASSERT(index <= 1);
ASSERT(m_pdf1 && m_pdf2);
double length = index ? m_lattice->length2() : m_lattice->length1();
int n = static_cast<int>(std::abs(m_domain_sizes[index] / length));
auto nd = static_cast<double>(n);
complex_t fp = FTPDF(qx, qy, xi, index);
if (n < 1)
return ((1.0 + fp) / (1.0 - fp)).real();
if (std::norm(1.0 - fp) < std::numeric_limits<double>::epsilon())
return nd;
// for (1-fp)*nd small, take the series expansion to second order in nd*(1-fp)
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;
return 1.0 + 2.0 * intermediate.real();
}
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);
return 1.0 + 2.0 * intermediate.real();
}
complex_t Interference2DParacrystal::FTPDF(double qx, double qy, double xi, size_t index) const
{
ASSERT(m_validated);
double length = (index ? m_lattice->length2() : m_lattice->length1());
const IProfile2D* pdf = (index ? m_pdf2.get() : m_pdf1.get());
double qa = qx * length * std::cos(xi) + qy * length * std::sin(xi);
complex_t phase = exp_I(qa);
// transform q to principal axes:
double qp1, qp2;
double gamma = xi + pdf->gamma();
double delta = pdf->delta();
transformToPrincipalAxes(qx, qy, gamma, delta, qp1, qp2);
double amplitude = pdf->standardizedFT2D(qp1, qp2);
complex_t result = phase * amplitude;
if (m_damping_length != 0.0)
result *= std::exp(-length / m_damping_length);
return result;
}
std::vector<double> Interference2DParacrystal::domainSizes() const
{
return {m_domain_sizes[0], m_domain_sizes[1]};
}
//! Enables/disables averaging over the lattice rotation angle.
//! @param integrate_xi: integration flag
void Interference2DParacrystal::setIntegrationOverXi(bool integrate_xi)
{
m_integrate_xi = integrate_xi;
m_lattice->setRotationEnabled(!m_integrate_xi); // deregister Xi in the case of integration
}
const Lattice2D& Interference2DParacrystal::lattice() const
{
ASSERT(m_lattice);
return *m_lattice;
}
std::string Interference2DParacrystal::validate() const
{
std::vector<std::string> errs;
requestGe0(errs, m_damping_length, "DampingLength");
requestGe0(errs, m_domain_sizes[0], "DomainSize1");
requestGe0(errs, m_domain_sizes[1], "DomainSize2");
if (!errs.empty())
return jointError(errs);
m_validated = true;
return "";
}
|
