Before starting scattering siulation check that refractive index has reasonable magnitude

Sample/HardParticle/Cylinder.cpp

@@ -18,8 +18,6 @@
 #include "Base/Util/Assert.h"
 #include "Sample/Shapes/DoubleEllipse.h"
 #include <numbers>
-#include <sstream>
-#include <stdexcept>
 using std::numbers::pi;
 
 Cylinder::Cylinder(const std::vector<double> P)
@@ -47,15 +45,7 @@ complex_t Cylinder::formfactor(C3 q) const
     const complex_t axial_part = H * Math::sinc(qH2);
     const complex_t radial_part = (2 * pi) * R * R * Math::Bessel::J1c(qR);
 
-    const complex_t result = radial_part * axial_part * exp_I(qH2);
-    if (!std::isfinite(result.real()) || !std::isfinite(result.imag())) {
-        std::stringstream msg;
-        msg << "Cylinder has infinite form factor " << result << " at q=" << q << " qH2=" << qH2
-            << " qR=" << qR << "; factors: radial=" << radial_part << " axial=" << axial_part
-            << " shift=" << exp_I(qH2);
-        throw std::runtime_error(msg.str());
-    }
-    return result;
+    return radial_part * axial_part * exp_I(qH2);
 }
 
 std::string Cylinder::validate() const

Sample/Material/Material.cpp

@@ -17,6 +17,8 @@
 #include "Base/Util/Assert.h"
 #include "Base/Vector/WavevectorInfo.h"
 #include <cmath>
+#include <sstream>
+#include <stdexcept>
 #include <typeinfo>
 
 Material::Material(std::unique_ptr<IMaterialImpl>&& material_impl)
@@ -75,6 +77,17 @@ MATERIAL_TYPES Material::typeID() const
     return m_material_impl->typeID();
 }
 
+void Material::checkRefractiveIndex(double wavelength) const
+{
+    const complex_t n = refractiveIndex(wavelength);
+    if (n.real() < 0.9 || n.real() > 1.1) {
+	std::stringstream msg;
+	msg << "Refractive index " << n << " at wavelength " << wavelength
+	    << " is too far from 1. Invalid material data?";
+	throw std::runtime_error(msg.str());
+    }
+}
+
 R3 Material::magnetization() const
 {
     return m_material_impl->magnetization();

Sample/Material/Material.h

@@ -69,6 +69,8 @@ public:
 #ifndef SWIG
     //! Returns the type of underlying material implementation
     MATERIAL_TYPES typeID() const;
+
+    void checkRefractiveIndex(double wavelength) const;
 #endif // SWIG
 
     //! Get the magnetization (in A/m)

Sample/Multilayer/MultiLayer.cpp

@@ -46,6 +46,12 @@ MultiLayer* MultiLayer::clone() const
     return result;
 }
 
+void MultiLayer::checkMaterials(double wavelength) const
+{
+    for (size_t i = 0; i < numberOfLayers(); ++i)
+	m_layers[i]->material()->checkRefractiveIndex(wavelength);
+}
+
 //! Adds layer with default (zero) roughness
 void MultiLayer::addLayer(const Layer& layer)
 {

Sample/Multilayer/MultiLayer.h

@@ -89,6 +89,8 @@ public:
         return m_ext_field;
     } //!< (units: A/m)
 
+    void checkMaterials(double wavelength) const;
+
 private:
     void addLayerExec(const Layer& layer, const LayerRoughness* roughness);
 

Sim/Simulation/ScatteringSimulation.cpp

@@ -25,6 +25,7 @@
 #include "Device/Detector/SimulationAreaIterator.h" // beginNonMaskedPoints
 #include "Param/Distrib/DistributionHandler.h"
 #include "Resample/Element/DiffuseElement.h"
+#include "Sample/Multilayer/MultiLayer.h"
 #include "Sim/Background/IBackground.h"
 #include "Sim/Computation/DWBAComputation.h"
 
@@ -71,6 +72,7 @@ void ScatteringSimulation::initDistributionHandler()
 
 void ScatteringSimulation::prepareSimulation()
 {
+    sample()->checkMaterials(beam().wavelength());
     m_active_indices = m_detector->active_indices();
     m_pixels.reserve(m_active_indices.size());
     for (auto detector_index : m_active_indices)