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
|
// ************************************************************************************************
//
// BornAgain: simulate and fit reflection and scattering
//
//! @file Resample/Flux/MatrixFlux.cpp
//! @brief Implements class MatrixFlux.
//!
//! @homepage http://www.bornagainproject.org
//! @license GNU General Public License v3 or higher (see COPYING)
//! @copyright Forschungszentrum Jülich GmbH 2020
//! @authors Scientific Computing Group at MLZ (see CITATION, AUTHORS)
//
// ************************************************************************************************
#include "Resample/Flux/MatrixFlux.h"
#include "Base/Util/Assert.h"
namespace {
const auto eps = std::numeric_limits<double>::epsilon() * 10.;
complex_t GetImExponential(complex_t exponent)
{
if (exponent.imag() > -std::log(std::numeric_limits<double>::min()))
return 0.0;
return std::exp(I * exponent);
}
} // namespace
MatrixFlux::MatrixFlux(double kz_sign, const Spinor& k_eigen, const R3& b, double magnetic_SLD)
: m_k_eigen(k_eigen)
, m_T(+1, 0, 0, +1)
, m_R(-1, 0, 0, -1)
, m_b(b)
, m_kz_sign(kz_sign)
, m_magnetic_SLD(magnetic_SLD)
{
ASSERT(std::abs(m_b.mag() - 1) < eps || (m_b.mag() < eps && magnetic_SLD < eps));
}
//! See PhysRef, chapter "Polarized", section "Wavenumber operator".
SpinMatrix MatrixFlux::computeKappa() const
{
const complex_t alpha = m_k_eigen.u + m_k_eigen.v;
const complex_t beta = m_k_eigen.u - m_k_eigen.v;
return SpinMatrix(alpha + beta * m_b.z(), beta * (m_b.x() - I * m_b.y()),
beta * (m_b.x() + I * m_b.y()), alpha - beta * m_b.z())
/ 2.;
}
//! See PhysRef, chapter "Polarized", section "Wavenumber operator".
SpinMatrix MatrixFlux::computeInverseKappa() const
{
const complex_t alpha = m_k_eigen.u + m_k_eigen.v;
const complex_t beta = m_k_eigen.u - m_k_eigen.v;
if (std::abs(alpha * alpha - beta * beta) == 0.)
throw std::runtime_error("Inverse wavenumber operator is singular for given B field."
" Add some absorption to the nuclear SLD.");
return 2. / (alpha * alpha - beta * beta)
* SpinMatrix(alpha - beta * m_b.z(), beta * (-m_b.x() + I * m_b.y()),
-beta * (m_b.x() + I * m_b.y()), alpha + beta * m_b.z());
}
//! See PhysRef, chapter "Polarized", section "Eigendecomposition"
SpinMatrix MatrixFlux::eigenToMatrix(const Spinor& diagonal) const
{
const SpinMatrix M(diagonal.u, 0, 0, diagonal.v);
if (std::abs(m_b.mag() - 1.) < eps) {
SpinMatrix Q(1. + m_b.z(), m_b.x() - I * m_b.y(), m_b.x() + I * m_b.y(), -m_b.z() - 1.);
return Q * M * Q.adjoint() / (2. * (1. + m_b.z()));
}
ASSERT(m_b.mag() < eps); // remaining case: no field
return M;
}
SpinMatrix MatrixFlux::computeDeltaMatrix(double thickness) const
{
return eigenToMatrix({GetImExponential(m_kz_sign * thickness * m_k_eigen.u),
GetImExponential(m_kz_sign * thickness * m_k_eigen.v)});
}
SpinMatrix MatrixFlux::T1Matrix() const
{
return eigenToMatrix({1., 0.});
}
SpinMatrix MatrixFlux::T2Matrix() const
{
return eigenToMatrix({0., 1.});
}
Spinor MatrixFlux::T1plus() const
{
return T1Matrix() * m_T.col0();
}
Spinor MatrixFlux::R1plus() const
{
return T1Matrix() * m_R.col0();
}
Spinor MatrixFlux::T2plus() const
{
return T2Matrix() * m_T.col0();
}
Spinor MatrixFlux::R2plus() const
{
return T2Matrix() * m_R.col0();
}
Spinor MatrixFlux::T1min() const
{
return T1Matrix() * m_T.col1();
}
Spinor MatrixFlux::R1min() const
{
return T1Matrix() * m_R.col1();
}
Spinor MatrixFlux::T2min() const
{
return T2Matrix() * m_T.col1();
}
Spinor MatrixFlux::R2min() const
{
return T2Matrix() * m_R.col1();
}
Spinor MatrixFlux::getKz() const
{
return m_kz_sign * m_k_eigen;
}
|
