Skip to content
Snippets Groups Projects
Commit a2a0f17f authored by Wuttke, Joachim's avatar Wuttke, Joachim
Browse files

replace reference data MagneticRotationReflectivityMM_NevotCroce etc...

replace reference data MagneticRotationReflectivityMM_NevotCroce etc (deviation up to 1e-4 in few points)
parent ddd35ea0
No related branches found
No related tags found
1 merge request!1687DevRef polarized ctd; described NC; simplified code
Pipeline #102447 passed
...@@ -32,30 +32,23 @@ std::pair<SpinMatrix, SpinMatrix> Compute::refractionMatrixBlocksNevot(const Mat ...@@ -32,30 +32,23 @@ std::pair<SpinMatrix, SpinMatrix> Compute::refractionMatrixBlocksNevot(const Mat
complex_t beta_b = coeff_b.k_eigen_up() - coeff_b.k_eigen_dn(); complex_t beta_b = coeff_b.k_eigen_up() - coeff_b.k_eigen_dn();
const complex_t alpha = alpha_b + sign * alpha_a; const complex_t alpha = alpha_b + sign * alpha_a;
C3 b_p_vec = beta_b * coeff_b.field() + sign * beta_a * coeff_a.field(); C3 b = beta_b * coeff_b.field() + sign * beta_a * coeff_a.field();
auto square = [](auto& v) { return v.x() * v.x() + v.y() * v.y() + v.z() * v.z(); }; auto square = [](auto& v) { return v.x() * v.x() + v.y() * v.y() + v.z() * v.z(); };
complex_t beta = std::sqrt(square(b_p_vec)); complex_t beta = std::sqrt(square(b));
if (std::abs(beta) < std::numeric_limits<double>::epsilon() * 10.) { if (std::abs(beta) < std::numeric_limits<double>::epsilon() * 10.) {
const complex_t alpha_pp = -(alpha * alpha) * sigma * sigma / 8.; const complex_t alpha_pp = -(alpha * alpha) * sigma * sigma / 8.;
return SpinMatrix(std::exp(alpha_pp), 0, 0, std::exp(alpha_pp)); return SpinMatrix(std::exp(alpha_pp), 0, 0, std::exp(alpha_pp));
} }
b_p_vec /= beta; b /= beta;
const complex_t alpha_pp = -(alpha * alpha + beta * beta) * sigma * sigma / 8.; const complex_t alpha_pp = -(alpha * alpha + beta * beta) * sigma * sigma / 8.;
const complex_t beta_pp = -alpha * beta * sigma * sigma / 4.; const complex_t beta_pp = -alpha * beta * sigma * sigma / 4.;
SpinMatrix Q(b.z() + 1., b.x() - I * b.y(), b.x() + I * b.y(), -1. - b.z());
const SpinMatrix M(std::exp(beta_pp), 0, 0, std::exp(-beta_pp));
const complex_t factor1 = std::sqrt(2. * (1. + b_p_vec.z())); return std::exp(alpha_pp) * Q * M * Q.adjoint() / (2. * (1. + b.z()));
const complex_t factor2 = std::sqrt(2. * (1. - b_p_vec.z()));
SpinMatrix Q((b_p_vec.z() + 1.) / factor1, (b_p_vec.z() - 1.) / factor2,
(b_p_vec.x() + I * b_p_vec.y()) / factor1,
(b_p_vec.x() + I * b_p_vec.y()) / factor2);
const SpinMatrix exp1(std::exp(alpha_pp), 0, 0, std::exp(alpha_pp));
const SpinMatrix exp2(std::exp(beta_pp), 0, 0, std::exp(-beta_pp));
return exp1 * Q * exp2 * Q.adjoint();
}; };
const auto kk = SpinMatrix(coeff_a.computeInverseKappa() * coeff_b.computeKappa()); const auto kk = SpinMatrix(coeff_a.computeInverseKappa() * coeff_b.computeKappa());
......
No preview for this file type
No preview for this file type
No preview for this file type
No preview for this file type
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment