CoordSystem1D.cpp

Go to the documentation of this file.

//  ************************************************************************************************
//
//  BornAgain: simulate and fit reflection and scattering
//
//! @file      Device/Coord/CoordSystem1D.cpp
//! @brief     Implements CoordSystem1D class and derived classes.
//!
//! @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 "Device/Coord/CoordSystem1D.h"
#include "Base/Axis/FixedBinAxis.h"
#include "Base/Axis/PointwiseAxis.h"
#include "Base/Const/Units.h"
#include "Base/Math/Constants.h"
#include "Base/Util/Assert.h"
#include "Device/Coord/AxisNames.h"
#include "Device/Data/Datafield.h"
#include <cmath>
 
namespace {
 
double getQ(double wavelength, double angle)
{
    return 4.0 * M_PI * std::sin(angle) / wavelength;
}
 
double getInvQ(double wavelength, double q)
{
    double sin_angle = q * wavelength / (4.0 * M_PI);
    return std::asin(sin_angle);
}
 
double backTransform(double value, Coords coords, double wavelength)
{
    switch (coords) {
    case Coords::RADIANS:
        return value;
    case Coords::DEGREES:
        return Units::deg2rad(value);
    case Coords::QSPACE:
        return getInvQ(wavelength, value);
    default:
        throw std::runtime_error("CoordSystem1D/backTransform: invalid coord system");
    }
}
 
std::string angularAxisName(const Coords units)
{
    const auto& name_map = DataUtils::AxisNames::specAxis;
    const auto& it = name_map.find(units == Coords::UNDEFINED ? Coords::DEGREES : units);
    ASSERT(it != name_map.cend());
    return it->second;
}
 
PointwiseAxis* createAxisFrom(const IAxis& axis, Coords coords, const std::string& name,
                              double wavelength)
{
    std::vector<double> ret;
    ret.reserve(axis.size());
    for (const double value : axis.binCenters())
        ret.emplace_back(backTransform(value, coords, wavelength));
    return new PointwiseAxis(name, ret);
}
 
} // namespace
 
 
//  ************************************************************************************************
//  class CoordSystem1D
//  ************************************************************************************************
 
CoordSystem1D::CoordSystem1D(const IAxis*&& axis)
    : m_axis(axis)
{
}
 
//! Returns the size of underlying axis.
size_t CoordSystem1D::axisSize(size_t i_axis) const
{
    ASSERT(i_axis == 0);
    return coordinateAxis()->size();
}
 
double CoordSystem1D::calculateMin(size_t i_axis, Coords units) const
{
    ASSERT(i_axis == 0);
    units = substituteDefaultUnits(units);
    if (units == Coords::NBINS)
        return 0.0;
    auto translator = getTraslatorTo(units);
    return translator(coordinateAxis()->binCenter(0));
}
 
double CoordSystem1D::calculateMax(size_t i_axis, Coords units) const
{
    ASSERT(i_axis == 0);
    units = substituteDefaultUnits(units);
    const IAxis* const coordinate_axis = coordinateAxis();
    if (units == Coords::NBINS)
        return static_cast<double>(coordinate_axis->size());
    auto translator = getTraslatorTo(units);
    return translator(coordinate_axis->binCenter(coordinate_axis->size() - 1));
}
 
IAxis* CoordSystem1D::createConvertedAxis(size_t i_axis, Coords units) const
{
    ASSERT(i_axis == 0);
    units = substituteDefaultUnits(units);
    if (units == Coords::NBINS)
        return new FixedBinAxis(axisName(0, units), coordinateAxis()->size(),
                                calculateMin(0, units), calculateMax(0, units));
 
    std::function<double(double)> translator = getTraslatorTo(units);
    auto coords = coordinateAxis()->binCenters();
    for (size_t i = 0, size = coords.size(); i < size; ++i)
        coords[i] = translator(coords[i]);
    return new PointwiseAxis(axisName(0, units), coords);
}
 
Datafield* CoordSystem1D::createConvertedData(const Datafield& data, Coords units) const
{
    ASSERT(data.rank() == 1);
 
    auto q_axis = createConvertedAxis(0, units);
    Datafield* result = new Datafield({q_axis});
 
    if (units != Coords::RQ4) {
        result->setVector(data.flatVector());
        return result;
    }
 
    for (size_t i = 0, size = result->size(); i < size; ++i)
        (*result)[i] = data[i] * std::pow((*q_axis)[i], 4);
    return result;
}
 
 
//  ************************************************************************************************
//  class AngularReflectometryCoordinates
//  ************************************************************************************************
 
AngularReflectometryCoordinates::AngularReflectometryCoordinates(double wavelength,
                                                                 const IAxis& axis,
                                                                 Coords axis_units)
    : CoordSystem1D(createAxisFrom(axis, axis_units, angularAxisName(axis_units), wavelength))
    , m_wavelength(wavelength)
{
    if (m_axis->min() < 0 || m_axis->max() > M_PI_2)
        throw std::runtime_error("Error in CoordSystem1D: input axis range is out of bounds");
}
 
AngularReflectometryCoordinates::AngularReflectometryCoordinates(
    const AngularReflectometryCoordinates& other)
    : CoordSystem1D(other.m_axis->clone())
    , m_wavelength(other.m_wavelength)
{
}
 
AngularReflectometryCoordinates::~AngularReflectometryCoordinates() = default;
 
AngularReflectometryCoordinates* AngularReflectometryCoordinates::clone() const
{
    return new AngularReflectometryCoordinates(*this);
}
 
std::vector<Coords> AngularReflectometryCoordinates::availableUnits() const
{
    return {Coords::NBINS, Coords::RADIANS, Coords::DEGREES, Coords::QSPACE, Coords::RQ4};
}
 
std::vector<std::map<Coords, std::string>> AngularReflectometryCoordinates::createNameMaps() const
{
    return {DataUtils::AxisNames::specAxis};
}
 
std::function<double(double)> AngularReflectometryCoordinates::getTraslatorTo(Coords units) const
{
    switch (units) {
    case Coords::RADIANS:
        return [](double value) { return value; };
    case Coords::DEGREES:
        return [](double value) { return Units::rad2deg(value); };
    case Coords::QSPACE:
        return [wl = m_wavelength](double value) { return getQ(wl, value); };
    case Coords::RQ4:
        return [wl = m_wavelength](double value) { return getQ(wl, value); };
    default:
        ASSERT(0);
    }
}
 
 
//  ************************************************************************************************
//  class WavenumberReflectometryCoordinates
//  ************************************************************************************************
 
WavenumberReflectometryCoordinates::WavenumberReflectometryCoordinates(const IAxis*&& axis)
    : CoordSystem1D(std::move(axis))
{
}
 
WavenumberReflectometryCoordinates::WavenumberReflectometryCoordinates(
    const WavenumberReflectometryCoordinates& other)
    : CoordSystem1D(other.coordinateAxis()->clone())
{
}
 
WavenumberReflectometryCoordinates::~WavenumberReflectometryCoordinates() = default;
 
WavenumberReflectometryCoordinates* WavenumberReflectometryCoordinates::clone() const
{
    return new WavenumberReflectometryCoordinates(*this);
}
 
//! Returns the list of all available units
std::vector<Coords> WavenumberReflectometryCoordinates::availableUnits() const
{
    return {Coords::NBINS, Coords::QSPACE, Coords::RQ4};
}
 
//! Creates name map for axis in various units
std::vector<std::map<Coords, std::string>>
WavenumberReflectometryCoordinates::createNameMaps() const
{
    return {DataUtils::AxisNames::specAxisQ};
}
 
//! Returns translating functional (inv. nm --> desired units)
std::function<double(double)> WavenumberReflectometryCoordinates::getTraslatorTo(Coords units) const
{
    switch (units) {
    case Coords::QSPACE:
        return [](double value) { return value; };
    case Coords::RQ4:
        return [](double value) { return value; };
    default:
        ASSERT(0);
    }
}