ComputeFluxScalar.cpp

Go to the documentation of this file.

//  ************************************************************************************************
//
//  BornAgain: simulate and fit reflection and scattering
//
//! @file      Resample/Specular/ComputeFluxScalar.cpp
//! @brief     Implements class SpecularScalarStrategy.
//!
//! @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/Specular/ComputeFluxScalar.h"
#include "Base/Math/Constants.h"
#include "Base/Math/Functions.h"
#include "Base/Util/Assert.h"
#include "Resample/Flux/ScalarFlux.h"
#include "Resample/Slice/KzComputation.h"
#include "Resample/Slice/SliceStack.h"
#include "Sample/Interface/LayerRoughness.h"
#include "Sample/Multilayer/Layer.h"
#include <stdexcept>
 
namespace {
 
std::pair<complex_t, complex_t> transition(complex_t kzi, complex_t kzi1, double sigma,
                                           const RoughnessModel& r_model)
{
    if (r_model == RoughnessModel::NEVOT_CROCE) {
        // Roughness is modelled by a Gaussian profile, i.e. Nevot-Croce factors for the
        // reflection coefficients.
        // Implementation follows A. Gibaud and G. Vignaud, in X-ray and Neutron Reflectivity,
        // edited by J. Daillant and A. Gibaud, volume 770 of Lecture Notes in Physics (2009)
        complex_t roughness_diff = 1;
        complex_t roughness_sum = 1;
        if (sigma > 0.0) {
            roughness_diff = std::exp(-(kzi1 - kzi) * (kzi1 - kzi) * sigma * sigma / 2.);
            roughness_sum = std::exp(-(kzi1 + kzi) * (kzi1 + kzi) * sigma * sigma / 2.);
        }
        const complex_t kz_ratio = kzi1 / kzi;
 
        const complex_t a00 = 0.5 * (1. + kz_ratio) * roughness_diff;
        const complex_t a01 = 0.5 * (1. - kz_ratio) * roughness_sum;
 
        return {a00, a01};
    }
    // TANH model assumed
    // Roughness is modelled by tanh profile
    // [e.g. Bahr, Press, et al, Phys. Rev. B, vol. 47 (8), p. 4385 (1993)].
    complex_t roughness = 1;
    if (sigma > 0.0) {
        const double sigeff = std::pow(M_PI_2, 1.5) * sigma;
        roughness = std::sqrt(Math::tanhc(sigeff * kzi1) / Math::tanhc(sigeff * kzi));
    }
    const complex_t inv_roughness = 1.0 / roughness;
    const complex_t kz_ratio = kzi1 / kzi * roughness;
 
    const complex_t a00 = 0.5 * (inv_roughness + kz_ratio);
    const complex_t a01 = 0.5 * (inv_roughness - kz_ratio);
 
    return {a00, a01};
}
 
std::vector<Spinor> computeTR(const SliceStack& slices, const std::vector<complex_t>& kz,
                              const RoughnessModel& r_model, bool top_exit)
{
    const size_t N = slices.size();
    std::vector<Spinor> TR(N, {1., 0.});
 
    if (N == 1) // If only one layer present, there's nothing left to calculate
        return TR;
 
    // Index reversal for bottom exit
    std::vector<size_t> X(N, 0);
    for (size_t i = 0; i < N; ++i)
        X[i] = top_exit ? i : N - 1 - i;
 
    if (kz[X[0]] == 0.0) { // If kz in layer 0 is zero, R0 = -T0 and all others equal to 0
        TR[X[0]] = {1.0, -1.0};
        for (size_t i = 1; i < N; ++i)
            TR[X[i]] = Spinor(0, 0);
        return TR;
    }
 
    // Calculate transmission/refraction coefficients t_r for each layer, from bottom to top.
    TR[X[N - 1]] = {1.0, 0.0};
    std::vector<complex_t> factors(N - 1);
    for (size_t i = N - 1; i > 0; i--) {
        size_t jthis = X[i - 1];
        size_t jlast = X[i];
        const auto* roughness = slices.bottomRoughness(jthis); // TODO verify
        const double sigma = roughness ? roughness->sigma() : 0.;
 
        const auto [mp, mm] = transition(kz[jthis], kz[jlast], sigma, r_model);
 
        const complex_t delta = exp_I(kz[jthis] * slices[jthis].thicknessOr0());
 
        complex_t S = delta / (mp + mm * TR[jlast].v);
        factors[i - 1] = S;
        TR[jthis].v = delta * (mm + mp * TR[jlast].v) * S;
    }
 
    // Now correct all amplitudes by dividing by the remaining factors in forward direction.
    // At some point this divison underflows; from there on all further amplitudes are set to zero.
    complex_t dampingFactor = 1;
    for (size_t i = 1; i < N; ++i) {
        dampingFactor = dampingFactor * factors[i - 1];
        TR[X[i]] = Spinor(dampingFactor, TR[X[i]].v * dampingFactor);
    }
 
    return TR;
}
 
} // namespace
 
//  ************************************************************************************************
//  namespace implementation
//  ************************************************************************************************
 
Fluxes Compute::SpecularScalar::fluxes(const SliceStack& slices, const R3& k)
{
    const bool top_exit = k.z() <= 0; // true if source or detector pixel are at z>=0
    std::vector<complex_t> kz = Compute::Kz::computeReducedKz(slices, k);
    ASSERT(slices.size() == kz.size());
 
    const std::vector<Spinor> TR = computeTR(slices, kz, slices.roughnessModel(), top_exit);
 
    Fluxes result;
    for (size_t i = 0; i < kz.size(); ++i)
        result.emplace_back(std::make_unique<const ScalarFlux>(kz[i], TR[i]));
    return result;
}
 
complex_t Compute::SpecularScalar::topLayerR(const SliceStack& slices,
                                             const std::vector<complex_t>& kz)
{
    const RoughnessModel& r_model = slices.roughnessModel();
 
    ASSERT(slices.size() == kz.size());
    const size_t N = slices.size();
    if (N == 1)
        return 0.; // only one layer present, there's nothing left to calculate
    if (kz[0] == 0.)
        return -1.;
 
    complex_t R_i1 = 0.;
 
    for (size_t i = N - 1; i > 0; i--) {
        double sigma = 0.0;
        if (const auto* const roughness = slices.bottomRoughness(i - 1))
            sigma = roughness->sigma();
 
        const auto [mp, mm] = transition(kz[i - 1], kz[i], sigma, r_model);
 
        const complex_t delta = exp_I(kz[i - 1] * slices[i - 1].thicknessOr0());
 
        R_i1 = pow(delta, 2) * (mm + mp * R_i1) / (mp + mm * R_i1);
    }
 
    return R_i1;
}