RoughMultiLayerContribution.cpp

Go to the documentation of this file.

//  ************************************************************************************************
//
//  BornAgain: simulate and fit reflection and scattering
//
//! @file      Sim/Contrib/RoughMultiLayerContribution.cpp
//! @brief     Implements class RoughMultiLayerContribution.
//!
//! @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 "Sim/Contrib/RoughMultiLayerContribution.h"
#include "Base/Math/Constants.h"
#include "Base/Util/Assert.h"
#include "Resample/Element/DiffuseElement.h"
#include "Resample/Flux/ScalarFlux.h"
#include "Resample/Processed/ReSample.h"
#include "Sample/Interface/LayerInterface.h"
#include "Sample/Interface/LayerRoughness.h"
#include "Sample/Multilayer/Layer.h"
#include "Sample/Multilayer/MultiLayer.h"
 
#include <cerfcpp.h>
 
// Diffuse scattering from rough interfaces is modelled after
// Phys. Rev. B, vol. 51 (4), p. 2311 (1995)
 
namespace {
 
complex_t h_plus(complex_t z)
{
    return 0.5 * cerfcx(-mul_I(z) / std::sqrt(2.0));
}
complex_t h_min(complex_t z)
{
    return 0.5 * cerfcx(mul_I(z) / std::sqrt(2.0));
}
 
} // namespace
 
 
RoughMultiLayerContribution::RoughMultiLayerContribution(const reSample& re_sample)
    : m_re_sample{re_sample}
{
}
 
void RoughMultiLayerContribution::compute(DiffuseElement& ele) const
{
    if (ele.alphaMean() < 0.0)
        return;
    const size_t n_slices = m_re_sample.numberOfSlices();
    R3 q = ele.meanQ();
    double wavelength = ele.wavelength();
    double autocorr(0.0);
    complex_t crosscorr(0.0, 0.0);
 
    std::vector<complex_t> rterm(n_slices - 1);
    std::vector<complex_t> sterm(n_slices - 1);
 
    for (size_t i = 0; i + 1 < n_slices; i++) {
        rterm[i] = get_refractive_term(i, wavelength);
        sterm[i] = get_sum8terms(i, ele);
    }
    for (size_t i = 0; i + 1 < n_slices; i++) {
        const LayerRoughness* rough = m_re_sample.avgeSlice(i + 1).topRoughness();
        if (rough)
            autocorr += std::norm(rterm[i]) * std::norm(sterm[i]) * rough->spectralFunction(q);
    }
    // cross correlation between layers
    if (m_re_sample.sample().crossCorrLength() != 0.0) {
        for (size_t j = 0; j < n_slices - 1; j++) {
            for (size_t k = 0; k < n_slices - 1; k++) {
                if (j == k)
                    continue;
                crosscorr += rterm[j] * sterm[j] * m_re_sample.crossCorrSpectralFun(q, j, k)
                             * std::conj(rterm[k]) * std::conj(sterm[k]);
            }
        }
    }
    //! @TODO clarify complex vs double
    ele.addIntensity((autocorr + crosscorr.real()) * M_PI / 4. / wavelength / wavelength);
}
 
complex_t RoughMultiLayerContribution::get_refractive_term(size_t i_layer, double wavelength) const
{
    return m_re_sample.avgeSlice(i_layer).material().refractiveIndex2(wavelength)
           - m_re_sample.avgeSlice(i_layer + 1).material().refractiveIndex2(wavelength);
}
 
complex_t RoughMultiLayerContribution::get_sum8terms(size_t i_layer,
                                                     const DiffuseElement& ele) const
{
    const auto* in_plus = dynamic_cast<const ScalarFlux*>(ele.fluxIn(i_layer));
    const auto* out_plus = dynamic_cast<const ScalarFlux*>(ele.fluxOut(i_layer));
    const auto* in_minus = dynamic_cast<const ScalarFlux*>(ele.fluxIn(i_layer + 1));
    const auto* out_minus = dynamic_cast<const ScalarFlux*>(ele.fluxOut(i_layer + 1));
    ASSERT(in_plus && out_plus && in_minus && out_minus);
 
    const complex_t kiz_plus = in_plus->getScalarKz();
    const complex_t kfz_plus = out_plus->getScalarKz();
    const complex_t qz1_plus = -kiz_plus - kfz_plus;
    const complex_t qz2_plus = -kiz_plus + kfz_plus;
    const complex_t qz3_plus = -qz2_plus;
    const complex_t qz4_plus = -qz1_plus;
 
    const double thickness = m_re_sample.avgeSlice(i_layer).thicknessOr0();
    const complex_t T_in_plus = in_plus->getScalarT() * exp_I(kiz_plus * thickness);
    const complex_t R_in_plus = in_plus->getScalarR() * exp_I(-kiz_plus * thickness);
    const complex_t T_out_plus = out_plus->getScalarT() * exp_I(kfz_plus * thickness);
    const complex_t R_out_plus = out_plus->getScalarR() * exp_I(-kfz_plus * thickness);
 
    const complex_t kiz_minus = in_minus->getScalarKz();
    const complex_t kfz_minus = out_minus->getScalarKz();
    const complex_t qz1_minus = -kiz_minus - kfz_minus;
    const complex_t qz2_minus = -kiz_minus + kfz_minus;
    const complex_t qz3_minus = -qz2_minus;
    const complex_t qz4_minus = -qz1_minus;
 
    const LayerRoughness* roughness = m_re_sample.averageSlices().bottomRoughness(i_layer);
    const double sigma = roughness ? roughness->sigma() : 0.;
    const complex_t term1 = T_in_plus * T_out_plus * h_plus(qz1_plus * sigma);
    const complex_t term2 = T_in_plus * R_out_plus * h_plus(qz2_plus * sigma);
    const complex_t term3 = R_in_plus * T_out_plus * h_plus(qz3_plus * sigma);
    const complex_t term4 = R_in_plus * R_out_plus * h_plus(qz4_plus * sigma);
    const complex_t term5 =
        in_minus->getScalarT() * out_minus->getScalarT() * h_min(qz1_minus * sigma);
    const complex_t term6 =
        in_minus->getScalarT() * out_minus->getScalarR() * h_min(qz2_minus * sigma);
    const complex_t term7 =
        in_minus->getScalarR() * out_minus->getScalarT() * h_min(qz3_minus * sigma);
    const complex_t term8 =
        in_minus->getScalarR() * out_minus->getScalarR() * h_min(qz4_minus * sigma);
 
    return term1 + term2 + term3 + term4 + term5 + term6 + term7 + term8;
}