ReSample.cpp

Go to the documentation of this file.

//  ************************************************************************************************
//
//  BornAgain: simulate and fit reflection and scattering
//
//! @file      Resample/Processed/ReSample.cpp
//! @brief     Implements class reSample.
//!
//! @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/Processed/ReSample.h"
#include "Base/Util/Assert.h"
#include "Resample/Coherence/CoherentFFSum.h"
#include "Resample/Coherence/SumDWBA.h"
#include "Resample/Flux/IFlux.h"
#include "Resample/Options/SimulationOptions.h"
#include "Resample/Processed/ParticleInSlice.h"
#include "Resample/Processed/ParticleRegions.h"
#include "Resample/Processed/ReLayout.h"
#include "Resample/Processed/Slicer.h"
#include "Resample/Specular/ComputeFluxMagnetic.h"
#include "Resample/Specular/ComputeFluxScalar.h"
#include "Sample/Aggregate/IInterference.h"
#include "Sample/Aggregate/ParticleLayout.h"
#include "Sample/Interface/LayerRoughness.h"
#include "Sample/Material/MaterialBySLDImpl.h"
#include "Sample/Material/MaterialUtils.h"
#include "Sample/Material/RefractiveMaterialImpl.h"
#include "Sample/Multilayer/Layer.h"
#include "Sample/Multilayer/MultiLayer.h"
#include "Sample/Multilayer/MultilayerUtils.h"
#include "Sample/Particle/IParticle.h"
 
namespace {
 
//! Returns the average of an intensive property.
//! Used for SLD (T=complex) and magnetization (T=vector).
 
template <class T>
T averageData(const Material& layer_mat, const Admixtures& admixtures,
              std::function<T(const Material&)> average)
{
    const T layer_data = average(layer_mat);
    T averaged_data = layer_data;
    for (const OneAdmixture& admix : admixtures)
        averaged_data += admix.fraction * (average(admix.material) - layer_data);
    return averaged_data;
}
 
// Tested by Tests/Unit/Sim/Sample/MaterialTest.cpp
 
Material createAveragedMaterial(const Material& layer_mat, const Admixtures& admixtures)
{
    const std::string avge_name = layer_mat.materialName() + "_avg";
 
    const R3 avge_magnetization = averageData<R3>(
        layer_mat, admixtures, [](const Material& mat) { return mat.magnetization(); });
 
    // Is there a uniform material type (refractive or SLD)?
    std::vector<const Material*> materials;
    materials.push_back(&layer_mat);
    for (const OneAdmixture& admix : admixtures)
        materials.push_back(&admix.material);
 
    const MATERIAL_TYPES common_type = MaterialUtils::checkMaterialTypes(materials);
    if (common_type == MATERIAL_TYPES::InvalidMaterialType)
        throw std::runtime_error("Error in createAveragedMaterial(): non-default materials of "
                                 "different types used simultaneously.");
 
    if (common_type == MATERIAL_TYPES::RefractiveMaterial) {
        // avrData returns (1 - mdc)^2 - 1, where mdc is material data conjugate
        auto avrData = [](const Material& mat) -> complex_t {
            const complex_t mdc = std::conj(mat.materialData());
            return mdc * mdc - 2.0 * mdc;
        };
        const complex_t avr_mat_data = std::conj(
            1.0 - std::sqrt(1.0 + averageData<complex_t>(layer_mat, admixtures, avrData)));
        return RefractiveMaterial(avge_name, avr_mat_data.real(), avr_mat_data.imag(),
                                  avge_magnetization);
    }
    if (common_type == MATERIAL_TYPES::MaterialBySLD) {
        complex_t (*avrData)(const Material&) = [](const Material& mat) {
            return mat.materialData();
        };
        const auto avr_mat_data = averageData<complex_t>(layer_mat, admixtures, avrData);
        return MaterialBySLD(avge_name, avr_mat_data.real(), avr_mat_data.imag(),
                             avge_magnetization);
    }
    ASSERT(0);
}
 
void checkVolumeFractions(const Admixtures& admixtures)
{
    double total_fraction = 0.0;
    for (const auto& admix : admixtures)
        total_fraction += admix.fraction;
    ASSERT(total_fraction >= 0);
    if (total_fraction > 1)
        throw std::runtime_error("average slices: "
                                 "total volumetric fraction of particles exceeds 1!");
}
 
SliceStack refineStack(const SliceStack& stack, const std::vector<reLayout>& layouts)
{
    SliceStack result = stack;
    for (const reLayout& layout : layouts) {
        for (const auto& entry : layout.regionMap()) {
            const size_t i_slice = entry.first;
            if (i_slice == 0 || i_slice == result.size() - 1)
                continue; // skip semi-infinite layers
            const auto slice_mat = result[i_slice].material();
            checkVolumeFractions(entry.second);
            result[i_slice].setMaterial(createAveragedMaterial(slice_mat, entry.second));
        }
    }
    ASSERT(result.size() == stack.size());
    return result;
}
 
//! Returns a SliceStack that refines the layer structure of sample,
//! taking into account the average SLD of layer material and particle decoration.
//!
//! For use in reSample constructor.
 
SliceStack slicify(const MultiLayer& sample, const SimulationOptions& options)
{
    SliceStack result(sample.roughnessModel());
    if (sample.numberOfLayers() == 0)
        return result;
    const bool use_slicing = options.useAvgMaterials() && sample.numberOfLayers() > 1;
    const auto layer_limits = Compute::Slicing::particleRegions(sample, use_slicing);
    for (size_t i = 0; i < sample.numberOfLayers(); ++i) {
        const Layer* const layer = sample.layer(i);
        const size_t N = layer->numberOfSlices();
        if (N == 0)
            throw std::runtime_error("Cannot process layer with numberOfSlices=0");
        const ZLimits& slice_limits = layer_limits[i];
        double tl = layer->thickness();
        const Material* const material = layer->material();
        const LayerRoughness* roughness = SampleUtils::Multilayer::LayerTopRoughness(sample, i);
        if (roughness && roughness->sigma() <= 0)
            roughness = nullptr;
        // if no slicing is needed, create single slice for the layer
        if (!slice_limits.isFinite() || N == 0) {
            if (i == sample.numberOfLayers() - 1)
                tl = 0.0;
            if (i == 0)
                result.addTopSlice(tl, *material);
            else
                result.addSlice(tl, *material, roughness);
            continue;
        }
        const double top = slice_limits.zTop();
        const double bottom = slice_limits.zBottom();
        // top layer
        if (i == 0) {
            if (top <= 0)
                throw std::runtime_error("reSample::reSample: "
                                         "top limit for top layer must be > 0.");
            result.addTopSlice(top, *material); // semiinfinite top slice
            result.addNSlices(N, top - bottom, *material);
            if (bottom > 0)
                result.addSlice(bottom, *material);
        }
        // middle or bottom layer
        else {
            if (top < 0) {
                result.addSlice(-top, *material, roughness);
                result.addNSlices(N, top - bottom, *material);
            } else {
                result.addNSlices(N, top - bottom, *material, roughness);
            }
            // middle layer
            if (i < sample.numberOfLayers() - 1 && bottom > -tl)
                result.addSlice(bottom + tl, *material);
            // bottom layer
            if (i == sample.numberOfLayers() - 1)
                result.addSlice(0.0, *material); // semiinfinite bottom slice
        }
    }
    return result;
}
 
size_t zToSliceIndex(double z, const SliceStack& slices)
{
    auto n_layers = slices.size();
    if (n_layers < 2 || z > 0.0)
        return 0;
 
    double z0 = slices[0].zBottom();
    size_t result = 0;
    while (result < n_layers - 1) {
        ++result;
        if (slices[result].zBottom() - z0 < z)
            break;
    }
    return result;
}
 
std::pair<size_t, size_t> SliceIndexSpan(const IParticle& particle, const SliceStack& slices,
                                         double z_ref)
{
    const ZLimits zSpan = Compute::Slicing::zSpan(&particle);
    const double zbottom = zSpan.zBottom();
    const double ztop = zSpan.zTop();
    const double eps =
        (ztop - zbottom) * 1e-6; // allow for relatively small crossing
                                 // due to numerical approximations (like rotation over 180 degrees)
    const double zmax = ztop + z_ref - eps;
    const double zmin = zbottom + z_ref + eps;
    size_t top_index = zToSliceIndex(zmax, slices);
    const size_t bottom_index = zToSliceIndex(zmin, slices);
    if (top_index > bottom_index) // happens for zero size particles
        top_index = bottom_index;
    return {top_index, bottom_index};
}
 
reLayout makereLayout(const ParticleLayout& layout, const SliceStack& slices, double z_ref,
                      bool polarized)
{
    const double layout_abundance = layout.totalAbundance();
    ASSERT(layout_abundance > 0);
    const double surface_density = layout.weightedParticleSurfaceDensity();
 
    std::vector<std::unique_ptr<const CoherentFFSum>> formfactors;
    std::map<size_t, Admixtures> slice2admixtures;
 
    for (const auto& particle : layout.particles()) {
 
        std::vector<std::shared_ptr<const SumDWBA>> terms;
 
        for (const auto& subparticle : particle->decompose()) {
            const auto slice_indices = SliceIndexSpan(*subparticle, slices, z_ref);
            bool single_layer = (slice_indices.first == slice_indices.second);
            double z0 = slices.at(0).zBottom();
            for (size_t iSlice = slice_indices.first; iSlice < slice_indices.second + 1; ++iSlice) {
                const Slice& slice = slices[iSlice];
                double z_top_i = iSlice == 0 ? 0 : slice.zTop() - z0;
                R3 translation(0.0, 0.0, z_ref - z_top_i);
                z_ref = z_top_i; // subparticle now has coordinates relative to z_top_i
                subparticle->translate(translation);
                // if subparticle is contained in this layer, set limits to infinite:
                ZLimits limits;
                if (!single_layer) {
                    double z1 = slice.zTopOr0();
                    limits = {slice.zBottom() - z1, slice.zTop() - z1};
                }
                ParticleInSlice pis =
                    Compute::Slicing::createParticleInSlice(subparticle.get(), limits);
                pis.particleSlice->setAmbientMaterial(slices.at(iSlice).material());
 
                std::unique_ptr<SumDWBA> computer;
                if (slices.size() > 1)
                    computer = std::make_unique<SumDWBA>(*pis.particleSlice, iSlice);
                else
                    computer = std::make_unique<SumDWBA>(*pis.particleSlice);
                terms.emplace_back(computer.release());
 
                double factor = particle->abundance() / layout_abundance * surface_density;
                if (double thickness = slice.thicknessOr0(); thickness > 0.0)
                    factor /= thickness;
                for (auto& admixture : pis.admixtures)
                    admixture.fraction *= factor;
 
                slice2admixtures[iSlice].insert(slice2admixtures[iSlice].begin(),
                                                pis.admixtures.begin(), pis.admixtures.end());
            }
        }
 
        formfactors.emplace_back(
            std::make_unique<CoherentFFSum>(particle->abundance() / layout_abundance, terms));
    }
 
    const auto* iff = layout.interferenceFunction();
 
    return reLayout(polarized, surface_density, std::move(formfactors),
                    iff ? iff->clone() : nullptr, std::move(slice2admixtures));
}
 
//! Returns collection of all reLayout%s, for use in reSample constructor.
 
std::vector<reLayout> collectLayouts(const MultiLayer& sample, const SliceStack& slices,
                                     bool polarized)
{
    std::vector<reLayout> result;
 
    double z_ref = -slices.front().zBottom();
 
    for (size_t i = 0; i < sample.numberOfLayers(); ++i) {
        if (i > 1)
            z_ref -= sample.layer(i - 1)->thickness();
        for (const auto* layout : sample.layer(i)->layouts())
            result.emplace_back(makereLayout(*layout, slices, z_ref, polarized));
    }
    return result;
}
 
} // namespace
 
 
//  ************************************************************************************************
//  class implementation
//  ************************************************************************************************
 
reSample reSample::make(const MultiLayer& sample, const SimulationOptions& options,
                        bool forcePolarized)
{
    const bool polarized = forcePolarized || sample.isMagnetic() || sample.externalField() != R3();
 
    // slices1: accounting for roughness, but not yet for particles
    const SliceStack slices1 = slicify(sample, options).setBField(sample.externalField());
 
    std::vector<reLayout> layouts = collectLayouts(sample, slices1, polarized);
 
    const SliceStack refined_stack =
        options.useAvgMaterials() ? refineStack(slices1, layouts) : slices1;
 
    return reSample(sample, polarized, std::move(layouts), refined_stack);
}
 
reSample::reSample(const MultiLayer& sample, bool polarized, std::vector<reLayout>&& layouts,
                   const SliceStack& refined_stack)
    : m_sample(sample)
    , m_polarized(polarized)
    , m_layouts(std::move(layouts))
    , m_stack(refined_stack)
{
}
 
reSample::~reSample() = default;
 
size_t reSample::numberOfSlices() const
{
    return m_stack.size();
}
 
const SliceStack& reSample::averageSlices() const
{
    return m_stack;
}
 
const Slice& reSample::avgeSlice(size_t i) const
{
    return m_stack.at(i);
}
 
const std::vector<reLayout>& reSample::layouts() const
{
    return m_layouts;
}
 
double reSample::sliceTopZ(size_t i) const
{
    return m_stack.at(i).zTop();
}
 
double reSample::sliceBottomZ(size_t i) const
{
    return m_stack.at(i).zBottom();
}
 
bool reSample::containsMagneticMaterial() const
{
    return m_polarized;
}
 
bool reSample::polarizing() const
{
    return containsMagneticMaterial() || m_sample.externalField() != R3{};
}
 
bool reSample::hasRoughness() const
{
    for (const auto& slice : m_stack)
        if (slice.topRoughness())
            return true;
    return false;
}
 
Fluxes reSample::fluxesIn(const R3& k) const
{
    if (m_polarized)
        return Compute::SpecularMagnetic::fluxes(m_stack, k, true);
    return Compute::SpecularScalar::fluxes(m_stack, k);
}
 
Fluxes reSample::fluxesOut(const R3& k) const
{
    if (m_polarized)
        return Compute::SpecularMagnetic::fluxes(m_stack, -k, false);
    return Compute::SpecularScalar::fluxes(m_stack, -k);
}
 
double reSample::crossCorrSpectralFun(const R3 kvec, size_t j, size_t k) const
{
    const double xCorrLength = m_sample.crossCorrLength();
    if (xCorrLength <= 0.0)
        return 0.0;
    const double z_j = sliceBottomZ(j);
    const double z_k = sliceBottomZ(k);
    const LayerRoughness* rough_j = m_stack.at(j + 1).topRoughness();
    const LayerRoughness* rough_k = m_stack.at(k + 1).topRoughness();
    if (!rough_j || !rough_k)
        return 0.0;
    const double sigma_j = rough_j->sigma();
    const double sigma_k = rough_k->sigma();
    if (sigma_j <= 0 || sigma_k <= 0)
        return 0.0;
    return 0.5
           * ((sigma_k / sigma_j) * rough_j->spectralFunction(kvec)
              + (sigma_j / sigma_k) * rough_k->spectralFunction(kvec))
           * std::exp(-1 * std::abs(z_j - z_k) / xCorrLength);
}