ComputeFluxMagnetic.cpp

Go to the documentation of this file.

//  ************************************************************************************************
//
//  BornAgain: simulate and fit reflection and scattering
//
//! @file      Resample/Specular/ComputeFluxMagnetic.cpp
//! @brief     Implements class SpecularMagneticStrategy.
//!
//! @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/ComputeFluxMagnetic.h"
#include "Base/Const/PhysicalConstants.h"
#include "Base/Util/Assert.h"
#include "Resample/Flux/MatrixFlux.h"
#include "Resample/Slice/KzComputation.h"
#include "Resample/Slice/SliceStack.h"
#include "Resample/Specular/TransitionMagneticNevot.h"
#include "Resample/Specular/TransitionMagneticTanh.h"
#include "Sample/Interface/LayerRoughness.h"
#include <algorithm>
 
namespace {
 
// The factor 1e-18 is here to have unit: 1/T*nm^-2
constexpr double magnetic_prefactor = PhysConsts::m_n * PhysConsts::g_factor_n * PhysConsts::mu_N
                                      / PhysConsts::h_bar / PhysConsts::h_bar / 4. / M_PI * 1e-18;
const auto eps = std::numeric_limits<double>::epsilon() * 10.;
 
double magneticSLD(R3 B_field)
{
    return magnetic_prefactor * B_field.mag();
}
 
Spinor eigenvalues(complex_t kz, double magnetic_SLD)
{
    const complex_t a = kz * kz;
    return {std::sqrt(a - 4. * M_PI * magnetic_SLD), std::sqrt(a + 4. * M_PI * magnetic_SLD)};
}
 
Spinor checkForUnderflow(const Spinor& eigenvs)
{
    auto lambda = [](complex_t value) { return std::abs(value) < 1e-40 ? 1e-40 : value; };
    return {lambda(eigenvs.u), lambda(eigenvs.v)};
}
 
std::pair<SpinMatrix, SpinMatrix> computeBackwardsSubmatrices(const MatrixFlux& coeff_i,
                                                              const MatrixFlux& coeff_i1,
                                                              double sigma,
                                                              const RoughnessModel& r_model)
{
    if (r_model == RoughnessModel::NEVOT_CROCE)
        return Compute::MagneticNevotCroceTransition::backwardsSubmatrices(coeff_i, coeff_i1,
                                                                           sigma);
    return Compute::TransitionMagneticTanh::backwardsSubmatrices(coeff_i, coeff_i1, sigma);
}
 
void calculateUpwards(std::vector<MatrixFlux>& coeff, const SliceStack& slices,
                      const RoughnessModel& r_model)
{
    const auto N = slices.size();
    std::vector<SpinMatrix> SMatrices(N - 1);
    std::vector<complex_t> Normalization(N - 1);
 
    // bottom boundary condition
    coeff.back().m_T = SpinMatrix::One();
    coeff.back().m_R = SpinMatrix();
 
    for (size_t i = N - 1; i > 0; --i) {
        double sigma = 0.;
        if (const auto* const roughness = slices.bottomRoughness(i))
            sigma = roughness->sigma();
 
        // compute the 2x2 submatrices in the write-up denoted as P+, P- and delta
        const auto [mp, mm] = computeBackwardsSubmatrices(coeff[i - 1], coeff[i], sigma, r_model);
 
        const SpinMatrix delta = coeff[i - 1].computeDeltaMatrix(slices[i - 1].thicknessOr0());
 
        // compute the rotation matrix
        SpinMatrix Si = mp + mm * coeff[i].m_R;
        SpinMatrix S(Si.d, -Si.b, -Si.c, Si.a);
        const complex_t norm = S.determinant();
        S = S * delta;
 
        // store the rotation matrix and normalization constant in order to rotate
        // the coefficients for all lower slices at the end of the computation
        SMatrices[i - 1] = S;
        Normalization[i - 1] = norm;
 
        // compute the reflection matrix and
        // rotate the polarization such that we have pure incoming states (T = I)
        S /= norm;
 
        // T is always equal to the identity at this point, no need to store
        coeff[i - 1].m_R = delta * (mm + mp * coeff[i].m_R) * S;
    }
 
    // now correct all amplitudes in forward direction by dividing with the remaining
    // 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
    complex_t dumpingFactor = 1;
    SpinMatrix S = SpinMatrix::One();
    for (size_t i = 1; i < N; ++i) {
        dumpingFactor = dumpingFactor * Normalization[i - 1];
        S = SMatrices[i - 1] * S;
 
        if (std::isinf(std::norm(dumpingFactor))) {
            std::for_each(coeff.begin() + i, coeff.end(), [](auto& coeff) {
                coeff.m_T = SpinMatrix();
                coeff.m_R = SpinMatrix();
            });
            break;
        }
 
        coeff[i].m_T = S / dumpingFactor; // T * S omitted, since T is always I
        coeff[i].m_R *= S / dumpingFactor;
    }
}
 
std::vector<MatrixFlux> computeFlux(const SliceStack& slices, const std::vector<complex_t>& kzs,
                                    const RoughnessModel& r_model, bool forward)
{
    const size_t N = slices.size();
 
    if (slices.size() != kzs.size())
        throw std::runtime_error(
            "Error in SpecularMagnetic_::execute: kz vector and slices size shall coinside.");
    if (slices.empty())
        return {};
 
    std::vector<MatrixFlux> result;
    result.reserve(N);
 
    const double kz_sign = kzs.front().real() >= 0.0 ? 1.0 : -1.0; // save sign to restore it later
 
    auto B_0 = slices.front().bField();
    result.emplace_back(kz_sign, eigenvalues(kzs.front(), 0.0), R3{0.0, 0.0, 0.0}, 0.0);
    for (size_t i = 1; i < slices.size(); ++i) {
        auto B = (forward ? +1 : -1) * (slices[i].bField() - B_0);
        auto magnetic_SLD = magneticSLD(B);
        result.emplace_back(kz_sign, checkForUnderflow(eigenvalues(kzs[i], magnetic_SLD)),
                            B.mag() > eps ? B / B.mag() : R3{0.0, 0.0, 0.0}, magnetic_SLD);
    }
 
    if (N == 1) {
        result[0].m_T = SpinMatrix::One();
        result[0].m_R = SpinMatrix();
        return result;
    }
 
    if (kzs[0] == 0.) {
        result[0].m_T = SpinMatrix::One();
        result[0].m_R = SpinMatrix::Diag(-1, -1);
        for (size_t i = 1; i < N; ++i) {
            result[i].m_T = SpinMatrix();
            result[i].m_R = SpinMatrix();
        }
        return result;
    }
 
    calculateUpwards(result, slices, r_model);
 
    return result;
}
 
} // namespace
 
//  ************************************************************************************************
//  namespace implementation
//  ************************************************************************************************
 
Fluxes Compute::SpecularMagnetic::fluxes(const SliceStack& slices, const R3& k, bool forward)
{
    if (slices.size() > 1 && k.z() > 0)
        throw std::runtime_error(
            "source or detector below horizon not yet implemented for polarized scattering");
    std::vector<complex_t> kz = Compute::Kz::computeReducedKz(slices, k);
    ASSERT(slices.size() == kz.size());
 
    Fluxes result;
    for (auto& coeff : computeFlux(slices, kz, slices.roughnessModel(), forward))
        result.emplace_back(std::make_unique<const MatrixFlux>(coeff));
 
    return result;
}
 
SpinMatrix Compute::SpecularMagnetic::topLayerR(const SliceStack& slices,
                                                const std::vector<complex_t>& kzs, bool forward)
{
    const RoughnessModel& r_model = slices.roughnessModel();
 
    ASSERT(slices.size() == kzs.size());
    const auto N = slices.size();
    if (N == 1)
        return SpinMatrix();
    if (kzs[0] == 0.)
        return -SpinMatrix::One();
 
    auto B_0 = slices.front().bField();
    const double kz_sign = kzs.front().real() >= 0.0 ? 1.0 : -1.0; // save sign to restore it later
 
    auto createCoeff = [&slices, &kzs, kz_sign, B_0, forward](size_t i) {
        const auto B = (forward ? +1 : -1) * (slices[i].bField() - B_0);
        const auto magnetic_SLD = magneticSLD(B);
 
        return MatrixFlux(kz_sign, checkForUnderflow(eigenvalues(kzs[i], magnetic_SLD)),
                          B.mag() > eps ? B / B.mag() : R3{0.0, 0.0, 0.0}, magnetic_SLD);
    };
 
    auto c_i1 = createCoeff(N - 1);
 
    // bottom boundary condition
    c_i1.m_R = SpinMatrix();
 
    for (size_t i = N - 1; i > 0; --i) {
        auto c_i = createCoeff(i - 1);
 
        double sigma = 0.;
        if (const auto* const roughness = slices.bottomRoughness(i - 1))
            sigma = roughness->sigma();
 
        // compute the 2x2 submatrices in the write-up denoted as P+, P- and delta
        const auto [mp, mm] = computeBackwardsSubmatrices(c_i, c_i1, sigma, r_model);
 
        const SpinMatrix delta = c_i.computeDeltaMatrix(slices[i - 1].thicknessOr0());
 
        // compute the rotation matrix
        SpinMatrix Si = mp + mm * c_i1.m_R;
        SpinMatrix S(Si.d, -Si.b, -Si.c, Si.a);
        const complex_t norm = S.determinant();
        S = S * (delta / norm);
 
        c_i.m_R = delta * (mm + mp * c_i1.m_R) * S;
        c_i1 = c_i;
    }
    return c_i1.m_R;
}