MatrixFlux.cpp

Go to the documentation of this file.

//  ************************************************************************************************
//
//  BornAgain: simulate and fit reflection and scattering
//
//! @file     Resample/Flux/MatrixFlux.cpp
//! @brief    Implements class MatrixFlux.
//!
//! @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/Flux/MatrixFlux.h"
#include "Base/Util/Assert.h"
 
namespace {
 
const auto eps = std::numeric_limits<double>::epsilon() * 10.;
 
complex_t GetImExponential(complex_t exponent)
{
    if (exponent.imag() > -std::log(std::numeric_limits<double>::min()))
        return 0.0;
    return std::exp(I * exponent);
}
 
} // namespace
 
 
MatrixFlux::MatrixFlux(double kz_sign, const Spinor& eigenvalues, const R3& b, double magnetic_SLD)
    : m_lambda(std::move(eigenvalues))
    , m_T(+1, 0, 0, +1)
    , m_R(-1, 0, 0, -1)
    , m_b(b)
    , m_kz_sign(kz_sign)
    , m_magnetic_SLD(magnetic_SLD)
{
    ASSERT(std::abs(m_b.mag() - 1) < eps || (m_b.mag() < eps && magnetic_SLD < eps));
}
 
SpinMatrix MatrixFlux::TransformationMatrix(const Spinor& selection) const
{
    const SpinMatrix exp2(selection.u, 0, 0, selection.v);
 
    if (std::abs(m_b.mag() - 1.) < eps) {
        const double factor1 = 2. * (1. + m_b.z());
        SpinMatrix Q(1. + m_b.z(), I * m_b.y() - m_b.x(), m_b.x() + I * m_b.y(), m_b.z() + 1.);
        return Q * exp2 * Q.adjoint() / factor1;
    }
    if (m_b.mag() < eps)
        return exp2;
 
    throw std::runtime_error("Broken magnetic field vector");
}
 
SpinMatrix MatrixFlux::T1Matrix() const
{
    return TransformationMatrix({0., 1.});
}
 
SpinMatrix MatrixFlux::T2Matrix() const
{
    return TransformationMatrix({1., 0.});
}
 
Spinor MatrixFlux::T1plus() const
{
    return T1Matrix() * m_T.col0();
}
 
Spinor MatrixFlux::R1plus() const
{
    return T1Matrix() * m_R.col0();
}
 
Spinor MatrixFlux::T2plus() const
{
    return T2Matrix() * m_T.col0();
}
 
Spinor MatrixFlux::R2plus() const
{
    return T2Matrix() * m_R.col0();
}
 
Spinor MatrixFlux::T1min() const
{
    return T1Matrix() * m_T.col1();
}
 
Spinor MatrixFlux::R1min() const
{
    return T1Matrix() * m_R.col1();
}
 
Spinor MatrixFlux::T2min() const
{
    return T2Matrix() * m_T.col1();
}
 
Spinor MatrixFlux::R2min() const
{
    return T2Matrix() * m_R.col1();
}
 
Spinor MatrixFlux::getKz() const
{
    return m_kz_sign * m_lambda;
}
 
SpinMatrix MatrixFlux::pMatrixHelper(double sign) const
{
    const complex_t alpha = m_lambda.v + m_lambda.u;
    const complex_t beta = m_lambda.v - m_lambda.u;
 
    SpinMatrix result(alpha + sign * beta * m_b.z(), sign * beta * (m_b.x() - I * m_b.y()),
                      sign * beta * (m_b.x() + I * m_b.y()), alpha - sign * beta * m_b.z());
 
    return m_kz_sign * result;
}
 
SpinMatrix MatrixFlux::computeP() const
{
    SpinMatrix result = pMatrixHelper(1.);
    result *= 0.5;
 
    return result;
}
 
SpinMatrix MatrixFlux::computeInverseP() const
{
    const complex_t alpha = m_lambda.v + m_lambda.u;
    const complex_t beta = m_lambda.v - m_lambda.u;
 
    if (std::abs(alpha * alpha - beta * beta) == 0.)
        return SpinMatrix();
 
    SpinMatrix result = pMatrixHelper(-1.);
    result *= 2. / (alpha * alpha - beta * beta);
 
    return result;
}
 
SpinMatrix MatrixFlux::computeDeltaMatrix(double thickness)
{
    const complex_t alpha = 0.5 * thickness * (m_lambda.v + m_lambda.u);
 
    // 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) {
        const SpinMatrix exp2(GetImExponential(m_kz_sign * thickness * m_lambda.v), 0, 0,
                              GetImExponential(m_kz_sign * thickness * m_lambda.u));
        const double factor1 = 2. * (1. + m_b.z());
        SpinMatrix Q((1. + m_b.z()), (I * m_b.y() - m_b.x()), (m_b.x() + I * m_b.y()),
                     (m_b.z() + 1.));
 
        return Q * exp2 * Q.adjoint() / factor1;
    }
    if (m_b.mag() < eps)
        return SpinMatrix::One() * GetImExponential(m_kz_sign * alpha);
 
    throw std::runtime_error("Broken magnetic field vector");
}