Skip to content
Snippets Groups Projects
Commit b9bda0b7 authored by Van Herck, Walter's avatar Van Herck, Walter
Browse files

Took care of zero eigenvalues in magnetic matrix calculation

parent 8fec2e4f
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Core/Algorithms/inc/SpecularMagnetic.h
+ 38
− 1
View file @ b9bda0b7
......@@ -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;
}
......
Core/Algorithms/src/SpecularMagnetic.cpp
+ 148
− 89
View file @ b9bda0b7
......@@ -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));
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment