SumDWBA.cpp

//  ************************************************************************************************
//
//  BornAgain: simulate and fit reflection and scattering
//
//! @file      Resample/Coherence/SumDWBA.cpp
//! @brief     Implements class SumDWBA.
//!
//! @homepage  http://www.bornagainproject.org
//! @license   GNU General Public License v3 or higher (see COPYING)
//! @copyright Forschungszentrum Jülich GmbH 2018
//! @authors   Scientific Computing Group at MLZ (see CITATION, AUTHORS)
//
//  ************************************************************************************************

#include "Resample/Coherence/SumDWBA.h"
#include "Base/Util/Assert.h"
#include "Base/Vector/WavevectorInfo.h"
#include "Resample/Element/DiffuseElement.h"
#include "Resample/Flux/MatrixFlux.h"
#include "Resample/Flux/ScalarFlux.h"
#include "Sample/Scattering/IFormFactor.h"
#include <iostream>

namespace {

complex_t VecMatVecProduct(const Eigen::Vector2cd& vec1, const Eigen::Matrix2cd& ff,
                           const Eigen::Vector2cd& vec2)
{
    return vec1.transpose() * ff * vec2;
}

} // namespace


SumDWBA::SumDWBA(const IFormFactor& ff, size_t i_layer) : m_ff(ff.clone()), m_i_layer(i_layer) {}

SumDWBA::SumDWBA(const IFormFactor& ff) : m_ff(ff.clone()), m_i_layer() {}

SumDWBA::~SumDWBA() = default;

size_t SumDWBA::iLayer() const
{
    ASSERT(m_i_layer.has_value());
    return m_i_layer.value();
}

complex_t SumDWBA::coherentFF(const DiffuseElement& ele) const
{
    const WavevectorInfo& wavevectors = ele.wavevectorInfo();

    if (!m_i_layer.has_value())
        // no slicing, pure Born approximation
        return m_ff->theFF(wavevectors);

    const auto* inFlux = dynamic_cast<const ScalarFlux*>(ele.fluxIn(iLayer()));
    const auto* outFlux = dynamic_cast<const ScalarFlux*>(ele.fluxOut(iLayer()));

    // Retrieve the two different incoming wavevectors in the layer
    const C3& ki = wavevectors.getKi();
    const complex_t kiz = inFlux->getScalarKz();
    const C3 k_i_T{ki.x(), ki.y(), -kiz};
    const C3 k_i_R{ki.x(), ki.y(), +kiz};

    // Retrieve the two different outgoing wavevector bins in the layer
    const C3& kf = wavevectors.getKf();
    const complex_t kfz = outFlux->getScalarKz();
    const C3 k_f_T{kf.x(), kf.y(), +kfz};
    const C3 k_f_R{kf.x(), kf.y(), -kfz};

    // Construct the four different scattering contributions wavevector infos
    const double wavelength = wavevectors.vacuumLambda();
    const WavevectorInfo q_TT(k_i_T, k_f_T, wavelength);
    const WavevectorInfo q_RT(k_i_R, k_f_T, wavelength);
    const WavevectorInfo q_TR(k_i_T, k_f_R, wavelength);
    const WavevectorInfo q_RR(k_i_R, k_f_R, wavelength);

    // Get the four R,T coefficients
    const complex_t T_in = inFlux->getScalarT();
    const complex_t R_in = inFlux->getScalarR();
    const complex_t T_out = outFlux->getScalarT();
    const complex_t R_out = outFlux->getScalarR();

    // The four different scattering contributions;
    // S stands for scattering off the particle,
    // R for reflection off the layer interface.

    // Note that the order of multiplication matters:
    // If a prefactor is 0, then theFF() won't be called.

    std::cout << "DEBUG TERM q=" << q_TT.getQ() << std::endl;
    complex_t ff = m_ff->theFF(q_TT);
    std::cout << "DEBUG TERM ->" << ff << std::endl << std::endl;
    const complex_t term_S = T_in * T_out * ff;
    const complex_t term_RS = R_in * T_out * m_ff->theFF(q_RT);
    const complex_t term_SR = T_in * R_out * m_ff->theFF(q_TR);
    const complex_t term_RSR = R_in * R_out * m_ff->theFF(q_RR);
    std::cout << "DEBUG TERM done with remaining terms" << ff << std::endl << std::endl;

    return term_S + term_RS + term_SR + term_RSR;
}

Eigen::Matrix2cd SumDWBA::coherentPolFF(const DiffuseElement& ele) const
{
    const WavevectorInfo& wavevectors = ele.wavevectorInfo();

    if (!m_i_layer.has_value()) {
        // no slicing, pure Born approximation
        Eigen::Matrix2cd ff_BA = m_ff->thePolFF(wavevectors);
        Eigen::Matrix2cd result;
        result(0, 0) = -ff_BA(1, 0);
        result(0, 1) = +ff_BA(0, 0);
        result(1, 0) = -ff_BA(1, 1);
        result(1, 1) = +ff_BA(0, 1);
        return result;
    }

    const auto* inFlux = dynamic_cast<const MatrixFlux*>(ele.fluxIn(iLayer()));
    const auto* outFlux = dynamic_cast<const MatrixFlux*>(ele.fluxOut(iLayer()));

    // the required wavevectors inside the layer for
    // different eigenmodes and in- and outgoing wavevector;
    const complex_t kix = wavevectors.getKi().x();
    const complex_t kiy = wavevectors.getKi().y();
    const Eigen::Vector2cd& kiz = inFlux->getKz();
    const C3 ki_1R{kix, kiy, +kiz(0)};
    const C3 ki_1T{kix, kiy, -kiz(0)};
    const C3 ki_2R{kix, kiy, +kiz(1)};
    const C3 ki_2T{kix, kiy, -kiz(1)};

    const complex_t kfx = wavevectors.getKf().x();
    const complex_t kfy = wavevectors.getKf().y();
    const Eigen::Vector2cd& kfz = outFlux->getKz();
    const C3 kf_1R{kfx, kfy, -kfz(0)};
    const C3 kf_1T{kfx, kfy, +kfz(0)};
    const C3 kf_2R{kfx, kfy, -kfz(1)};
    const C3 kf_2T{kfx, kfy, +kfz(1)};

    // now each of the 16 matrix terms of the polarized DWBA is calculated:
    // NOTE: when the underlying reflection/transmission coefficients are
    // scalar, the eigenmodes have identical eigenvalues and spin polarization
    // is used as a basis; in this case however the matrices get mixed:
    //     real m_M11 = calculated m_M12
    //     real m_M12 = calculated m_M11
    //     real m_M21 = calculated m_M22
    //     real m_M22 = calculated m_M21
    // since both eigenvalues are identical, this does not influence the result.
    Eigen::Matrix2cd ff_BA;

    double wavelength = wavevectors.vacuumLambda();

    // The following matrices each contain the four polarization conditions
    // (p->p, p->m, m->p, m->m)
    // The first two indices indicate a scattering from the 1/2 eigenstate into
    // the 1/2 eigenstate, while the capital indices indicate a reflection
    // before and/or after the scattering event (first index is in-state,
    // second is out-state; this also applies to the internal matrix indices)
    Eigen::Matrix2cd M11_S, M11_RS, M11_SR, M11_RSR, M12_S, M12_RS, M12_SR, M12_RSR, M21_S, M21_RS,
        M21_SR, M21_RSR, M22_S, M22_RS, M22_SR, M22_RSR;

    // eigenmode 1 -> eigenmode 1: direct scattering
    ff_BA = m_ff->thePolFF({ki_1T, kf_1T, wavelength});
    M11_S(0, 0) = -VecMatVecProduct(outFlux->T1min(), ff_BA, inFlux->T1plus());
    M11_S(0, 1) = +VecMatVecProduct(outFlux->T1plus(), ff_BA, inFlux->T1plus());
    M11_S(1, 0) = -VecMatVecProduct(outFlux->T1min(), ff_BA, inFlux->T1min());
    M11_S(1, 1) = +VecMatVecProduct(outFlux->T1plus(), ff_BA, inFlux->T1min());
    // eigenmode 1 -> eigenmode 1: reflection and then scattering
    ff_BA = m_ff->thePolFF({ki_1R, kf_1T, wavelength});
    M11_RS(0, 0) = -VecMatVecProduct(outFlux->T1min(), ff_BA, inFlux->R1plus());
    M11_RS(0, 1) = +VecMatVecProduct(outFlux->T1plus(), ff_BA, inFlux->R1plus());
    M11_RS(1, 0) = -VecMatVecProduct(outFlux->T1min(), ff_BA, inFlux->R1min());
    M11_RS(1, 1) = +VecMatVecProduct(outFlux->T1plus(), ff_BA, inFlux->R1min());
    // eigenmode 1 -> eigenmode 1: scattering and then reflection
    ff_BA = m_ff->thePolFF({ki_1T, kf_1R, wavelength});
    M11_SR(0, 0) = -VecMatVecProduct(outFlux->R1min(), ff_BA, inFlux->T1plus());
    M11_SR(0, 1) = +VecMatVecProduct(outFlux->R1plus(), ff_BA, inFlux->T1plus());
    M11_SR(1, 0) = -VecMatVecProduct(outFlux->R1min(), ff_BA, inFlux->T1min());
    M11_SR(1, 1) = +VecMatVecProduct(outFlux->R1plus(), ff_BA, inFlux->T1min());
    // eigenmode 1 -> eigenmode 1: reflection, scattering and again reflection
    ff_BA = m_ff->thePolFF({ki_1R, kf_1R, wavelength});
    M11_RSR(0, 0) = -VecMatVecProduct(outFlux->R1min(), ff_BA, inFlux->R1plus());
    M11_RSR(0, 1) = +VecMatVecProduct(outFlux->R1plus(), ff_BA, inFlux->R1plus());
    M11_RSR(1, 0) = -VecMatVecProduct(outFlux->R1min(), ff_BA, inFlux->R1min());
    M11_RSR(1, 1) = +VecMatVecProduct(outFlux->R1plus(), ff_BA, inFlux->R1min());

    // eigenmode 1 -> eigenmode 2: direct scattering
    ff_BA = m_ff->thePolFF({ki_1T, kf_2T, wavelength});
    M12_S(0, 0) = -VecMatVecProduct(outFlux->T2min(), ff_BA, inFlux->T1plus());
    M12_S(0, 1) = +VecMatVecProduct(outFlux->T2plus(), ff_BA, inFlux->T1plus());
    M12_S(1, 0) = -VecMatVecProduct(outFlux->T2min(), ff_BA, inFlux->T1min());
    M12_S(1, 1) = +VecMatVecProduct(outFlux->T2plus(), ff_BA, inFlux->T1min());
    // eigenmode 1 -> eigenmode 2: reflection and then scattering
    ff_BA = m_ff->thePolFF({ki_1R, kf_2T, wavelength});
    M12_RS(0, 0) = -VecMatVecProduct(outFlux->T2min(), ff_BA, inFlux->R1plus());
    M12_RS(0, 1) = +VecMatVecProduct(outFlux->T2plus(), ff_BA, inFlux->R1plus());
    M12_RS(1, 0) = -VecMatVecProduct(outFlux->T2min(), ff_BA, inFlux->R1min());
    M12_RS(1, 1) = +VecMatVecProduct(outFlux->T2plus(), ff_BA, inFlux->R1min());
    // eigenmode 1 -> eigenmode 2: scattering and then reflection
    ff_BA = m_ff->thePolFF({ki_1T, kf_2R, wavelength});
    M12_SR(0, 0) = -VecMatVecProduct(outFlux->R2min(), ff_BA, inFlux->T1plus());
    M12_SR(0, 1) = +VecMatVecProduct(outFlux->R2plus(), ff_BA, inFlux->T1plus());
    M12_SR(1, 0) = -VecMatVecProduct(outFlux->R2min(), ff_BA, inFlux->T1min());
    M12_SR(1, 1) = +VecMatVecProduct(outFlux->R2plus(), ff_BA, inFlux->T1min());
    // eigenmode 1 -> eigenmode 2: reflection, scattering and again reflection
    ff_BA = m_ff->thePolFF({ki_1R, kf_2R, wavelength});
    M12_RSR(0, 0) = -VecMatVecProduct(outFlux->R2min(), ff_BA, inFlux->R1plus());
    M12_RSR(0, 1) = +VecMatVecProduct(outFlux->R2plus(), ff_BA, inFlux->R1plus());
    M12_RSR(1, 0) = -VecMatVecProduct(outFlux->R2min(), ff_BA, inFlux->R1min());
    M12_RSR(1, 1) = +VecMatVecProduct(outFlux->R2plus(), ff_BA, inFlux->R1min());

    // eigenmode 2 -> eigenmode 1: direct scattering
    ff_BA = m_ff->thePolFF({ki_2T, kf_1T, wavelength});
    M21_S(0, 0) = -VecMatVecProduct(outFlux->T1min(), ff_BA, inFlux->T2plus());
    M21_S(0, 1) = +VecMatVecProduct(outFlux->T1plus(), ff_BA, inFlux->T2plus());
    M21_S(1, 0) = -VecMatVecProduct(outFlux->T1min(), ff_BA, inFlux->T2min());
    M21_S(1, 1) = +VecMatVecProduct(outFlux->T1plus(), ff_BA, inFlux->T2min());
    // eigenmode 2 -> eigenmode 1: reflection and then scattering
    ff_BA = m_ff->thePolFF({ki_2R, kf_1T, wavelength});
    M21_RS(0, 0) = -VecMatVecProduct(outFlux->T1min(), ff_BA, inFlux->R2plus());
    M21_RS(0, 1) = +VecMatVecProduct(outFlux->T1plus(), ff_BA, inFlux->R2plus());
    M21_RS(1, 0) = -VecMatVecProduct(outFlux->T1min(), ff_BA, inFlux->R2min());
    M21_RS(1, 1) = +VecMatVecProduct(outFlux->T1plus(), ff_BA, inFlux->R2min());
    // eigenmode 2 -> eigenmode 1: scattering and then reflection
    ff_BA = m_ff->thePolFF({ki_2T, kf_1R, wavelength});
    M21_SR(0, 0) = -VecMatVecProduct(outFlux->R1min(), ff_BA, inFlux->T2plus());
    M21_SR(0, 1) = +VecMatVecProduct(outFlux->R1plus(), ff_BA, inFlux->T2plus());
    M21_SR(1, 0) = -VecMatVecProduct(outFlux->R1min(), ff_BA, inFlux->T2min());
    M21_SR(1, 1) = +VecMatVecProduct(outFlux->R1plus(), ff_BA, inFlux->T2min());
    // eigenmode 2 -> eigenmode 1: reflection, scattering and again reflection
    ff_BA = m_ff->thePolFF({ki_2R, kf_1R, wavelength});
    M21_RSR(0, 0) = -VecMatVecProduct(outFlux->R1min(), ff_BA, inFlux->R2plus());
    M21_RSR(0, 1) = +VecMatVecProduct(outFlux->R1plus(), ff_BA, inFlux->R2plus());
    M21_RSR(1, 0) = -VecMatVecProduct(outFlux->R1min(), ff_BA, inFlux->R2min());
    M21_RSR(1, 1) = +VecMatVecProduct(outFlux->R1plus(), ff_BA, inFlux->R2min());

    // eigenmode 2 -> eigenmode 2: direct scattering
    ff_BA = m_ff->thePolFF({ki_2T, kf_2T, wavelength});
    M22_S(0, 0) = -VecMatVecProduct(outFlux->T2min(), ff_BA, inFlux->T2plus());
    M22_S(0, 1) = +VecMatVecProduct(outFlux->T2plus(), ff_BA, inFlux->T2plus());
    M22_S(1, 0) = -VecMatVecProduct(outFlux->T2min(), ff_BA, inFlux->T2min());
    M22_S(1, 1) = +VecMatVecProduct(outFlux->T2plus(), ff_BA, inFlux->T2min());
    // eigenmode 2 -> eigenmode 2: reflection and then scattering
    ff_BA = m_ff->thePolFF({ki_2R, kf_2T, wavelength});
    M22_RS(0, 0) = -VecMatVecProduct(outFlux->T2min(), ff_BA, inFlux->R2plus());
    M22_RS(0, 1) = +VecMatVecProduct(outFlux->T2plus(), ff_BA, inFlux->R2plus());
    M22_RS(1, 0) = -VecMatVecProduct(outFlux->T2min(), ff_BA, inFlux->R2min());
    M22_RS(1, 1) = +VecMatVecProduct(outFlux->T2plus(), ff_BA, inFlux->R2min());
    // eigenmode 2 -> eigenmode 2: scattering and then reflection
    ff_BA = m_ff->thePolFF({ki_2T, kf_2R, wavelength});
    M22_SR(0, 0) = -VecMatVecProduct(outFlux->R2min(), ff_BA, inFlux->T2plus());
    M22_SR(0, 1) = +VecMatVecProduct(outFlux->R2plus(), ff_BA, inFlux->T2plus());
    M22_SR(1, 0) = -VecMatVecProduct(outFlux->R2min(), ff_BA, inFlux->T2min());
    M22_SR(1, 1) = +VecMatVecProduct(outFlux->R2plus(), ff_BA, inFlux->T2min());
    // eigenmode 2 -> eigenmode 2: reflection, scattering and again reflection
    ff_BA = m_ff->thePolFF({ki_2R, kf_2R, wavelength});
    M22_RSR(0, 0) = -VecMatVecProduct(outFlux->R2min(), ff_BA, inFlux->R2plus());
    M22_RSR(0, 1) = +VecMatVecProduct(outFlux->R2plus(), ff_BA, inFlux->R2plus());
    M22_RSR(1, 0) = -VecMatVecProduct(outFlux->R2min(), ff_BA, inFlux->R2min());
    M22_RSR(1, 1) = +VecMatVecProduct(outFlux->R2plus(), ff_BA, inFlux->R2min());

    return M11_S + M11_RS + M11_SR + M11_RSR + M12_S + M12_RS + M12_SR + M12_RSR + M21_S + M21_RS
           + M21_SR + M21_RSR + M22_S + M22_RS + M22_SR + M22_RSR;
}