[j.715a] Before starting scattering simulation check that refractive index has...

[j.715a] Before starting scattering simulation check that refractive index has reasonable magnitude (#715) (Closes #715)

Merging branch 'j.715a'  into 'main'.

See merge request !2001

CHANGELOG

@@ -21,6 +21,8 @@ BornAgain-22.0, in preparation
     * Simplified coordinate names and units, now both in Scale::name (#487, #614)
     * Remove ICoordSystem and children; functionality provided by Frame
     * Merge SimulationResult into Datafield
+    * In scattering simulation, restrict Re(n) to [0.9,1.1], to prevent unintended
+      computation with erroneous material data (e.g. accidental omission of factor 1e-6) (#715)
 
 BornAgain-21.0, released 2023.08.16
   > New functionality:

Device/IO/DiffUtil.cpp

@@ -30,7 +30,7 @@ double DiffUtil::meanRelVecDiff(const std::vector<double>& dat, const std::vecto
     for (size_t i = 0; i < dat.size(); ++i)
         diff += Numeric::relativeDifference(dat[i], ref[i]);
     diff /= dat.size();
-    ASSERT(!std::isnan(diff));
+    ASSERT(std::isfinite(diff));
     return diff;
 }
 

Resample/Interparticle/DecouplingApproximationStrategy.cpp

@@ -34,7 +34,7 @@ double DecouplingApproximationStrategy::scalarCalculation(const DiffuseElement&
     complex_t amplitude = 0;
     for (const auto& ffw : m_weighted_formfactors) {
         const complex_t ff = ffw->summedFF(ele);
-        ASSERT(!std::isnan(ff.real())); // numerical error in coherent sum?
+        ASSERT(std::isfinite(ff.real())); // numerical error in coherent sum?
         const double fraction = ffw->relativeAbundance();
         amplitude += fraction * ff;
         intensity += fraction * std::norm(ff);

Resample/Slice/Slice.cpp

@@ -13,13 +13,13 @@
 //  ************************************************************************************************
 
 #include "Resample/Slice/Slice.h"
-#include <numbers>
-using std::numbers::pi;
 #include "Base/Spin/SpinMatrix.h"
 #include "Resample/Slice/SliceStack.h"
 #include "Sample/Interface/LayerRoughness.h"
 #include "Sample/Material/MaterialUtil.h"
+#include <numbers>
 #include <utility>
+using std::numbers::pi;
 
 Slice::Slice(const ZLimits& zRange, Material material, const R3& B_field,
              const LayerRoughness* const roughness)

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);
 

Sample/Particle/IFormFactor.cpp

@@ -21,6 +21,8 @@
 #include "Sample/Particle/PolyhedralUtil.h"
 #include "Sample/Scattering/Rotations.h"
 #include "Sample/Shapes/IShape3D.h"
+#include <sstream>
+#include <stdexcept>
 
 IFormFactor::IFormFactor() = default;
 
@@ -58,7 +60,14 @@ std::string IFormFactor::pythonConstructor() const
 
 complex_t IFormFactor::theFF(const WavevectorInfo& wavevectors) const
 {
-    return formfactor(wavevectors.getQ());
+    const complex_t result = formfactor(wavevectors.getQ());
+    if (!std::isfinite(result.real()) || !std::isfinite(result.imag())) {
+        std::stringstream msg;
+        msg << "Infinite form factor " << result << " at ki " << wavevectors.getKi() << " kf "
+            << wavevectors.getKf() << " q " << wavevectors.getQ();
+        throw std::runtime_error(msg.str());
+    }
+    return result;
 }
 
 SpinMatrix IFormFactor::thePolFF(const WavevectorInfo& wavevectors) const

Sample/StandardSamples/ThickAbsorptiveSampleBuilder.h

@@ -20,7 +20,6 @@
 
 class MultiLayer;
 
-
 namespace ExemplarySamples {
 
 MultiLayer* createThickAbsorptiveSample();

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)

Tests/SimFactory/MakeSimulations.cpp

@@ -270,7 +270,7 @@ std::unique_ptr<ScatteringSimulation> test::makeSimulation::RectDetWithRoi(const
 std::unique_ptr<ScatteringSimulation>
 test::makeSimulation::ExtraLongWavelengthGISAS(const MultiLayer& sample)
 {
-    std::unique_ptr<Beam> beam(Beam(1e8, 13.52 * Units::nm, 0.2 * deg).clone());
+    std::unique_ptr<Beam> beam(Beam(1e8, 12 * Units::nm, 0.2 * deg).clone());
     FlatDetector detector(100, 100, 70, 70, stdBeam, FlatDetector::X, 2000, 0, -35);
     auto simulation = std::make_unique<ScatteringSimulation>(*beam, sample, detector);
     simulation->options().setIncludeSpecular(true);