fix/ue112-rml-materials: RML import and material tables

.gitignore

@@ -28,6 +28,7 @@ CMakeCache.txm
 cmake.check_cache
 imgui.ini
 *.h5
+*.testCsv.csv
 
 # VSCode Debugger files
 results.json

Extern/CMakeLists.txt

@@ -103,5 +103,26 @@ if(RAYX_BUILD_RAYX_UI)
     add_library(portable_file_dialogs INTERFACE)
     target_include_directories(portable_file_dialogs INTERFACE ${CMAKE_CURRENT_SOURCE_DIR}/portable-file-dialogs)
 
+    # rayx-ui only needs SDL for the window, input and the Vulkan surface. The imgui backend
+    # also polls gamepads, so joystick (and HIDAPI, which SDL 3.2 needs to link joystick) stays on.
+    set(SDL_AUDIO OFF)
+    set(SDL_CAMERA OFF)
+    set(SDL_GPU OFF)
+    set(SDL_RENDER OFF)
+    set(SDL_HAPTIC OFF)
+    set(SDL_POWER OFF)
+    set(SDL_SENSOR OFF)
+    set(SDL_DIALOG OFF)
+    set(SDL_TEST_LIBRARY OFF)
+
+    # hide SDL's configure checks and option summary; warnings and errors still show.
+    # `cmake --log-level=STATUS` takes precedence and brings them back.
+    set(CMAKE_MESSAGE_LOG_LEVEL NOTICE)
     add_subdirectory(SDL EXCLUDE_FROM_ALL)
+    unset(CMAKE_MESSAGE_LOG_LEVEL)
+
+    # third-party code, its warnings are not ours to fix
+    if(TARGET SDL3-shared AND NOT MSVC)
+        target_compile_options(SDL3-shared PRIVATE -w)
+    endif()
 endif()

Intern/rayx-core/src/Beamline/Beamline.cpp

@@ -1,12 +1,17 @@
 #include "Beamline.h"
 
+#include <algorithm>
+#include <limits>
+#include <optional>
 #include <sstream>
 #include <stack>
 #include <stdexcept>
+#include <type_traits>
 
 #include "Debug/Instrumentor.h"
 #include "Design/DesignElement.h"
 #include "Design/DesignSource.h"
+#include "Shader/RefractiveIndex.h"
 
 namespace rayx {
 
@@ -266,21 +271,15 @@ MaterialTables Group::calcMinimalMaterialTables() const {
         auto coating = elemPtr->getCoating();
         if (coating.is<detail::CoatingTypes::OneCoating>()) {
             int materialcoating = static_cast<int>(elemPtr->getMaterialCoating());
-            if (materialcoating >= 1 && materialcoating <= 133) {
-                relevantMaterials[materialcoating - 1] = true;
-            }
+            if (materialcoating >= 1 && materialcoating <= 133) { relevantMaterials[materialcoating - 1] = true; }
         } else if (coating.is<detail::CoatingTypes::MultilayerCoating>()) {
             auto mlCoating = coating.get<detail::CoatingTypes::MultilayerCoating>();
             for (const auto& mat : mlCoating.material) {
-                if (mat >= 1 && mat <= 133) {
-                    relevantMaterials[mat - 1] = true;
-                }
+                if (mat >= 1 && mat <= 133) { relevantMaterials[mat - 1] = true; }
             }
         }
         int material = static_cast<int>(elemPtr->getMaterial());  // assuming getMaterial() exists
-        if (material >= 1 && material <= 133) {
-            relevantMaterials[material - 1] = true;
-        }
+        if (material >= 1 && material <= 133) { relevantMaterials[material - 1] = true; }
     }
     return loadMaterialTables(relevantMaterials);
 }
@@ -327,6 +326,127 @@ std::vector<OpticalElementAndTransform> Group::compileElements() const {
     return elements;
 }
 
+namespace {
+
+/// returns the [min, max] energy covered by one energy distribution, if bounded.
+std::optional<std::pair<double, double>> energyRangeOf(const EnergyDistributionVariant& en) {
+    return std::visit(
+        [](const auto& value) -> std::optional<std::pair<double, double>> {
+            using T = std::decay_t<decltype(value)>;
+            if constexpr (std::is_same_v<T, DatFile>) {
+                if (value.m_Lines.empty()) return std::nullopt;
+                double lo = std::numeric_limits<double>::max();
+                double hi = std::numeric_limits<double>::lowest();
+                for (const auto& e : value.m_Lines) {
+                    lo = std::min(lo, e.m_energy);
+                    hi = std::max(hi, e.m_energy);
+                }
+                return std::make_pair(lo, hi);
+            } else if constexpr (std::is_same_v<T, HardEdge>) {
+                return std::make_pair(value.m_centerEnergy - value.m_energySpread / 2.0, value.m_centerEnergy + value.m_energySpread / 2.0);
+            } else if constexpr (std::is_same_v<T, SeparateEnergies>) {
+                return std::make_pair(value.m_centerEnergy - value.m_energySpread / 2.0, value.m_centerEnergy + value.m_energySpread / 2.0);
+            } else if constexpr (std::is_same_v<T, SoftEdge>) {
+                // a Gaussian is technically unbounded; use the central +/- 5 sigma window
+                const double sigma = value.m_sigma > 0 ? value.m_sigma : 1.0;
+                return std::make_pair(value.m_centerEnergy - 5.0 * sigma, value.m_centerEnergy + 5.0 * sigma);
+            }
+            return std::nullopt;
+        },
+        en);
+}
+
+/// whether the material tables of `material` cover the whole [lo, hi] energy range.
+bool materialCoversRange(int material, const MaterialTables& tables, double lo, double hi) {
+    if (material < 1 || material > 140) return true;  // vacuum/ideal reflectors need no data
+    const int* idx    = tables.indices.data();
+    const double* tab = tables.materials.data();
+
+    if (material <= 92) {
+        if (getPalikEntryCount(material, idx) > 0) {
+            auto first = getPalikEntry(0, material, idx, tab);
+            auto last  = getPalikEntry(getPalikEntryCount(material, idx) - 1, material, idx, tab);
+            if (first.m_energy <= lo && hi <= last.m_energy) return true;
+        }
+        if (getNffEntryCount(material, idx) > 0) {
+            auto first = getNffEntry(0, material, idx, tab);
+            auto last  = getNffEntry(getNffEntryCount(material, idx) - 1, material, idx, tab);
+            if (first.m_energy <= lo && hi <= last.m_energy) return true;
+        }
+        if (getCromerEntryCount(material, idx) > 0) {
+            auto first = getCromerEntry(0, material, idx, tab);
+            auto last  = getCromerEntry(getCromerEntryCount(material, idx) - 1, material, idx, tab);
+            if (first.m_energy <= lo && hi <= last.m_energy) return true;
+        }
+    } else {
+        if (getMolecEntryCount(material, idx) > 0) {
+            auto first = getMolecEntry(0, material, idx, tab);
+            auto last  = getMolecEntry(getMolecEntryCount(material, idx) - 1, material, idx, tab);
+            if (first.m_energy <= lo && hi <= last.m_energy) return true;
+        }
+    }
+    return false;
+}
+
+/// all materials that the given compiled element may perform lookups with.
+std::vector<int> materialsOf(const OpticalElement& element) {
+    std::vector<int> mats;
+
+    // only behaveMirror and behaveFoil call getRefractiveIndex; every other behaviour
+    // ignores the material entirely
+    if (!element.m_behaviour.is<Behaviour::Mirror>() && !element.m_behaviour.is<Behaviour::Foil>()) return mats;
+
+    if (element.m_material >= 1 && element.m_material <= 140) mats.push_back(element.m_material);
+
+    element.m_coating.visit([&](const auto& coating) {
+        using T = std::decay_t<decltype(coating)>;
+        if constexpr (std::is_same_v<T, Coating::OneCoating>) {
+            if (coating.material >= 1 && coating.material <= 140) mats.push_back(coating.material);
+        } else if constexpr (std::is_same_v<T, Coating::MultilayerCoating>) {
+            for (int i = 0; i < coating.numLayers; ++i) {
+                if (coating.material[i] >= 1 && coating.material[i] <= 140) mats.push_back(coating.material[i]);
+            }
+        }
+    });
+    return mats;
+}
+
+}  // unnamed namespace
+
+void Group::verifyMaterialCoverage(std::vector<OpticalElementAndTransform>& compiledElements, const MaterialTables& tables) const {
+    double energyLo = std::numeric_limits<double>::max();
+    double energyHi = std::numeric_limits<double>::lowest();
+    for (const auto* source : getSources()) {
+        // GenRays skips these too
+        if (source->getType() == ElementType::RayListSource) continue;
+        auto range = energyRangeOf(source->getEnergyDistribution());
+        if (!range) continue;
+        energyLo = std::min(energyLo, range->first);
+        energyHi = std::max(energyHi, range->second);
+    }
+    if (energyLo > energyHi) return;
+
+    for (size_t i = 0; i < compiledElements.size(); i++) {
+        auto& element = compiledElements[i].element;
+        bool covered  = true;
+        for (const auto material : materialsOf(element)) {
+            if (!materialCoversRange(material, tables, energyLo, energyHi)) {
+                covered = false;
+                break;
+            }
+        }
+        if (covered) continue;
+
+        bool isFoil = element.m_behaviour.is<Behaviour::Foil>();
+        RAYX_WARN << "material data for element " << i << " (" << (isFoil ? "foil" : "element") << ") does not cover the source energy range ["
+                  << energyLo << ", " << energyHi << "] eV; compiling it with ideal geometric behaviour "
+                  << (isFoil ? "(100% transmission)" : "(100% reflection)") << " instead";
+
+        element.m_material = static_cast<int>(Material::REFLECTIVE);
+        element.m_coating  = Coating::SubstrateOnly{};
+    }
+}
+
 std::vector<const DesignElement*> Group::getElements() const {
     std::vector<const DesignElement*> elements;
     ctraverse([&elements](const BeamlineNode& node) -> bool {

Intern/rayx-core/src/Beamline/Beamline.h

@@ -138,6 +138,10 @@ class RAYX_API Group : public BeamlineNode {
      */
     std::vector<OpticalElementAndTransform> compileElements() const;
 
+    /// Replaces elements whose materials lack table data over the source energy range with
+    /// ideal geometric behaviour. Only the passed compiled copy is modified.
+    void verifyMaterialCoverage(std::vector<OpticalElementAndTransform>& compiledElements, const MaterialTables& tables) const;
+
     // TODO: why would we need this? ray-ui uses this function
     /**
      * @brief Gathers the world positions of all light sources within a Group hierarchy.

Intern/rayx-core/src/Beamline/DatFile.h

@@ -26,8 +26,6 @@ struct DatFile {
 
     std::vector<DatEntry> m_Lines;
 
-    bool m_continuous;
-
     /** loads the .DAT file `filename` and writes it's contents to `out` */
     static bool load(const std::filesystem::path& filename, DatFile* out);
 

Intern/rayx-core/src/Beamline/Definitions.h

@@ -17,7 +17,7 @@ enum class CurvatureType {
     Quadric,
     RzpSphere
 };  // order is crucial for xml prser
-enum class BehaviourType { Mirror, Grating, Slit, Rzp, ImagePlane, Crystal, Foil };
+enum class BehaviourType { Mirror, Grating, Slit, Rzp, ImagePlane, Crystal, Foil, Count /* must stay last */ };
 enum class FigureRotation { Yes, Plane, A11 };
 
 // the direction of a plane, either XY or XZ. This is only used in the parsing.

Intern/rayx-core/src/Beamline/StringConversion.h

@@ -77,16 +77,12 @@ const std::map<std::string, CurvatureType> StringToCurvatureType = {
     {"Quadric", CurvatureType::Quadric},     {"RzpSphere", CurvatureType::RzpSphere}};
 
 // BehaviourType conversion
-const std::map<BehaviourType, std::string> BehaviourTypeToString = {{BehaviourType::Mirror, "Mirror"},
-                                                                    {BehaviourType::Grating, "Grating"},
-                                                                    {BehaviourType::Slit, "Slit"},
-                                                                    {BehaviourType::Rzp, "Rzp"},
-                                                                    {BehaviourType::ImagePlane, "ImagePlane"}};
-const std::map<std::string, BehaviourType> StringToBehaviourType = {{"Mirror", BehaviourType::Mirror},
-                                                                    {"Grating", BehaviourType::Grating},
-                                                                    {"Slit", BehaviourType::Slit},
-                                                                    {"Rzp", BehaviourType::Rzp},
-                                                                    {"ImagePlane", BehaviourType::ImagePlane}};
+const std::map<BehaviourType, std::string> BehaviourTypeToString = {
+    {BehaviourType::Mirror, "Mirror"},         {BehaviourType::Grating, "Grating"}, {BehaviourType::Slit, "Slit"}, {BehaviourType::Rzp, "Rzp"},
+    {BehaviourType::ImagePlane, "ImagePlane"}, {BehaviourType::Crystal, "Crystal"}, {BehaviourType::Foil, "Foil"}};
+const std::map<std::string, BehaviourType> StringToBehaviourType = {
+    {"Mirror", BehaviourType::Mirror},         {"Grating", BehaviourType::Grating}, {"Slit", BehaviourType::Slit}, {"Rzp", BehaviourType::Rzp},
+    {"ImagePlane", BehaviourType::ImagePlane}, {"Crystal", BehaviourType::Crystal}, {"Foil", BehaviourType::Foil}};
 // FigureRotation conversion
 const std::map<FigureRotation, std::string> FigureRotationToString = {
     {FigureRotation::Yes, "Yes"}, {FigureRotation::Plane, "Plane"}, {FigureRotation::A11, "A11"}};
@@ -139,7 +135,9 @@ const std::map<ElementType, std::string> ElementTypeToString = {{ElementType::Ci
                                                                 {ElementType::PixelSource, "Pixel Source"},
                                                                 {ElementType::EllipsoidMirror, "Ellipsoid"},
                                                                 {ElementType::Crystal, "Crystal"},
-                                                                {ElementType::Foil, "Foil"}};
+                                                                {ElementType::Foil, "Foil"},
+                                                                {ElementType::RayListSource, "Ray List Source"},
+                                                                {ElementType::Undefined, "Undefined"}};
 const std::map<std::string, ElementType> StringToElementType = {{"Circle Source", ElementType::CircleSource},
                                                                 {"Cylinder", ElementType::CylinderMirror},
                                                                 {"ImagePlane", ElementType::ImagePlane},
@@ -163,7 +161,9 @@ const std::map<std::string, ElementType> StringToElementType = {{"Circle Source"
                                                                 {"Pixel Source", ElementType::PixelSource},
                                                                 {"Ellipsoid", ElementType::EllipsoidMirror},
                                                                 {"Crystal", ElementType::Crystal},
-                                                                {"Foil", ElementType::Foil}};
+                                                                {"Foil", ElementType::Foil},
+                                                                {"Ray List Source", ElementType::RayListSource},
+                                                                {"Undefined", ElementType::Undefined}};
 
 const std::map<std::string, EventType> StringToEventType = {
     {"Uninitialized", EventType::Uninitialized}, {"Emitted", EventType::Emitted},   {"HitElement", EventType::HitElement},

Intern/rayx-core/src/Design/DesignSource.cpp

@@ -185,18 +185,14 @@ double DesignSource::getEnergySpread() const { return m_elementParameters["energ
 
 void DesignSource::setEnergySpreadUnit(EnergySpreadUnit value) { m_elementParameters["energySpreadUnit"] = value; }
 EnergySpreadUnit DesignSource::getEnergySpreadUnit() const {
-    if (!m_elementParameters.hasKey("energySpreadUnit")) {
-        return EnergySpreadUnit::EU_eV;
-    }
+    if (!m_elementParameters.hasKey("energySpreadUnit")) { return EnergySpreadUnit::EU_eV; }
 
     return m_elementParameters["energySpreadUnit"].as_energySpreadUnit();
 }
 
 double DesignSource::getEnergySpreadInEv() const {
     const double energySpread = getEnergySpread();
-    if (getEnergySpreadUnit() == EnergySpreadUnit::EU_PERCENT) {
-        return getEnergy() * energySpread / 100.0;
-    }
+    if (getEnergySpreadUnit() == EnergySpreadUnit::EU_PERCENT) { return getEnergy() * energySpread / 100.0; }
 
     return energySpread;
 }
@@ -219,7 +215,6 @@ double DesignSource::getPhotonFlux() const { return m_elementParameters["photonF
 
 EnergyDistributionVariant DesignSource::getEnergyDistribution() const {
     EnergyDistributionVariant en;
-    SpreadType spreadType                         = m_elementParameters["energyDistribution"].as_energySpreadType();
     EnergyDistributionType energyDistributionType = m_elementParameters["energyDistributionType"].as_energyDistributionType();
 
     if (energyDistributionType == EnergyDistributionType::File) {
@@ -228,11 +223,11 @@ EnergyDistributionVariant DesignSource::getEnergyDistribution() const {
         DatFile df;
         DatFile::load(filename, &df);
 
-        df.m_continuous = (spreadType == SpreadType::SoftEdge ? true : false);
-        en              = EnergyDistributionVariant(df);
+        en = EnergyDistributionVariant(df);
     } else if (energyDistributionType == EnergyDistributionType::Values) {
-        double photonEnergy = m_elementParameters["energy"].as_double();
-        double energySpread = getEnergySpreadInEv();
+        SpreadType spreadType = m_elementParameters["energyDistribution"].as_energySpreadType();
+        double photonEnergy   = m_elementParameters["energy"].as_double();
+        double energySpread   = getEnergySpreadInEv();
 
         if (spreadType == SpreadType::SoftEdge) {
             if (energySpread == 0) { energySpread = 1; }

Intern/rayx-core/src/Element/Element.h

@@ -43,6 +43,7 @@ enum class ElementType {
     CircleSource,
     SimpleUndulatorSource,
     RayListSource,
+    Count,  // must stay last
 };
 
 /**

Intern/rayx-core/src/Element/Surface.cpp

@@ -244,20 +244,32 @@ Surface makeParaboloid(const DesignElement& dele) {
     y0 = ArmLength * sin1;
     z0 = ArmLength * cos1 * sign;
 
-    a24 = -y0;
-    a34 = -parameterP;
     a44 = pow(y0, 2) - 2 * parameterP * z0 - pow(parameterP, 2);
+
+    // The paraboloid is built with its axis along z, which leaves the surface tilted by the
+    // grazing angle at the pole. Rotate it about x so the surface is tangent to the element's
+    // xz plane there, as makeEllipsoid does. RAY-UI applies the equivalent rotation through its
+    // tangentAngle (= sign * grazingIncAngle) in the element transform instead.
+    const double alpha = -sign * grazingIncAngle.rad;
+    const double sa    = sin(alpha);
+    const double ca    = cos(alpha);
+
+    const double a22 = ca * ca;
+    const double a23 = sa * ca;
+    a24              = -y0 * ca + parameterP * sa;
+    const double a33 = sa * sa;
+    a34              = -y0 * sa - parameterP * ca;
     //---------------------------- Serialization -------------------------------
     return Surface::Quadric{
         .m_icurv = 1,
         .m_a11   = a11,
         .m_a12   = 0,
         .m_a13   = 0,
         .m_a14   = 0,
-        .m_a22   = 1.0,
-        .m_a23   = 0,
+        .m_a22   = a22,
+        .m_a23   = a23,
         .m_a24   = a24,
-        .m_a33   = 0,
+        .m_a33   = a33,
         .m_a34   = a34,
         .m_a44   = a44,
     };