TransitionMagneticNevot.cpp

Go to the documentation of this file.

//  ************************************************************************************************
//
//  BornAgain: simulate and fit reflection and scattering
//
//! @file      Resample/Specular/TransitionMagneticNevot.cpp
//! @brief     Implements namespace MagneticNevotCroceTransition.
//!
//! @homepage  http://www.bornagainproject.org
//! @license   GNU General Public License v3 or higher (see COPYING)
//! @copyright Forschungszentrum Jülich GmbH 2020
//! @authors   Scientific Computing Group at MLZ (see CITATION, AUTHORS)
//
//  ************************************************************************************************
 
#include "Resample/Specular/TransitionMagneticNevot.h"
#include "Resample/Flux/MatrixFlux.h"
 
namespace {
 
complex_t checkForUnderflow(complex_t val)
{
    return std::abs(val.imag()) < 1e-80 && val.real() < 0 ? complex_t(val.real(), 1e-40) : val;
}
 
std::pair<SpinMatrix, SpinMatrix> computeRoughnessMatrices(const MatrixFlux& coeff_i,
                                                           const MatrixFlux& coeff_i1, double sigma)
{
    complex_t beta_i = coeff_i.m_lambda.v - coeff_i.m_lambda.u;
    complex_t beta_i1 = coeff_i1.m_lambda.v - coeff_i1.m_lambda.u;
 
    auto roughness_matrix = [sigma, &coeff_i, &coeff_i1, beta_i, beta_i1](double sign) {
        const complex_t alpha_p = coeff_i1.m_lambda.u + coeff_i1.m_lambda.v
                                  + sign * (coeff_i.m_lambda.u + coeff_i.m_lambda.v);
        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(); };
        complex_t beta_p = std::sqrt(checkForUnderflow(square(b_p_vec)));
        b_p_vec /= beta_p;
 
        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.;
 
        const complex_t factor1 = std::sqrt(2. * (1. + b_p_vec.z()));
        const complex_t factor2 = std::sqrt(2. * (1. - b_p_vec.z()));
        SpinMatrix 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);
        SpinMatrix 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);
 
        const SpinMatrix exp1(std::exp(alpha_pp), 0, 0, std::exp(alpha_pp));
 
        if (std::abs(beta_p) > std::numeric_limits<double>::epsilon() * 10.) {
            SpinMatrix exp2(std::exp(beta_pp), 0, 0, std::exp(-beta_pp));
            return SpinMatrix{exp1 * QL * exp2 * QR};
        }
 
        return exp1;
    };
 
    const SpinMatrix roughness_sum = roughness_matrix(1.);
    const SpinMatrix roughness_diff = roughness_matrix(-1.);
 
    return {roughness_sum, roughness_diff};
}
 
} // namespace
 
 
std::pair<SpinMatrix, SpinMatrix> Compute::MagneticNevotCroceTransition::backwardsSubmatrices(
    const MatrixFlux& coeff_i, const MatrixFlux& coeff_i1, double sigma)
{
    SpinMatrix roughness_sum = SpinMatrix::One();
    SpinMatrix roughness_diff = SpinMatrix::One();
    if (sigma != 0.) {
        const auto ret = computeRoughnessMatrices(coeff_i, coeff_i1, sigma);
        roughness_sum = std::get<0>(ret);
        roughness_diff = std::get<1>(ret);
    }
 
    const auto P = SpinMatrix(coeff_i.computeInverseP() * coeff_i1.computeP());
    const auto mp = 0.5 * (SpinMatrix::One() + P) * roughness_diff;
    const auto mm = 0.5 * (SpinMatrix::One() - P) * roughness_sum;
 
    return {mp, mm};
}