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
|
// ************************************************************************************************
//
// BornAgain: simulate and fit reflection and scattering
//
//! @file Resample/Slice/ProfileHelper.cpp
//! @brief Implements class ProfileHelper.
//!
//! @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 "Resample/Slice/ProfileHelper.h"
#include "Base/Math/Functions.h"
#include "Base/Util/Assert.h"
#include "Sample/Interface/Roughness.h"
#include <numbers>
using std::numbers::pi;
namespace {
double Transition(double x, const Roughness* roughness, double rms)
{
if (!roughness)
return Math::Heaviside(x);
return roughness->transient()->transient(x, rms);
}
const std::string SLD = "SLD";
const std::string X = "X";
const std::string Y = "Y";
const std::string Z = "Z";
complex_t quantity(const Material& mat, std::string q)
{
if (q == SLD)
return mat.refractiveIndex_or_SLD();
if (q == X)
return mat.magnetization().x();
if (q == Y)
return mat.magnetization().y();
if (q == Z)
return mat.magnetization().z();
ASSERT_NEVER;
}
std::vector<double> vec2real(const std::vector<complex_t>& v)
{
std::vector<double> result(v.size());
for (size_t i = 0; i < v.size(); i++)
result[i] = real(v[i]);
return result;
}
} // namespace
ProfileHelper::ProfileHelper(const SliceStack& stack)
: m_stack(stack)
{
}
std::vector<complex_t> ProfileHelper::profile(const std::vector<double>& z_values,
std::string component) const
{
const complex_t top_value =
!m_stack.empty() ? quantity(m_stack.at(0).material(), component) : 0.0;
std::vector<complex_t> result(z_values.size(), top_value);
for (size_t i = 1; i < m_stack.size(); ++i) {
const Slice& slice = m_stack.at(i);
const Slice& sliceAbove = m_stack.at(i - 1);
const complex_t diff =
quantity(slice.material(), component) - quantity(sliceAbove.material(), component);
for (size_t j = 0; j < z_values.size(); ++j) {
const double arg = (slice.hig() - z_values[j]);
const double t = Transition(arg, slice.topRoughness(), slice.topRMS());
result[j] += diff * t;
}
}
return result;
}
// Note: for refractive index materials, the material interpolation actually happens at the level
// of n^2. To first order in delta and beta, this implies the same smooth interpolation of delta
// and beta, as is done here.
std::vector<complex_t> ProfileHelper::calculateSLDProfile(const std::vector<double>& z_values) const
{
return profile(z_values, SLD);
}
std::vector<double>
ProfileHelper::calculateMagnetizationProfile(const std::vector<double>& z_values,
std::string xyz) const
{
if (xyz != X && xyz != Y && xyz != Z)
throw std::runtime_error("Incorrect magnetization component \"" + xyz + "\".\nOnly \"" + X
+ "\", \"" + Y + "\" or \"" + Z + "\" are allowed.");
return ::vec2real(profile(z_values, xyz));
}
std::pair<double, double> ProfileHelper::defaultLimits() const
{
if (m_stack.size() < 2)
return {0.0, 0.0};
double interface_span = m_stack.front().low() - m_stack.back().hig();
double default_margin = interface_span > 0.0 ? interface_span / 20.0 : 10.0;
const double top_rms = m_stack.at(1).topRMS();
const double bottom_rms = m_stack.back().topRMS();
double top_margin = top_rms > 0 ? 5.0 * top_rms : default_margin;
double bottom_margin = bottom_rms > 0 ? 5.0 * bottom_rms : default_margin;
double z_min = m_stack.back().hig() - bottom_margin;
double z_max = m_stack.front().low() + top_margin;
return {z_min, z_max};
}
|
