diff --git a/Core/Algorithms/inc/SpecularMagnetic.h b/Core/Algorithms/inc/SpecularMagnetic.h
index 376328000e442075260813d03f1a9c3f77670dda..94b3433d6e55c7745d77d4989dea65fe0e170d89 100644
--- a/Core/Algorithms/inc/SpecularMagnetic.h
+++ b/Core/Algorithms/inc/SpecularMagnetic.h
@@ -33,7 +33,7 @@ public:
//! layer coefficients for matrix formalism
class LayerMatrixCoeff {
public:
- LayerMatrixCoeff() {}
+ LayerMatrixCoeff() : m_kt(0.0) {}
~LayerMatrixCoeff() {}
Eigen::Vector2cd T1plus() const;
Eigen::Vector2cd R1plus() const;
@@ -57,6 +57,7 @@ public:
complex_t m_a; // polarization independent part
complex_t m_b_mag; // magnitude of polarization part
complex_t m_bz; // z-part of polarization scattering
+ double m_kt; // wavevector length times thickness of layer for use when lambda=0
void calculateTRMatrices();
void initializeBottomLayerPhiPsi();
private:
@@ -100,6 +101,9 @@ inline Eigen::Vector2cd SpecularMagnetic::LayerMatrixCoeff::T1plus() const {
Eigen::Vector2cd result;
result(0) = T1m.row(2).dot(phi_psi_plus);
result(1) = T1m.row(3).dot(phi_psi_plus);
+ if (lambda(0)==0.0 && result==Eigen::Vector2cd::Zero()) {
+ result(0) = 0.5;
+ }
return result;
}
@@ -107,6 +111,12 @@ inline Eigen::Vector2cd SpecularMagnetic::LayerMatrixCoeff::R1plus() const {
Eigen::Vector2cd result;
result(0) = R1m.row(2).dot(phi_psi_plus);
result(1) = R1m.row(3).dot(phi_psi_plus);
+ if (lambda(0)==0.0) {
+ if (T1m.row(2).dot(phi_psi_plus)==0.0
+ && T1m.row(3).dot(phi_psi_plus)==0.0) {
+ result(0) = -0.5;
+ }
+ }
return result;
}
@@ -114,6 +124,9 @@ inline Eigen::Vector2cd SpecularMagnetic::LayerMatrixCoeff::T2plus() const {
Eigen::Vector2cd result;
result(0) = T2m.row(2).dot(phi_psi_plus);
result(1) = T2m.row(3).dot(phi_psi_plus);
+ if (lambda(1)==0.0 && result==Eigen::Vector2cd::Zero()) {
+ result(0) = 0.5;
+ }
return result;
}
@@ -121,6 +134,12 @@ inline Eigen::Vector2cd SpecularMagnetic::LayerMatrixCoeff::R2plus() const {
Eigen::Vector2cd result;
result(0) = R2m.row(2).dot(phi_psi_plus);
result(1) = R2m.row(3).dot(phi_psi_plus);
+ if (lambda(1)==0.0) {
+ if (T2m.row(2).dot(phi_psi_plus)==0.0
+ && T2m.row(3).dot(phi_psi_plus)==0.0) {
+ result(0) = -0.5;
+ }
+ }
return result;
}
@@ -128,6 +147,9 @@ inline Eigen::Vector2cd SpecularMagnetic::LayerMatrixCoeff::T1min() const {
Eigen::Vector2cd result;
result(0) = T1m.row(2).dot(phi_psi_min);
result(1) = T1m.row(3).dot(phi_psi_min);
+ if (lambda(0)==0.0 && result==Eigen::Vector2cd::Zero()) {
+ result(1) = 0.5;
+ }
return result;
}
@@ -135,6 +157,12 @@ inline Eigen::Vector2cd SpecularMagnetic::LayerMatrixCoeff::R1min() const {
Eigen::Vector2cd result;
result(0) = R1m.row(2).dot(phi_psi_min);
result(1) = R1m.row(3).dot(phi_psi_min);
+ if (lambda(0)==0.0) {
+ if (T1m.row(2).dot(phi_psi_min)==0.0
+ && T1m.row(3).dot(phi_psi_min)==0.0) {
+ result(1) = -0.5;
+ }
+ }
return result;
}
@@ -142,6 +170,9 @@ inline Eigen::Vector2cd SpecularMagnetic::LayerMatrixCoeff::T2min() const {
Eigen::Vector2cd result;
result(0) = T2m.row(2).dot(phi_psi_min);
result(1) = T2m.row(3).dot(phi_psi_min);
+ if (lambda(1)==0.0 && result==Eigen::Vector2cd::Zero()) {
+ result(1) = 0.5;
+ }
return result;
}
@@ -149,6 +180,12 @@ inline Eigen::Vector2cd SpecularMagnetic::LayerMatrixCoeff::R2min() const {
Eigen::Vector2cd result;
result(0) = R2m.row(2).dot(phi_psi_min);
result(1) = R2m.row(3).dot(phi_psi_min);
+ if (lambda(1)==0.0) {
+ if (T2m.row(2).dot(phi_psi_min)==0.0
+ && T2m.row(3).dot(phi_psi_min)==0.0) {
+ result(1) = -0.5;
+ }
+ }
return result;
}
diff --git a/Core/Algorithms/src/SpecularMagnetic.cpp b/Core/Algorithms/src/SpecularMagnetic.cpp
index 0b83c6170b3d823d287c4ce228ce9df6cbe2f692..14f91e782b3c5b3eb8500bdbd07eaa52eb303d9d 100644
--- a/Core/Algorithms/src/SpecularMagnetic.cpp
+++ b/Core/Algorithms/src/SpecularMagnetic.cpp
@@ -35,6 +35,7 @@ void SpecularMagnetic::calculateEigenvalues(const MultiLayer& sample,
for(size_t i=0; i<coeff.size(); ++i) {
coeff[i].m_scatt_matrix = sample.getLayer(i)->getMaterial()->
getScatteringMatrix(k);
+ coeff[i].m_kt = mag_k*sample.getLayer(i)->getThickness();
coeff[i].m_a = coeff[i].m_scatt_matrix.trace()/2.0;
coeff[i].m_b_mag = std::sqrt(coeff[i].m_a*coeff[i].m_a -
coeff[i].m_scatt_matrix.determinant());
@@ -62,8 +63,8 @@ void SpecularMagnetic::calculateTransferAndBoundary(const MultiLayer& sample,
coeff[0].calculateTRMatrices();
coeff[0].l.setIdentity();
for (int i=(int)N-2; i>0; --i) {
- coeff[i].calculateTRMatrices();
double t = sample.getLayer(i)->getThickness();
+ coeff[i].calculateTRMatrices();
coeff[i].l =
coeff[i].R1m * getImExponential((complex_t)(coeff[i].kz(0)*t)) +
coeff[i].T1m * getImExponential((complex_t)(-coeff[i].kz(0)*t)) +
@@ -135,93 +136,137 @@ void SpecularMagnetic::LayerMatrixCoeff::calculateTRMatrices()
return;
}
- // T1m:
- // row 0:
- T1m(0,0) = (1.0 - m_bz/m_b_mag)/4.0;
- T1m(0,1) = - m_scatt_matrix(0,1)/4.0/m_b_mag;
- T1m(0,2) = lambda(0)*(m_bz/m_b_mag - 1.0)/4.0;
- T1m(0,3) = m_scatt_matrix(0,1)*lambda(0)/4.0/m_b_mag;
- // row 1:
- T1m(1,0) = - m_scatt_matrix(1,0)/4.0/m_b_mag;
- T1m(1,1) = (1.0 + m_bz/m_b_mag)/4.0;
- T1m(1,2) = m_scatt_matrix(1,0)*lambda(0)/4.0/m_b_mag;
- T1m(1,3) = - lambda(0)*(m_bz/m_b_mag + 1.0)/4.0;
- // row 2:
- T1m(2,0) = -(1.0 - m_bz/m_b_mag)/4.0/lambda(0);
- T1m(2,1) = m_scatt_matrix(0,1)/4.0/m_b_mag/lambda(0);
- T1m(2,2) = -(m_bz/m_b_mag - 1.0)/4.0;
- T1m(2,3) = - m_scatt_matrix(0,1)/4.0/m_b_mag;
- // row 3:
- T1m(3,0) = m_scatt_matrix(1,0)/4.0/m_b_mag/lambda(0);
- T1m(3,1) = -(1.0 + m_bz/m_b_mag)/4.0/lambda(0);
- T1m(3,2) = - m_scatt_matrix(1,0)/4.0/m_b_mag;
- T1m(3,3) = (m_bz/m_b_mag + 1.0)/4.0;
+ if (lambda(0)==0.0) {
+ complex_t ikt = complex_t(0.0, 1.0) * m_kt;
+ // Lambda1 component contained only in T1m (R1m=0)
+ // row 0:
+ T1m(0,0) = (1.0 - m_bz/m_b_mag)/2.0;
+ T1m(0,1) = - m_scatt_matrix(0,1)/2.0/m_b_mag;
+ // row 1:
+ T1m(1,0) = - m_scatt_matrix(1,0)/2.0/m_b_mag;
+ T1m(1,1) = (1.0 + m_bz/m_b_mag)/2.0;
+ // row 2:
+ T1m(2,0) = ikt*(1.0 - m_bz/m_b_mag)/2.0;
+ T1m(2,1) = -ikt*m_scatt_matrix(0,1)/2.0/m_b_mag;
+ T1m(2,2) = T1m(0,0);
+ T1m(2,3) = T1m(0,1);
+ // row 3:
+ T1m(3,0) = -ikt*m_scatt_matrix(1,0)/2.0/m_b_mag;
+ T1m(3,1) = ikt*(1.0 + m_bz/m_b_mag)/2.0;
+ T1m(3,2) = T1m(1,0);
+ T1m(3,3) = T1m(1,1);
+ }
+ else {
+ // T1m:
+ // row 0:
+ T1m(0,0) = (1.0 - m_bz/m_b_mag)/4.0;
+ T1m(0,1) = - m_scatt_matrix(0,1)/4.0/m_b_mag;
+ T1m(0,2) = lambda(0)*(m_bz/m_b_mag - 1.0)/4.0;
+ T1m(0,3) = m_scatt_matrix(0,1)*lambda(0)/4.0/m_b_mag;
+ // row 1:
+ T1m(1,0) = - m_scatt_matrix(1,0)/4.0/m_b_mag;
+ T1m(1,1) = (1.0 + m_bz/m_b_mag)/4.0;
+ T1m(1,2) = m_scatt_matrix(1,0)*lambda(0)/4.0/m_b_mag;
+ T1m(1,3) = - lambda(0)*(m_bz/m_b_mag + 1.0)/4.0;
+ // row 2:
+ T1m(2,0) = -(1.0 - m_bz/m_b_mag)/4.0/lambda(0);
+ T1m(2,1) = m_scatt_matrix(0,1)/4.0/m_b_mag/lambda(0);
+ T1m(2,2) = -(m_bz/m_b_mag - 1.0)/4.0;
+ T1m(2,3) = - m_scatt_matrix(0,1)/4.0/m_b_mag;
+ // row 3:
+ T1m(3,0) = m_scatt_matrix(1,0)/4.0/m_b_mag/lambda(0);
+ T1m(3,1) = -(1.0 + m_bz/m_b_mag)/4.0/lambda(0);
+ T1m(3,2) = - m_scatt_matrix(1,0)/4.0/m_b_mag;
+ T1m(3,3) = (m_bz/m_b_mag + 1.0)/4.0;
- // R1m:
- // row 0:
- R1m(0,0) = T1m(0,0);
- R1m(0,1) = T1m(0,1);
- R1m(0,2) = -T1m(0,2);
- R1m(0,3) = -T1m(0,3);
- // row 1:
- R1m(1,0) = T1m(1,0);
- R1m(1,1) = T1m(1,1);
- R1m(1,2) = -T1m(1,2);
- R1m(1,3) = -T1m(1,3);
- // row 2:
- R1m(2,0) = -T1m(2,0);
- R1m(2,1) = -T1m(2,1);
- R1m(2,2) = T1m(2,2);
- R1m(2,3) = T1m(2,3);
- // row 3:
- R1m(3,0) = -T1m(3,0);
- R1m(3,1) = -T1m(3,1);
- R1m(3,2) = T1m(3,2);
- R1m(3,3) = T1m(3,3);
+ // R1m:
+ // row 0:
+ R1m(0,0) = T1m(0,0);
+ R1m(0,1) = T1m(0,1);
+ R1m(0,2) = -T1m(0,2);
+ R1m(0,3) = -T1m(0,3);
+ // row 1:
+ R1m(1,0) = T1m(1,0);
+ R1m(1,1) = T1m(1,1);
+ R1m(1,2) = -T1m(1,2);
+ R1m(1,3) = -T1m(1,3);
+ // row 2:
+ R1m(2,0) = -T1m(2,0);
+ R1m(2,1) = -T1m(2,1);
+ R1m(2,2) = T1m(2,2);
+ R1m(2,3) = T1m(2,3);
+ // row 3:
+ R1m(3,0) = -T1m(3,0);
+ R1m(3,1) = -T1m(3,1);
+ R1m(3,2) = T1m(3,2);
+ R1m(3,3) = T1m(3,3);
+ }
- // T2m:
- // row 0:
- T2m(0,0) = (1.0 + m_bz/m_b_mag)/4.0;
- T2m(0,1) = m_scatt_matrix(0,1)/4.0/m_b_mag;
- T2m(0,2) = - lambda(1)*(m_bz/m_b_mag + 1.0)/4.0;
- T2m(0,3) = - m_scatt_matrix(0,1)*lambda(1)/4.0/m_b_mag;
- // row 1:
- T2m(1,0) = m_scatt_matrix(1,0)/4.0/m_b_mag;
- T2m(1,1) = (1.0 - m_bz/m_b_mag)/4.0;
- T2m(1,2) = - m_scatt_matrix(1,0)*lambda(1)/4.0/m_b_mag;
- T2m(1,3) = lambda(1)*(m_bz/m_b_mag - 1.0)/4.0;
- // row 2:
- T2m(2,0) = -(1.0 + m_bz/m_b_mag)/4.0/lambda(1);
- T2m(2,1) = - m_scatt_matrix(0,1)/4.0/m_b_mag/lambda(1);
- T2m(2,2) = (m_bz/m_b_mag + 1.0)/4.0;
- T2m(2,3) = m_scatt_matrix(0,1)/4.0/m_b_mag;
- // row 3:
- T2m(3,0) = - m_scatt_matrix(1,0)/4.0/m_b_mag/lambda(1);
- T2m(3,1) = -(1.0 - m_bz/m_b_mag)/4.0/lambda(1);
- T2m(3,2) = m_scatt_matrix(1,0)/4.0/m_b_mag;
- T2m(3,3) = (1.0 - m_bz/m_b_mag)/4.0;
+ if (lambda(1)==0.0) {
+ complex_t ikt = complex_t(0.0, 1.0) * m_kt;
+ // Lambda2 component contained only in T2m (R2m=0)
+ // row 0:
+ T2m(0,0) = (1.0 + m_bz/m_b_mag)/2.0;
+ T2m(0,1) = m_scatt_matrix(0,1)/2.0/m_b_mag;
+ // row 1:
+ T2m(1,0) = m_scatt_matrix(1,0)/2.0/m_b_mag;
+ T2m(1,1) = (1.0 - m_bz/m_b_mag)/2.0;
+ // row 2:
+ T2m(2,0) = ikt*(1.0 + m_bz/m_b_mag)/2.0;
+ T2m(2,1) = ikt*m_scatt_matrix(0,1)/2.0/m_b_mag;
+ T2m(2,2) = T2m(0,0);
+ T2m(2,3) = T2m(0,1);
+ // row 3:
+ T2m(3,0) = ikt*m_scatt_matrix(1,0)/2.0/m_b_mag;
+ T2m(3,1) = ikt*(1.0 - m_bz/m_b_mag)/2.0;
+ T2m(3,2) = T2m(1,0);
+ T2m(3,3) = T2m(1,1);
+ }
+ else {
+ // T2m:
+ // row 0:
+ T2m(0,0) = (1.0 + m_bz/m_b_mag)/4.0;
+ T2m(0,1) = m_scatt_matrix(0,1)/4.0/m_b_mag;
+ T2m(0,2) = - lambda(1)*(m_bz/m_b_mag + 1.0)/4.0;
+ T2m(0,3) = - m_scatt_matrix(0,1)*lambda(1)/4.0/m_b_mag;
+ // row 1:
+ T2m(1,0) = m_scatt_matrix(1,0)/4.0/m_b_mag;
+ T2m(1,1) = (1.0 - m_bz/m_b_mag)/4.0;
+ T2m(1,2) = - m_scatt_matrix(1,0)*lambda(1)/4.0/m_b_mag;
+ T2m(1,3) = lambda(1)*(m_bz/m_b_mag - 1.0)/4.0;
+ // row 2:
+ T2m(2,0) = -(1.0 + m_bz/m_b_mag)/4.0/lambda(1);
+ T2m(2,1) = - m_scatt_matrix(0,1)/4.0/m_b_mag/lambda(1);
+ T2m(2,2) = (m_bz/m_b_mag + 1.0)/4.0;
+ T2m(2,3) = m_scatt_matrix(0,1)/4.0/m_b_mag;
+ // row 3:
+ T2m(3,0) = - m_scatt_matrix(1,0)/4.0/m_b_mag/lambda(1);
+ T2m(3,1) = -(1.0 - m_bz/m_b_mag)/4.0/lambda(1);
+ T2m(3,2) = m_scatt_matrix(1,0)/4.0/m_b_mag;
+ T2m(3,3) = (1.0 - m_bz/m_b_mag)/4.0;
- // R2m:
- // row 0:
- R2m(0,0) = T2m(0,0);
- R2m(0,1) = T2m(0,1);
- R2m(0,2) = -T2m(0,2);
- R2m(0,3) = -T2m(0,3);
- // row 1:
- R2m(1,0) = T2m(1,0);
- R2m(1,1) = T2m(1,1);
- R2m(1,2) = -T2m(1,2);
- R2m(1,3) = -T2m(1,3);
- // row 2:
- R2m(2,0) = -T2m(2,0);
- R2m(2,1) = -T2m(2,1);
- R2m(2,2) = T2m(2,2);
- R2m(2,3) = T2m(2,3);
- // row 3:
- R2m(3,0) = -T2m(3,0);
- R2m(3,1) = -T2m(3,1);
- R2m(3,2) = T2m(3,2);
- R2m(3,3) = T2m(3,3);
+ // R2m:
+ // row 0:
+ R2m(0,0) = T2m(0,0);
+ R2m(0,1) = T2m(0,1);
+ R2m(0,2) = -T2m(0,2);
+ R2m(0,3) = -T2m(0,3);
+ // row 1:
+ R2m(1,0) = T2m(1,0);
+ R2m(1,1) = T2m(1,1);
+ R2m(1,2) = -T2m(1,2);
+ R2m(1,3) = -T2m(1,3);
+ // row 2:
+ R2m(2,0) = -T2m(2,0);
+ R2m(2,1) = -T2m(2,1);
+ R2m(2,2) = T2m(2,2);
+ R2m(2,3) = T2m(2,3);
+ // row 3:
+ R2m(3,0) = -T2m(3,0);
+ R2m(3,1) = -T2m(3,1);
+ R2m(3,2) = T2m(3,2);
+ R2m(3,3) = T2m(3,3);
+ }
}
void SpecularMagnetic::LayerMatrixCoeff::initializeBottomLayerPhiPsi()
@@ -250,29 +295,43 @@ void SpecularMagnetic::LayerMatrixCoeff::initializeBottomLayerPhiPsi()
void SpecularMagnetic::LayerMatrixCoeff::calculateTRWithoutMagnetization()
{
- // T1m:
T1m.setZero();
+ R1m.setZero();
+ T2m.setZero();
+ R2m.setZero();
+
+ if (m_a==0.0) {
+ // Spin down component contained only in T1 (R1=0)
+ T1m(1,1) = 1.0;
+ T1m(3,1) = complex_t(0.0, 1.0)*m_kt;
+ T1m(3,3) = 1.0;
+
+ // Spin up component contained only in T2 (R2=0)
+ T2m(0,0) = 1.0;
+ T2m(2,0) = complex_t(0.0, 1.0)*m_kt;
+ T2m(2,2) = 1.0;
+ return;
+ }
+
+ // T1m:
T1m(1,1) = 0.5;
T1m(1,3) = -std::sqrt(m_a)/2.0;
T1m(3,1) = -1.0/(2.0*std::sqrt(m_a));
T1m(3,3) = 0.5;
// R1m:
- R1m.setZero();
R1m(1,1) = 0.5;
R1m(1,3) = std::sqrt(m_a)/2.0;
R1m(3,1) = 1.0/(2.0*std::sqrt(m_a));
R1m(3,3) = 0.5;
// T2m:
- T2m.setZero();
T2m(0,0) = 0.5;
T2m(0,2) = -std::sqrt(m_a)/2.0;
T2m(2,0) = -1.0/(2.0*std::sqrt(m_a));
T2m(2,2) = 0.5;
// R2m:
- R2m.setZero();
R2m(0,0) = 0.5;
R2m(0,2) = std::sqrt(m_a)/2.0;
R2m(2,0) = 1.0/(2.0*std::sqrt(m_a));