diff --git a/Core/Computation/ProcessedLayout.cpp b/Core/Computation/ProcessedLayout.cpp
index 1a8947691810df66176382414c69c3298070ae63..045e881448c39562c0501a8ec5a8bd28548fc9ab 100644
--- a/Core/Computation/ProcessedLayout.cpp
+++ b/Core/Computation/ProcessedLayout.cpp
@@ -81,7 +81,7 @@ void ProcessedLayout::collectFormFactors(const ILayout& layout, const std::vecto
 {
     double layout_abundance = layout.getTotalAbundance();
     for (auto p_particle : layout.particles()) {
-        auto ff_coh = ProcessParticle(*p_particle, slices, z_ref);
+        auto ff_coh = processParticle(*p_particle, slices, z_ref);
         ff_coh.scaleRelativeAbundance(layout_abundance);
         m_formfactors.push_back(std::move(ff_coh));
     }
@@ -91,7 +91,7 @@ void ProcessedLayout::collectFormFactors(const ILayout& layout, const std::vecto
     ScaleRegionMap(m_region_map, scale_factor);
 }
 
-FormFactorCoherentSum ProcessedLayout::ProcessParticle(const IParticle& particle,
+FormFactorCoherentSum ProcessedLayout::processParticle(const IParticle& particle,
                                                        const std::vector<Slice>& slices,
                                                        double z_ref)
 {
diff --git a/Core/Computation/ProcessedLayout.h b/Core/Computation/ProcessedLayout.h
index e4071c55998140a38e6783ff4e90b92bced09081..54c01e3aa1f0619aa16b8e4c3f6542fb4543bd08 100644
--- a/Core/Computation/ProcessedLayout.h
+++ b/Core/Computation/ProcessedLayout.h
@@ -50,7 +50,7 @@ public:
 
 private:
     void collectFormFactors(const ILayout& layout, const std::vector<Slice>& slices, double z_ref);
-    FormFactorCoherentSum ProcessParticle(const IParticle& particle,
+    FormFactorCoherentSum processParticle(const IParticle& particle,
                                           const std::vector<Slice>& slices, double z_ref);
     void mergeRegionMap(const std::map<size_t, std::vector<HomogeneousRegion>>& region_map);
     const IFresnelMap* m_fresnel_map;
diff --git a/Device/Histo/SimulationResult.cpp b/Device/Histo/SimulationResult.cpp
index 6fab60ed2a726085f2132d3243bea1780ddf5a38..5108a31b2ef35f8d2d706b778b7d1d844e49e23c 100644
--- a/Device/Histo/SimulationResult.cpp
+++ b/Device/Histo/SimulationResult.cpp
@@ -89,26 +89,20 @@ std::vector<AxisInfo> SimulationResult::axisInfo(Axes::Units units) const
 
 const IUnitConverter& SimulationResult::converter() const
 {
-    if (!m_unit_converter)
-        throw std::runtime_error(
-            "Error in SimulationResult::converter: unit converter was not initialized");
+    ASSERT(m_unit_converter);
     return *m_unit_converter;
 }
 
 double& SimulationResult::operator[](size_t i)
 {
-    if (m_data)
-        return (*m_data)[i];
-    throw std::runtime_error("Error in SimulationResult::operator[]: "
-                             "no data initialized");
+    ASSERT(m_data);
+    return (*m_data)[i];
 }
 
 const double& SimulationResult::operator[](size_t i) const
 {
-    if (m_data)
-        return (*m_data)[i];
-    throw std::runtime_error("Error in SimulationResult::operator[]: "
-                             "no data initialized");
+    ASSERT(m_data);
+    return (*m_data)[i];
 }
 
 size_t SimulationResult::size() const
@@ -116,6 +110,16 @@ size_t SimulationResult::size() const
     return m_data ? m_data->getAllocatedSize() : 0;
 }
 
+double SimulationResult::max() const
+{
+    ASSERT(m_data);
+    double result = 0;
+    for (size_t i=0; i<size(); ++i)
+        if ((*m_data)[i]>result)
+            result = (*m_data)[i];
+    return result;
+}
+
 #ifdef BORNAGAIN_PYTHON
 PyObject* SimulationResult::array(Axes::Units units) const
 {
diff --git a/Device/Histo/SimulationResult.h b/Device/Histo/SimulationResult.h
index dee399c5d593150de6f4c23b49e5bc08e9432d3c..d5011cb88841acf02016e524708f65fecf619bfe 100644
--- a/Device/Histo/SimulationResult.h
+++ b/Device/Histo/SimulationResult.h
@@ -44,9 +44,7 @@ public:
     SimulationResult(const OutputData<double>& data, const IUnitConverter& unit_converter);
 
     SimulationResult(const SimulationResult& other);
-#ifndef SWIG
     SimulationResult(SimulationResult&& other);
-#endif
 
     SimulationResult& operator=(const SimulationResult& other);
     SimulationResult& operator=(SimulationResult&& other);
@@ -67,6 +65,7 @@ public:
     double& operator[](size_t i);
     const double& operator[](size_t i) const;
     size_t size() const;
+    double max() const;
     bool empty() const { return size() == 0; }
 
     //! returns intensity data as Python numpy array
diff --git a/Sample/Interference/IInterferenceFunctionStrategy.cpp b/Sample/Interference/IInterferenceFunctionStrategy.cpp
index 89cc44c9af0a52caa79ede41fdd599ae7af04752..8e48c4c7fc18ceaf32b1d7b51beeaf540d74b62d 100644
--- a/Sample/Interference/IInterferenceFunctionStrategy.cpp
+++ b/Sample/Interference/IInterferenceFunctionStrategy.cpp
@@ -15,6 +15,7 @@
 #include "Sample/Interference/IInterferenceFunctionStrategy.h"
 #include "Base/Pixel/SimulationElement.h"
 #include "Base/Types/Exceptions.h"
+#include "Base/Utils/Assert.h"
 #include "Base/Utils/IntegratorMCMiser.h"
 #include "Sample/Aggregate/InterferenceFunctionNone.h"
 #include "Sample/Fresnel/FormFactorCoherentSum.h"
@@ -22,18 +23,13 @@
 IInterferenceFunctionStrategy::IInterferenceFunctionStrategy(
     const std::vector<FormFactorCoherentSum>& weighted_formfactors,
     const IInterferenceFunction* p_iff, const SimulationOptions& sim_params, bool polarized)
-    : m_iff(nullptr), m_options(sim_params), m_polarized(polarized),
+    : m_formfactor_wrappers(weighted_formfactors),
+      m_iff(p_iff ? p_iff->clone() : new InterferenceFunctionNone()), m_options(sim_params),
+      m_polarized(polarized),
       m_integrator(
           make_integrator_miser(this, &IInterferenceFunctionStrategy::evaluate_for_fixed_angles, 2))
 {
-    if (weighted_formfactors.empty())
-        throw Exceptions::ClassInitializationException(
-            "IInterferenceFunctionStrategy::init: strategy gets no form factors.");
-    m_formfactor_wrappers = weighted_formfactors;
-    if (p_iff)
-        m_iff.reset(p_iff->clone());
-    else
-        m_iff.reset(new InterferenceFunctionNone());
+    ASSERT(!m_formfactor_wrappers.empty());
 }
 
 IInterferenceFunctionStrategy::~IInterferenceFunctionStrategy() = default;
@@ -50,8 +46,7 @@ IInterferenceFunctionStrategy::evaluateSinglePoint(const SimulationElement& sim_
 {
     if (!m_polarized)
         return scalarCalculation(sim_element);
-    else
-        return polarizedCalculation(sim_element);
+    return polarizedCalculation(sim_element);
 }
 
 //! Performs a Monte Carlo integration over the bin for the evaluation of the intensity.
diff --git a/Sample/Material/MaterialBySLDImpl.h b/Sample/Material/MaterialBySLDImpl.h
index bb62940fe65d48fe1fb8001d78fc7581b1720cc1..bdd97cd85a5382f71c06b8b9667a7dc016bf672e 100644
--- a/Sample/Material/MaterialBySLDImpl.h
+++ b/Sample/Material/MaterialBySLDImpl.h
@@ -59,8 +59,8 @@ private:
     //! Returns the scattering length density
     complex_t sld() const;
 
-    double m_sld_real; //!< complex-valued scattering length density
-    double m_sld_imag; //!< imaginary part of scattering length density (negative by default)
+    const double m_sld_real; //!< complex-valued scattering length density
+    const double m_sld_imag; //!< imaginary part of scattering length density (negative by default)
 };
 
 #endif // BORNAGAIN_SAMPLE_MATERIAL_MATERIALBYSLDIMPL_H
diff --git a/Sample/Material/RefractiveMaterialImpl.cpp b/Sample/Material/RefractiveMaterialImpl.cpp
index d8faf1d9c94cbb6ad444b213b4d65265c6b16f2d..b9295e7029fdadbf79357e08c8d5faef80fb9b12 100644
--- a/Sample/Material/RefractiveMaterialImpl.cpp
+++ b/Sample/Material/RefractiveMaterialImpl.cpp
@@ -48,8 +48,7 @@ complex_t RefractiveMaterialImpl::materialData() const
 complex_t RefractiveMaterialImpl::scalarSubtrSLD(const WavevectorInfo& wavevectors) const
 {
     double wavelength = wavevectors.getWavelength();
-    double prefactor = M_PI / wavelength / wavelength;
-    return prefactor * refractiveIndex2(wavelength);
+    return M_PI / wavelength / wavelength * refractiveIndex2(wavelength);
 }
 
 void RefractiveMaterialImpl::print(std::ostream& ostr) const
diff --git a/Sample/Material/RefractiveMaterialImpl.h b/Sample/Material/RefractiveMaterialImpl.h
index 005cd1684e09932d770964f8082260604750f9c7..201230486eefa452b2462f5eac77d27c925f9cc1 100644
--- a/Sample/Material/RefractiveMaterialImpl.h
+++ b/Sample/Material/RefractiveMaterialImpl.h
@@ -60,9 +60,10 @@ private:
     RefractiveMaterialImpl(const std::string& name, double delta, double beta,
                            kvector_t magnetization);
 
-    double
-        m_delta;   //!< \f$\delta\f$ coefficient for refractive index \f$n = 1 - \delta + i \beta\f$
-    double m_beta; //!< \f$\beta\f$ coefficient for refractive index \f$n = 1 - \delta + i \beta\f$
+    const double
+        m_delta; //!< \f$\delta\f$ coefficient for refractive index \f$n = 1 - \delta + i \beta\f$
+    const double
+        m_beta; //!< \f$\beta\f$ coefficient for refractive index \f$n = 1 - \delta + i \beta\f$
 };
 
 #endif // BORNAGAIN_SAMPLE_MATERIAL_REFRACTIVEMATERIALIMPL_H
diff --git a/auto/Wrap/doxygenDevice.i b/auto/Wrap/doxygenDevice.i
index a5bdaf8d91d8103424946a7530a74b51cb9d97c8..0896fd73285617fc92988420c74bf8191db50fbc 100644
--- a/auto/Wrap/doxygenDevice.i
+++ b/auto/Wrap/doxygenDevice.i
@@ -2870,6 +2870,9 @@ Returns underlying unit converter.
 %feature("docstring")  SimulationResult::size "size_t SimulationResult::size() const
 ";
 
+%feature("docstring")  SimulationResult::max "double & SimulationResult::max() const
+";
+
 %feature("docstring")  SimulationResult::empty "bool SimulationResult::empty() const
 ";
 
diff --git a/auto/Wrap/libBornAgainDevice.py b/auto/Wrap/libBornAgainDevice.py
index 73e00a8ef26b5711a22e35ac5803cb36b0df0703..6a4529adb44be29552a1703d0a1bb676eaaf2a74 100644
--- a/auto/Wrap/libBornAgainDevice.py
+++ b/auto/Wrap/libBornAgainDevice.py
@@ -6082,6 +6082,7 @@ class SimulationResult(object):
         __init__(SimulationResult self) -> SimulationResult
         __init__(SimulationResult self, IntensityData data, IUnitConverter unit_converter) -> SimulationResult
         __init__(SimulationResult self, SimulationResult other) -> SimulationResult
+        __init__(SimulationResult self, SimulationResult other) -> SimulationResult
         SimulationResult::SimulationResult(SimulationResult &&other)
 
         """
@@ -6123,6 +6124,14 @@ class SimulationResult(object):
         """
         return _libBornAgainDevice.SimulationResult_size(self)
 
+    def max(self):
+        r"""
+        max(SimulationResult self) -> double
+        double & SimulationResult::max() const
+
+        """
+        return _libBornAgainDevice.SimulationResult_max(self)
+
     def empty(self):
         r"""
         empty(SimulationResult self) -> bool
diff --git a/auto/Wrap/libBornAgainDevice_wrap.cpp b/auto/Wrap/libBornAgainDevice_wrap.cpp
index ef74f231d2158e04ef12f1b037fd360640c03b3b..cf6146414700ecf6c3a811ff41b7129b4732df12 100644
--- a/auto/Wrap/libBornAgainDevice_wrap.cpp
+++ b/auto/Wrap/libBornAgainDevice_wrap.cpp
@@ -44907,6 +44907,30 @@ fail:
 }
 
 
+SWIGINTERN PyObject *_wrap_new_SimulationResult__SWIG_3(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
+  PyObject *resultobj = 0;
+  SimulationResult *arg1 = 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  SimulationResult *result = 0 ;
+  
+  if ((nobjs < 1) || (nobjs > 1)) SWIG_fail;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1, SWIGTYPE_p_SimulationResult,  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "new_SimulationResult" "', argument " "1"" of type '" "SimulationResult &&""'"); 
+  }
+  if (!argp1) {
+    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "new_SimulationResult" "', argument " "1"" of type '" "SimulationResult &&""'"); 
+  }
+  arg1 = reinterpret_cast< SimulationResult * >(argp1);
+  result = (SimulationResult *)new SimulationResult((SimulationResult &&)*arg1);
+  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_SimulationResult, SWIG_POINTER_NEW |  0 );
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
 SWIGINTERN PyObject *_wrap_new_SimulationResult(PyObject *self, PyObject *args) {
   Py_ssize_t argc;
   PyObject *argv[3] = {
@@ -44926,6 +44950,15 @@ SWIGINTERN PyObject *_wrap_new_SimulationResult(PyObject *self, PyObject *args)
       return _wrap_new_SimulationResult__SWIG_2(self, argc, argv);
     }
   }
+  if (argc == 1) {
+    int _v;
+    void *vptr = 0;
+    int res = SWIG_ConvertPtr(argv[0], &vptr, SWIGTYPE_p_SimulationResult, SWIG_POINTER_NO_NULL);
+    _v = SWIG_CheckState(res);
+    if (_v) {
+      return _wrap_new_SimulationResult__SWIG_3(self, argc, argv);
+    }
+  }
   if (argc == 2) {
     int _v;
     int res = SWIG_ConvertPtr(argv[0], 0, SWIGTYPE_p_OutputDataT_double_t, SWIG_POINTER_NO_NULL | 0);
@@ -44944,7 +44977,8 @@ fail:
     "  Possible C/C++ prototypes are:\n"
     "    SimulationResult::SimulationResult()\n"
     "    SimulationResult::SimulationResult(OutputData< double > const &,IUnitConverter const &)\n"
-    "    SimulationResult::SimulationResult(SimulationResult const &)\n");
+    "    SimulationResult::SimulationResult(SimulationResult const &)\n"
+    "    SimulationResult::SimulationResult(SimulationResult &&)\n");
   return 0;
 }
 
@@ -45179,6 +45213,29 @@ fail:
 }
 
 
+SWIGINTERN PyObject *_wrap_SimulationResult_max(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+  PyObject *resultobj = 0;
+  SimulationResult *arg1 = (SimulationResult *) 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  PyObject *swig_obj[1] ;
+  double result;
+  
+  if (!args) SWIG_fail;
+  swig_obj[0] = args;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_SimulationResult, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "SimulationResult_max" "', argument " "1"" of type '" "SimulationResult const *""'"); 
+  }
+  arg1 = reinterpret_cast< SimulationResult * >(argp1);
+  result = (double)((SimulationResult const *)arg1)->max();
+  resultobj = SWIG_From_double(static_cast< double >(result));
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
 SWIGINTERN PyObject *_wrap_SimulationResult_empty(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
   SimulationResult *arg1 = (SimulationResult *) 0 ;
@@ -48625,6 +48682,7 @@ static PyMethodDef SwigMethods[] = {
 	 { "new_SimulationResult", _wrap_new_SimulationResult, METH_VARARGS, "\n"
 		"SimulationResult()\n"
 		"SimulationResult(IntensityData data, IUnitConverter unit_converter)\n"
+		"SimulationResult(SimulationResult other)\n"
 		"new_SimulationResult(SimulationResult other) -> SimulationResult\n"
 		"SimulationResult::SimulationResult(SimulationResult &&other)\n"
 		"\n"
@@ -48653,6 +48711,11 @@ static PyMethodDef SwigMethods[] = {
 		"size_t SimulationResult::size() const\n"
 		"\n"
 		""},
+	 { "SimulationResult_max", _wrap_SimulationResult_max, METH_O, "\n"
+		"SimulationResult_max(SimulationResult self) -> double\n"
+		"double & SimulationResult::max() const\n"
+		"\n"
+		""},
 	 { "SimulationResult_empty", _wrap_SimulationResult_empty, METH_O, "\n"
 		"SimulationResult_empty(SimulationResult self) -> bool\n"
 		"bool SimulationResult::empty() const\n"