diff --git a/Core/Multilayer/SpecularMagneticNewNCStrategy.cpp b/Core/Multilayer/SpecularMagneticNewNCStrategy.cpp index da03ec2b0ecd6143883ef4c3e0b94dee8ca37b8e..426fc0c91ffcdb1b3a51b2ebea4e8b895599d1da 100644 --- a/Core/Multilayer/SpecularMagneticNewNCStrategy.cpp +++ b/Core/Multilayer/SpecularMagneticNewNCStrategy.cpp @@ -24,51 +24,49 @@ SpecularMagneticNewNCStrategy::computeRoughnessMatrices(const MatrixRTCoefficien const MatrixRTCoefficients_v3& coeff_i1, double sigma) const { - auto beta_i = coeff_i.m_lambda(1) - coeff_i.m_lambda(0); - auto beta_i1 = coeff_i1.m_lambda(1) - coeff_i1.m_lambda(0); + complex_t beta_i = coeff_i.m_lambda(1) - coeff_i.m_lambda(0); + complex_t beta_i1 = coeff_i1.m_lambda(1) - coeff_i1.m_lambda(0); auto roughness_matrix = [sigma, &coeff_i, &coeff_i1, beta_i, beta_i1](double sign){ - auto alpha_p = coeff_i1.m_lambda(0) + coeff_i1.m_lambda(1) + + const complex_t alpha_p = coeff_i1.m_lambda(0) + coeff_i1.m_lambda(1) + sign * (coeff_i.m_lambda(0) + coeff_i.m_lambda(1)); auto b_p_vec = beta_i1 * coeff_i1.m_b + sign * beta_i * coeff_i.m_b; auto square = [](auto & v){ return v.x() * v.x() + v.y() * v.y() + v.z() * v.z(); }; - auto beta_p = std::sqrt(checkForUnderflow(square(b_p_vec))); + complex_t beta_p = std::sqrt(checkForUnderflow(square(b_p_vec))); b_p_vec /= beta_p; - auto alpha_pp = -(alpha_p * alpha_p + beta_p * beta_p) * sigma * sigma / 8.; - auto beta_pp = - alpha_p * beta_p * sigma * sigma / 4.; + const complex_t alpha_pp = -(alpha_p * alpha_p + beta_p * beta_p) * sigma * sigma / 8.; + const complex_t beta_pp = - alpha_p * beta_p * sigma * sigma / 4.; Eigen::Matrix2cd QL, QR; - auto factor1 = std::sqrt(2. * (1. + b_p_vec.z())); - auto factor2 = std::sqrt(2. * (1. - b_p_vec.z())); + const complex_t factor1 = std::sqrt(2. * (1. + b_p_vec.z())); + const complex_t factor2 = std::sqrt(2. * (1. - b_p_vec.z())); QL << (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; QR << (b_p_vec.z() + 1.)/factor1, (b_p_vec.x() - I * b_p_vec.y())/factor1, (b_p_vec.z() - 1.)/factor2, (b_p_vec.x() - I * b_p_vec.y())/factor2; - auto exp1 = Eigen::Matrix2cd(Eigen::DiagonalMatrix<complex_t, 2>( - {std::exp(alpha_pp), std::exp(alpha_pp)})); + const Eigen::Matrix2cd exp1 = Eigen::DiagonalMatrix<complex_t, 2>( + {std::exp(alpha_pp), std::exp(alpha_pp)}); - Eigen::Matrix2cd roughnessMatrix; if( std::abs(beta_p) > std::numeric_limits<double>::epsilon() * 10. ) { Eigen::Matrix2cd exp2 = Eigen::Matrix2cd(Eigen::DiagonalMatrix<complex_t, 2>( {std::exp(beta_pp), std::exp(-beta_pp)})); - roughnessMatrix = exp1 * QL * exp2 * QR; + return Eigen::Matrix2cd{exp1 * QL * exp2 * QR}; + } - }else - roughnessMatrix = exp1; - return roughnessMatrix; + return exp1; }; - auto roughness_sum = roughness_matrix(1.); - auto roughness_diff = roughness_matrix(-1.); + const Eigen::Matrix2cd roughness_sum = roughness_matrix(1.); + const Eigen::Matrix2cd roughness_diff = roughness_matrix(-1.); return {roughness_sum, roughness_diff}; } @@ -81,14 +79,14 @@ SpecularMagneticNewNCStrategy::computeBackwardsSubmatrices(const MatrixRTCoeffic Eigen::Matrix2cd roughness_sum{Eigen::Matrix2cd::Identity()}; Eigen::Matrix2cd roughness_diff{Eigen::Matrix2cd::Identity()}; if (sigma != 0.) { - auto ret = computeRoughnessMatrices(coeff_i, coeff_i1, sigma); + const auto ret = computeRoughnessMatrices(coeff_i, coeff_i1, sigma); roughness_sum = std::get<0>(ret); roughness_diff = std::get<1>(ret); } - auto P = Eigen::Matrix2cd(coeff_i.computeInverseP() * coeff_i1.computeP()); - auto mp = 0.5 * (Eigen::Matrix2cd::Identity() + P) * roughness_diff; - auto mm = 0.5 * (Eigen::Matrix2cd::Identity() - P) * roughness_sum; + const auto P = Eigen::Matrix2cd(coeff_i.computeInverseP() * coeff_i1.computeP()); + const auto mp = 0.5 * (Eigen::Matrix2cd::Identity() + P) * roughness_diff; + const auto mm = 0.5 * (Eigen::Matrix2cd::Identity() - P) * roughness_sum; return {mp, mm}; } diff --git a/Core/Multilayer/SpecularMagneticNewStrategy.cpp b/Core/Multilayer/SpecularMagneticNewStrategy.cpp index 53d41f45cf0f5ed3a877b414976dbc8af7ebe1a0..5462c6643906ff369098079cdbb6e8d7e6e60675 100644 --- a/Core/Multilayer/SpecularMagneticNewStrategy.cpp +++ b/Core/Multilayer/SpecularMagneticNewStrategy.cpp @@ -97,7 +97,7 @@ SpecularMagneticNewStrategy::computeTR(const std::vector<Slice>& slices, void SpecularMagneticNewStrategy::calculateUpwards( std::vector<MatrixRTCoefficients_v3>& coeff, const std::vector<Slice>& slices) const { - auto N = slices.size(); + const auto N = slices.size(); std::vector<Eigen::Matrix2cd> SMatrices(N-1); std::vector<complex_t> Normalization(N-1); @@ -111,7 +111,7 @@ void SpecularMagneticNewStrategy::calculateUpwards( sigma = roughness->getSigma(); // compute the 2x2 submatrices in the write-up denoted as P+, P- and delta - auto mpmm = computeBackwardsSubmatrices(coeff[i], coeff[i + 1], sigma); + const auto mpmm = computeBackwardsSubmatrices(coeff[i], coeff[i + 1], sigma); const Eigen::Matrix2cd mp = mpmm.first; const Eigen::Matrix2cd mm = mpmm.second; @@ -122,7 +122,7 @@ void SpecularMagneticNewStrategy::calculateUpwards( Si = mp + mm * coeff[i+1].m_R; S << Si(1, 1) , -Si(0, 1), -Si(1, 0), Si(0, 0); - auto norm = S(0, 0) * S(1, 1) - S(0, 1) * S(1, 0); + const complex_t norm = S(0, 0) * S(1, 1) - S(0, 1) * S(1, 0); S = S * delta; // store the rotation matrix and normalization constant in order to rotate @@ -142,7 +142,7 @@ void SpecularMagneticNewStrategy::calculateUpwards( // normalization constants. In addition rotate the polarization by the amount // that was rotated above the current interface // if the normalization overflows, all amplitudes below that point are set to zero - auto dumpingFactor = complex_t(1, 0); + complex_t dumpingFactor = 1; Eigen::Matrix2cd S = Eigen::Matrix2cd::Identity(); for (size_t i = 1; i < N; ++i) { dumpingFactor = dumpingFactor * Normalization[i - 1]; diff --git a/Core/Multilayer/SpecularMagneticNewTanhStrategy.cpp b/Core/Multilayer/SpecularMagneticNewTanhStrategy.cpp index 4725681b3321799c794051e2fc4b8f129fb9383b..7ca12455e59b85d4f940dbef5e382f672fbbd354 100644 --- a/Core/Multilayer/SpecularMagneticNewTanhStrategy.cpp +++ b/Core/Multilayer/SpecularMagneticNewTanhStrategy.cpp @@ -28,41 +28,38 @@ SpecularMagneticNewTanhStrategy::computeRoughnessMatrix(const MatrixRTCoefficien if (sigma < 10 * std::numeric_limits<double>::epsilon()) return Eigen::Matrix2cd{Eigen::Matrix2cd::Identity()}; - const auto sigeff = pi2_15 * sigma; - auto b = coeff.m_b; + const double sigeff = pi2_15 * sigma; + const auto b = coeff.m_b; - Eigen::Matrix2cd R; if (std::abs(b.mag() - 1.) < std::numeric_limits<double>::epsilon() * 10.) { Eigen::Matrix2cd Q; - auto factor1 = 2. * (1. + b.z()); + const double factor1 = 2. * (1. + b.z()); Q << (1. + b.z()), (I * b.y() - b.x()), (b.x() + I * b.y()), (b.z() + 1.); - auto l1 = std::sqrt(MathFunctions::tanhc(sigeff * coeff.m_lambda(1))); - auto l2 = std::sqrt(MathFunctions::tanhc(sigeff * coeff.m_lambda(0))); + complex_t l1 = std::sqrt(MathFunctions::tanhc(sigeff * coeff.m_lambda(1))); + complex_t l2 = std::sqrt(MathFunctions::tanhc(sigeff * coeff.m_lambda(0))); if (inverse) { l1 = 1. / l1; l2 = 1. / l2; } - auto lambda = Eigen::Matrix2cd(Eigen::DiagonalMatrix<complex_t, 2>({l1, l2})); + const Eigen::Matrix2cd lambda = Eigen::DiagonalMatrix<complex_t, 2>({l1, l2}); - R = Q * lambda * Q.adjoint() / factor1; + return Q * lambda * Q.adjoint() / factor1; } else if (b.mag() < 10 * std::numeric_limits<double>::epsilon()) { - auto alpha = + complex_t alpha = std::sqrt(MathFunctions::tanhc(0.5 * sigeff * (coeff.m_lambda(1) + coeff.m_lambda(0)))); if (inverse) alpha = 1. / alpha; - auto lambda = Eigen::Matrix2cd(Eigen::DiagonalMatrix<complex_t, 2>({alpha, alpha})); + const Eigen::Matrix2cd lambda = Eigen::DiagonalMatrix<complex_t, 2>({alpha, alpha}); - R = lambda; - - } else - throw std::runtime_error("Broken magnetic field vector"); + return lambda; + } - return R; + throw std::runtime_error("Broken magnetic field vector"); } std::pair<Eigen::Matrix2cd, Eigen::Matrix2cd> @@ -81,9 +78,9 @@ SpecularMagneticNewTanhStrategy::computeBackwardsSubmatrices( * computeRoughnessMatrix(coeff_i1, sigma, true); } - auto mproduct = Eigen::Matrix2cd(coeff_i.computeInverseP() * coeff_i1.computeP()); - auto mp = 0.5 * (RInv + mproduct * R); - auto mm = 0.5 * (RInv - mproduct * R); + const Eigen::Matrix2cd mproduct = coeff_i.computeInverseP() * coeff_i1.computeP(); + const Eigen::Matrix2cd mp = 0.5 * (RInv + mproduct * R); + const Eigen::Matrix2cd mm = 0.5 * (RInv - mproduct * R); return {mp, mm}; } diff --git a/Core/RT/MatrixRTCoefficients_v3.cpp b/Core/RT/MatrixRTCoefficients_v3.cpp index 32e2e74b8c896fc826b63d45f293e4a50c12b52a..18b4ec09c8c749ee42ecd529e209429a698104ec 100644 --- a/Core/RT/MatrixRTCoefficients_v3.cpp +++ b/Core/RT/MatrixRTCoefficients_v3.cpp @@ -14,7 +14,6 @@ #include "Core/RT/MatrixRTCoefficients_v3.h" #include "Core/Basics/Assert.h" -#include <iostream> namespace { @@ -45,25 +44,21 @@ MatrixRTCoefficients_v3* MatrixRTCoefficients_v3::clone() const Eigen::Matrix2cd MatrixRTCoefficients_v3::TransformationMatrix(complex_t eigenvalue, Eigen::Vector2d selection) const { - Eigen::Matrix2cd result; - - auto exp2 = Eigen::Matrix2cd(Eigen::DiagonalMatrix<complex_t, 2>(selection)); + const Eigen::Matrix2cd exp2 = Eigen::DiagonalMatrix<complex_t, 2>(selection); if (std::abs(m_b.mag() - 1.) < eps) { Eigen::Matrix2cd Q; - auto factor1 = 2. * (1. + m_b.z()); + const double factor1 = 2. * (1. + m_b.z()); Q << (1. + m_b.z()), (I * m_b.y() - m_b.x()), (m_b.x() + I * m_b.y()), (m_b.z() + 1.); - result = Q * exp2 * Q.adjoint() / factor1; + return Q * exp2 * Q.adjoint() / factor1; } else if (m_b.mag() < eps && eigenvalue != 0.) - result = exp2; + return exp2; // result = Eigen::Matrix2cd(Eigen::DiagonalMatrix<complex_t, 2>({0.5, 0.5})); else if (m_b.mag() < eps && eigenvalue == 0.) - result = Eigen::Matrix2cd(Eigen::DiagonalMatrix<complex_t, 2>({0.5, 0.5})); - else - throw std::runtime_error("Broken magnetic field vector"); + return Eigen::Matrix2cd(Eigen::DiagonalMatrix<complex_t, 2>({0.5, 0.5})); - return result; + throw std::runtime_error("Broken magnetic field vector"); } Eigen::Matrix2cd MatrixRTCoefficients_v3::T1Matrix() const @@ -78,52 +73,42 @@ Eigen::Matrix2cd MatrixRTCoefficients_v3::T2Matrix() const Eigen::Vector2cd MatrixRTCoefficients_v3::T1plus() const { - auto mat = T1Matrix(); - return mat * m_T.col(0); - ; + return T1Matrix() * m_T.col(0); } Eigen::Vector2cd MatrixRTCoefficients_v3::R1plus() const { - auto mat = T1Matrix(); - return mat * m_R.col(0); - ; + return T1Matrix() * m_R.col(0); } Eigen::Vector2cd MatrixRTCoefficients_v3::T2plus() const { - auto mat = T2Matrix(); - return mat * m_T.col(0); + return T2Matrix() * m_T.col(0); } Eigen::Vector2cd MatrixRTCoefficients_v3::R2plus() const { - auto mat = T2Matrix(); - return mat * m_R.col(0); + return T2Matrix() * m_R.col(0); } Eigen::Vector2cd MatrixRTCoefficients_v3::T1min() const { - auto mat = T1Matrix(); - return mat * m_T.col(1); + return T1Matrix() * m_T.col(1); } Eigen::Vector2cd MatrixRTCoefficients_v3::R1min() const { - auto mat = T1Matrix(); - return mat * m_R.col(1); + return T1Matrix() * m_R.col(1); } Eigen::Vector2cd MatrixRTCoefficients_v3::T2min() const { - auto mat = T2Matrix(); - return mat * m_T.col(1); + return T2Matrix() * m_T.col(1); } Eigen::Vector2cd MatrixRTCoefficients_v3::R2min() const { - auto mat = T2Matrix(); - return mat * m_R.col(1); + return T2Matrix() * m_R.col(1); } Eigen::Vector2cd MatrixRTCoefficients_v3::getKz() const @@ -171,10 +156,10 @@ Eigen::Matrix2cd MatrixRTCoefficients_v3::computeInverseP() const Eigen::Matrix2cd MatrixRTCoefficients_v3::computeDeltaMatrix(double thickness) { Eigen::Matrix2cd result; - const auto alpha = 0.5 * thickness * (m_lambda(1) + m_lambda(0)); + const complex_t alpha = 0.5 * thickness * (m_lambda(1) + m_lambda(0)); - const auto exp2 = Eigen::Matrix2cd(Eigen::DiagonalMatrix<complex_t, 2>( - {GetImExponential(thickness * m_lambda(1)), GetImExponential(thickness * m_lambda(0))})); + const Eigen::Matrix2cd exp2 = Eigen::DiagonalMatrix<complex_t, 2>( + {GetImExponential(thickness * m_lambda(1)), GetImExponential(thickness * m_lambda(0))}); // Compute resulting phase matrix according to exp(i p_m d_m) = exp1 * Q * exp2 * Q.adjoint(); if (std::abs(m_b.mag() - 1.) < eps) { @@ -182,15 +167,12 @@ Eigen::Matrix2cd MatrixRTCoefficients_v3::computeDeltaMatrix(double thickness) const double factor1 = 2. * (1. + m_b.z()); Q << (1. + m_b.z()), (I * m_b.y() - m_b.x()), (m_b.x() + I * m_b.y()), (m_b.z() + 1.); - result = Q * exp2 * Q.adjoint() / factor1; - - } else if (m_b.mag() < eps) { - result = Eigen::Matrix2cd::Identity() * GetImExponential(alpha); + return Q * exp2 * Q.adjoint() / factor1; - } else - throw std::runtime_error("Broken magnetic field vector"); + } else if (m_b.mag() < eps) + return Eigen::Matrix2cd::Identity() * GetImExponential(alpha); - return result; + throw std::runtime_error("Broken magnetic field vector"); } namespace