BornAgain  1.18.0
Simulate and fit neutron and x-ray scattering at grazing incidence
MatrixRTCoefficients.cpp
Go to the documentation of this file.
1 // ************************************************************************** //
2 //
3 // BornAgain: simulate and fit scattering at grazing incidence
4 //
5 //! @file Sample/RT/MatrixRTCoefficients.cpp
6 //! @brief Implements class MatrixRTCoefficients.
7 //!
8 //! @homepage http://www.bornagainproject.org
9 //! @license GNU General Public License v3 or higher (see COPYING)
10 //! @copyright Forschungszentrum Jülich GmbH 2018
11 //! @authors Scientific Computing Group at MLZ (see CITATION, AUTHORS)
12 //
13 // ************************************************************************** //
14 
15 #include "Sample/RT/MatrixRTCoefficients.h"
16 
17 MatrixRTCoefficients* MatrixRTCoefficients::clone() const
18 {
19  return new MatrixRTCoefficients(*this);
20 }
21 
22 Eigen::Vector2cd MatrixRTCoefficients::T1plus() const
23 {
24  Eigen::Vector2cd result;
25  Eigen::Matrix<complex_t, 4, 1> m = T1m * phi_psi_plus;
26  result(0) = m(2);
27  result(1) = m(3);
28  if (lambda(0) == 0.0 && result == Eigen::Vector2cd::Zero())
29  result(0) = 0.5;
30  return result;
31 }
32 
33 Eigen::Vector2cd MatrixRTCoefficients::R1plus() const
34 {
35  Eigen::Vector2cd result;
36  Eigen::Matrix<complex_t, 4, 1> m = R1m * phi_psi_plus;
37  result(0) = m(2);
38  result(1) = m(3);
39  Eigen::Matrix<complex_t, 4, 1> mT = T1m * phi_psi_plus;
40  if (lambda(0) == 0.0 && mT(2) == 0.0 && mT(3) == 0.0)
41  result(0) = -0.5;
42  return result;
43 }
44 
45 Eigen::Vector2cd MatrixRTCoefficients::T2plus() const
46 {
47  Eigen::Vector2cd result;
48  Eigen::Matrix<complex_t, 4, 1> m = T2m * phi_psi_plus;
49  result(0) = m(2);
50  result(1) = m(3);
51  if (lambda(1) == 0.0 && result == Eigen::Vector2cd::Zero())
52  result(0) = 0.5;
53  return result;
54 }
55 
56 Eigen::Vector2cd MatrixRTCoefficients::R2plus() const
57 {
58  Eigen::Vector2cd result;
59  Eigen::Matrix<complex_t, 4, 1> m = R2m * phi_psi_plus;
60  result(0) = m(2);
61  result(1) = m(3);
62  Eigen::Matrix<complex_t, 4, 1> mT = T2m * phi_psi_plus;
63  if (lambda(1) == 0.0 && mT(2) == 0.0 && mT(3) == 0.0)
64  result(0) = -0.5;
65  return result;
66 }
67 
68 Eigen::Vector2cd MatrixRTCoefficients::T1min() const
69 {
70  Eigen::Vector2cd result;
71  Eigen::Matrix<complex_t, 4, 1> m = T1m * phi_psi_min;
72  result(0) = m(2);
73  result(1) = m(3);
74  if (lambda(0) == 0.0 && result == Eigen::Vector2cd::Zero())
75  result(1) = 0.5;
76  return result;
77 }
78 
79 Eigen::Vector2cd MatrixRTCoefficients::R1min() const
80 {
81  Eigen::Vector2cd result;
82  Eigen::Matrix<complex_t, 4, 1> m = R1m * phi_psi_min;
83  result(0) = m(2);
84  result(1) = m(3);
85  Eigen::Matrix<complex_t, 4, 1> mT = T1m * phi_psi_min;
86  if (lambda(0) == 0.0 && mT(2) == 0.0 && mT(3) == 0.0)
87  result(1) = -0.5;
88  return result;
89 }
90 
91 Eigen::Vector2cd MatrixRTCoefficients::T2min() const
92 {
93  Eigen::Vector2cd result;
94  Eigen::Matrix<complex_t, 4, 1> m = T2m * phi_psi_min;
95  result(0) = m(2);
96  result(1) = m(3);
97  if (lambda(1) == 0.0 && result == Eigen::Vector2cd::Zero())
98  result(1) = 0.5;
99  return result;
100 }
101 
102 Eigen::Vector2cd MatrixRTCoefficients::R2min() const
103 {
104  Eigen::Vector2cd result;
105  Eigen::Matrix<complex_t, 4, 1> m = R2m * phi_psi_min;
106  result(0) = m(2);
107  result(1) = m(3);
108  Eigen::Matrix<complex_t, 4, 1> mT = T2m * phi_psi_min;
109  if (lambda(1) == 0.0 && mT(2) == 0.0 && mT(3) == 0.0)
110  result(1) = -0.5;
111  return result;
112 }
113 
114 void MatrixRTCoefficients::calculateTRMatrices()
115 {
116  if (m_b_mag == 0.0) {
117  calculateTRWithoutMagnetization();
118  return;
119  }
120 
121  if (lambda(0) == 0.0) {
122  complex_t ikt = mul_I(m_kt);
123  // Lambda1 component contained only in T1m (R1m=0)
124  // row 0:
125  T1m(0, 0) = (1.0 - m_bz / m_b_mag) / 2.0;
126  T1m(0, 1) = -m_scatt_matrix(0, 1) / 2.0 / m_b_mag;
127  // row 1:
128  T1m(1, 0) = -m_scatt_matrix(1, 0) / 2.0 / m_b_mag;
129  T1m(1, 1) = (1.0 + m_bz / m_b_mag) / 2.0;
130  // row 2:
131  T1m(2, 0) = ikt * (1.0 - m_bz / m_b_mag) / 2.0;
132  T1m(2, 1) = -ikt * m_scatt_matrix(0, 1) / 2.0 / m_b_mag;
133  T1m(2, 2) = T1m(0, 0);
134  T1m(2, 3) = T1m(0, 1);
135  // row 3:
136  T1m(3, 0) = -ikt * m_scatt_matrix(1, 0) / 2.0 / m_b_mag;
137  T1m(3, 1) = ikt * (1.0 + m_bz / m_b_mag) / 2.0;
138  T1m(3, 2) = T1m(1, 0);
139  T1m(3, 3) = T1m(1, 1);
140  } else {
141  // T1m:
142  // row 0:
143  T1m(0, 0) = (1.0 - m_bz / m_b_mag) / 4.0;
144  T1m(0, 1) = -m_scatt_matrix(0, 1) / 4.0 / m_b_mag;
145  T1m(0, 2) = lambda(0) * (m_bz / m_b_mag - 1.0) / 4.0;
146  T1m(0, 3) = m_scatt_matrix(0, 1) * lambda(0) / 4.0 / m_b_mag;
147  // row 1:
148  T1m(1, 0) = -m_scatt_matrix(1, 0) / 4.0 / m_b_mag;
149  T1m(1, 1) = (1.0 + m_bz / m_b_mag) / 4.0;
150  T1m(1, 2) = m_scatt_matrix(1, 0) * lambda(0) / 4.0 / m_b_mag;
151  T1m(1, 3) = -lambda(0) * (m_bz / m_b_mag + 1.0) / 4.0;
152  // row 2:
153  T1m(2, 0) = -(1.0 - m_bz / m_b_mag) / 4.0 / lambda(0);
154  T1m(2, 1) = m_scatt_matrix(0, 1) / 4.0 / m_b_mag / lambda(0);
155  T1m(2, 2) = -(m_bz / m_b_mag - 1.0) / 4.0;
156  T1m(2, 3) = -m_scatt_matrix(0, 1) / 4.0 / m_b_mag;
157  // row 3:
158  T1m(3, 0) = m_scatt_matrix(1, 0) / 4.0 / m_b_mag / lambda(0);
159  T1m(3, 1) = -(1.0 + m_bz / m_b_mag) / 4.0 / lambda(0);
160  T1m(3, 2) = -m_scatt_matrix(1, 0) / 4.0 / m_b_mag;
161  T1m(3, 3) = (m_bz / m_b_mag + 1.0) / 4.0;
162 
163  // R1m:
164  // row 0:
165  R1m(0, 0) = T1m(0, 0);
166  R1m(0, 1) = T1m(0, 1);
167  R1m(0, 2) = -T1m(0, 2);
168  R1m(0, 3) = -T1m(0, 3);
169  // row 1:
170  R1m(1, 0) = T1m(1, 0);
171  R1m(1, 1) = T1m(1, 1);
172  R1m(1, 2) = -T1m(1, 2);
173  R1m(1, 3) = -T1m(1, 3);
174  // row 2:
175  R1m(2, 0) = -T1m(2, 0);
176  R1m(2, 1) = -T1m(2, 1);
177  R1m(2, 2) = T1m(2, 2);
178  R1m(2, 3) = T1m(2, 3);
179  // row 3:
180  R1m(3, 0) = -T1m(3, 0);
181  R1m(3, 1) = -T1m(3, 1);
182  R1m(3, 2) = T1m(3, 2);
183  R1m(3, 3) = T1m(3, 3);
184  }
185 
186  if (lambda(1) == 0.0) {
187  complex_t ikt = mul_I(m_kt);
188  // Lambda2 component contained only in T2m (R2m=0)
189  // row 0:
190  T2m(0, 0) = (1.0 + m_bz / m_b_mag) / 2.0;
191  T2m(0, 1) = m_scatt_matrix(0, 1) / 2.0 / m_b_mag;
192  // row 1:
193  T2m(1, 0) = m_scatt_matrix(1, 0) / 2.0 / m_b_mag;
194  T2m(1, 1) = (1.0 - m_bz / m_b_mag) / 2.0;
195  // row 2:
196  T2m(2, 0) = ikt * (1.0 + m_bz / m_b_mag) / 2.0;
197  T2m(2, 1) = ikt * m_scatt_matrix(0, 1) / 2.0 / m_b_mag;
198  T2m(2, 2) = T2m(0, 0);
199  T2m(2, 3) = T2m(0, 1);
200  // row 3:
201  T2m(3, 0) = ikt * m_scatt_matrix(1, 0) / 2.0 / m_b_mag;
202  T2m(3, 1) = ikt * (1.0 - m_bz / m_b_mag) / 2.0;
203  T2m(3, 2) = T2m(1, 0);
204  T2m(3, 3) = T2m(1, 1);
205  } else {
206  // T2m:
207  // row 0:
208  T2m(0, 0) = (1.0 + m_bz / m_b_mag) / 4.0;
209  T2m(0, 1) = m_scatt_matrix(0, 1) / 4.0 / m_b_mag;
210  T2m(0, 2) = -lambda(1) * (m_bz / m_b_mag + 1.0) / 4.0;
211  T2m(0, 3) = -m_scatt_matrix(0, 1) * lambda(1) / 4.0 / m_b_mag;
212  // row 1:
213  T2m(1, 0) = m_scatt_matrix(1, 0) / 4.0 / m_b_mag;
214  T2m(1, 1) = (1.0 - m_bz / m_b_mag) / 4.0;
215  T2m(1, 2) = -m_scatt_matrix(1, 0) * lambda(1) / 4.0 / m_b_mag;
216  T2m(1, 3) = lambda(1) * (m_bz / m_b_mag - 1.0) / 4.0;
217  // row 2:
218  T2m(2, 0) = -(1.0 + m_bz / m_b_mag) / 4.0 / lambda(1);
219  T2m(2, 1) = -m_scatt_matrix(0, 1) / 4.0 / m_b_mag / lambda(1);
220  T2m(2, 2) = (m_bz / m_b_mag + 1.0) / 4.0;
221  T2m(2, 3) = m_scatt_matrix(0, 1) / 4.0 / m_b_mag;
222  // row 3:
223  T2m(3, 0) = -m_scatt_matrix(1, 0) / 4.0 / m_b_mag / lambda(1);
224  T2m(3, 1) = -(1.0 - m_bz / m_b_mag) / 4.0 / lambda(1);
225  T2m(3, 2) = m_scatt_matrix(1, 0) / 4.0 / m_b_mag;
226  T2m(3, 3) = (1.0 - m_bz / m_b_mag) / 4.0;
227 
228  // R2m:
229  // row 0:
230  R2m(0, 0) = T2m(0, 0);
231  R2m(0, 1) = T2m(0, 1);
232  R2m(0, 2) = -T2m(0, 2);
233  R2m(0, 3) = -T2m(0, 3);
234  // row 1:
235  R2m(1, 0) = T2m(1, 0);
236  R2m(1, 1) = T2m(1, 1);
237  R2m(1, 2) = -T2m(1, 2);
238  R2m(1, 3) = -T2m(1, 3);
239  // row 2:
240  R2m(2, 0) = -T2m(2, 0);
241  R2m(2, 1) = -T2m(2, 1);
242  R2m(2, 2) = T2m(2, 2);
243  R2m(2, 3) = T2m(2, 3);
244  // row 3:
245  R2m(3, 0) = -T2m(3, 0);
246  R2m(3, 1) = -T2m(3, 1);
247  R2m(3, 2) = T2m(3, 2);
248  R2m(3, 3) = T2m(3, 3);
249  }
250 }
251 
252 void MatrixRTCoefficients::calculateTRWithoutMagnetization()
253 {
254  T1m.setZero();
255  R1m.setZero();
256  T2m.setZero();
257  R2m.setZero();
258 
259  if (m_a == 0.0) {
260  // Spin down component contained only in T1 (R1=0)
261  T1m(1, 1) = 1.0;
262  T1m(3, 1) = mul_I(m_kt);
263  T1m(3, 3) = 1.0;
264 
265  // Spin up component contained only in T2 (R2=0)
266  T2m(0, 0) = 1.0;
267  T2m(2, 0) = mul_I(m_kt);
268  T2m(2, 2) = 1.0;
269  return;
270  }
271 
272  // T1m:
273  T1m(1, 1) = 0.5;
274  T1m(1, 3) = -std::sqrt(m_a) / 2.0;
275  T1m(3, 1) = -1.0 / (2.0 * std::sqrt(m_a));
276  T1m(3, 3) = 0.5;
277 
278  // R1m:
279  R1m(1, 1) = 0.5;
280  R1m(1, 3) = std::sqrt(m_a) / 2.0;
281  R1m(3, 1) = 1.0 / (2.0 * std::sqrt(m_a));
282  R1m(3, 3) = 0.5;
283 
284  // T2m:
285  T2m(0, 0) = 0.5;
286  T2m(0, 2) = -std::sqrt(m_a) / 2.0;
287  T2m(2, 0) = -1.0 / (2.0 * std::sqrt(m_a));
288  T2m(2, 2) = 0.5;
289 
290  // R2m:
291  R2m(0, 0) = 0.5;
292  R2m(0, 2) = std::sqrt(m_a) / 2.0;
293  R2m(2, 0) = 1.0 / (2.0 * std::sqrt(m_a));
294  R2m(2, 2) = 0.5;
295 }
296 
297 void MatrixRTCoefficients::initializeBottomLayerPhiPsi()
298 {
299  if (m_b_mag == 0.0) {
300  phi_psi_min << 0.0, -std::sqrt(m_a), 0.0, 1.0;
301  phi_psi_plus << -std::sqrt(m_a), 0.0, 1.0, 0.0;
302  return;
303  }
304  // First basis vector that has no upward going wave amplitude
305  phi_psi_min(0) = m_scatt_matrix(0, 1) * (lambda(0) - lambda(1)) / 2.0 / m_b_mag;
306  phi_psi_min(1) = (m_bz * (lambda(1) - lambda(0)) / m_b_mag - lambda(1) - lambda(0)) / 2.0;
307  phi_psi_min(2) = 0.0;
308  phi_psi_min(3) = 1.0;
309 
310  // Second basis vector that has no upward going wave amplitude
311  phi_psi_plus(0) = -(m_scatt_matrix(0, 0) + lambda(0) * lambda(1)) / (lambda(0) + lambda(1));
312  phi_psi_plus(1) = m_scatt_matrix(1, 0) * (lambda(0) - lambda(1)) / 2.0 / m_b_mag;
313  phi_psi_plus(2) = 1.0;
314  phi_psi_plus(3) = 0.0;
315 }
mul_I
complex_t mul_I(complex_t z)
Returns product I*z, where I is the imaginary unit.
Definition: Complex.h:24
complex_t
std::complex< double > complex_t
Definition: Complex.h:20
MatrixRTCoefficients.h
Defines class MatrixRTCoefficients.
MatrixRTCoefficients
Specular reflection and transmission coefficients in a layer in case of 2x2 matrix interactions betwe...
Definition: MatrixRTCoefficients.h:25
MatrixRTCoefficients::m_bz
complex_t m_bz
z-part of magnetic interaction term
Definition: MatrixRTCoefficients.h:65
MatrixRTCoefficients::T2m
Eigen::Matrix4cd T2m
matrix selecting the transmitted part of the second eigenmode
Definition: MatrixRTCoefficients.h:58
MatrixRTCoefficients::m_scatt_matrix
Eigen::Matrix2cd m_scatt_matrix
scattering matrix
Definition: MatrixRTCoefficients.h:62
MatrixRTCoefficients::lambda
Eigen::Vector2cd lambda
positive eigenvalues of transfer matrix
Definition: MatrixRTCoefficients.h:51
MatrixRTCoefficients::m_a
complex_t m_a
polarization independent part of scattering matrix
Definition: MatrixRTCoefficients.h:63
MatrixRTCoefficients::R2m
Eigen::Matrix4cd R2m
matrix selecting the reflected part of the second eigenmode
Definition: MatrixRTCoefficients.h:60
MatrixRTCoefficients::phi_psi_min
Eigen::Vector4cd phi_psi_min
boundary values for down-polarization
Definition: MatrixRTCoefficients.h:53
MatrixRTCoefficients::T1plus
virtual Eigen::Vector2cd T1plus() const
The following functions return the transmitted and reflected amplitudes for different incoming beam p...
Definition: MatrixRTCoefficients.cpp:22
MatrixRTCoefficients::phi_psi_plus
Eigen::Vector4cd phi_psi_plus
boundary values for up-polarization
Definition: MatrixRTCoefficients.h:52
MatrixRTCoefficients::m_kt
double m_kt
wavevector length times thickness of layer for use when lambda=0
Definition: MatrixRTCoefficients.h:66
MatrixRTCoefficients::m_b_mag
complex_t m_b_mag
magnitude of magnetic interaction term
Definition: MatrixRTCoefficients.h:64
MatrixRTCoefficients::T1min
virtual Eigen::Vector2cd T1min() const
Definition: MatrixRTCoefficients.cpp:68
MatrixRTCoefficients::R2plus
virtual Eigen::Vector2cd R2plus() const
Definition: MatrixRTCoefficients.cpp:56
MatrixRTCoefficients::R1min
virtual Eigen::Vector2cd R1min() const
Definition: MatrixRTCoefficients.cpp:79
MatrixRTCoefficients::T2min
virtual Eigen::Vector2cd T2min() const
Definition: MatrixRTCoefficients.cpp:91
MatrixRTCoefficients::T1m
Eigen::Matrix4cd T1m
matrix selecting the transmitted part of the first eigenmode
Definition: MatrixRTCoefficients.h:54
MatrixRTCoefficients::R1m
Eigen::Matrix4cd R1m
matrix selecting the reflected part of the first eigenmode
Definition: MatrixRTCoefficients.h:56
MatrixRTCoefficients::R2min
virtual Eigen::Vector2cd R2min() const
Definition: MatrixRTCoefficients.cpp:102
MatrixRTCoefficients::clone
virtual MatrixRTCoefficients * clone() const
Definition: MatrixRTCoefficients.cpp:17
MatrixRTCoefficients::T2plus
virtual Eigen::Vector2cd T2plus() const
Definition: MatrixRTCoefficients.cpp:45
MatrixRTCoefficients::R1plus
virtual Eigen::Vector2cd R1plus() const
Definition: MatrixRTCoefficients.cpp:33