diff --git a/Core/Aggregate/IInterferenceFunction.cpp b/Core/Aggregate/IInterferenceFunction.cpp
index 4e31f7ed11c5486a0de32a48312e8d2ef891ae65..ad8bc3cab94e1b78b3f3d1cf4c286c08dca40df1 100644
--- a/Core/Aggregate/IInterferenceFunction.cpp
+++ b/Core/Aggregate/IInterferenceFunction.cpp
@@ -29,8 +29,6 @@ IInterferenceFunction::IInterferenceFunction(double position_var) : m_position_v
     registerParameter("PositionVariance", &m_position_var).setUnit("nm^2").setNonnegative();
 }
 
-IInterferenceFunction::~IInterferenceFunction() = default;
-
 // Default implementation of evaluate assumes no inner structure
 // It is only to be overriden in case of the presence of such inner structure. See for example
 // InterferenceFunction2DSuperLattice for such a case.
diff --git a/Core/Aggregate/IInterferenceFunction.h b/Core/Aggregate/IInterferenceFunction.h
index 2657a08b7d88ae27f8a4354914a554f226a804d0..a18434af14283ba9fafb42a9fb90498a79ec8e20 100644
--- a/Core/Aggregate/IInterferenceFunction.h
+++ b/Core/Aggregate/IInterferenceFunction.h
@@ -26,7 +26,6 @@ class BA_CORE_API_ IInterferenceFunction : public ISample
 public:
     IInterferenceFunction(const NodeMeta& meta, const std::vector<double>& PValues);
     IInterferenceFunction(double position_var);
-    virtual ~IInterferenceFunction();
 
     virtual IInterferenceFunction* clone() const = 0;
 
@@ -51,13 +50,13 @@ public:
     double DWfactor(kvector_t q) const;
 
 protected:
+    double m_position_var;
+
     //! Calculates the structure factor in the absence of extra inner structure
     double iff_no_inner(const kvector_t q, double outer_iff) const;
 
     //! Calculates the structure factor without Debye-Waller factor
     virtual double iff_without_dw(const kvector_t q) const = 0;
-
-    double m_position_var;
 };
 
 #endif // BORNAGAIN_CORE_AGGREGATE_IINTERFERENCEFUNCTION_H
diff --git a/Core/Aggregate/InterferenceFunction2DLattice.cpp b/Core/Aggregate/InterferenceFunction2DLattice.cpp
index 6877abaef1034e3103368caaee87e8e4c6cfd743..1c16d41a6f0b05a1d8b31914ea99def7c5563113 100644
--- a/Core/Aggregate/InterferenceFunction2DLattice.cpp
+++ b/Core/Aggregate/InterferenceFunction2DLattice.cpp
@@ -16,6 +16,7 @@
 #include "Core/Basics/Exceptions.h"
 #include "Core/Basics/MathConstants.h"
 #include "Core/Parametrization/RealParameter.h"
+#include "Core/Tools/Integrator.h"
 #include <algorithm>
 
 namespace
@@ -123,8 +124,8 @@ double InterferenceFunction2DLattice::iff_without_dw(const kvector_t q) const
     m_qy = q.y();
     if (!m_integrate_xi)
         return interferenceForXi(m_lattice->rotationAngle());
-    return m_integrator.integrate([&](double xi) -> double { return interferenceForXi(xi); }, 0.0,
-                                  M_TWOPI)
+    return RealIntegrator().integrate([&](double xi) -> double { return interferenceForXi(xi); },
+                                      0.0, M_TWOPI)
            / M_TWOPI;
 }
 
diff --git a/Core/Aggregate/InterferenceFunction2DLattice.h b/Core/Aggregate/InterferenceFunction2DLattice.h
index a86dd8c6cdcff69b8fe54bff1120d2759d22db62..b6c9cbdf57653edba4ccbaae84438e1368589dc7 100644
--- a/Core/Aggregate/InterferenceFunction2DLattice.h
+++ b/Core/Aggregate/InterferenceFunction2DLattice.h
@@ -19,7 +19,6 @@
 #include "Core/Correlations/FTDecay1D.h"
 #include "Core/Correlations/FTDecay2D.h"
 #include "Core/Lattice/Lattice2D.h"
-#include "Core/Tools/Integrator.h"
 
 //! Interference function of a 2D lattice.
 //! @ingroup interference
@@ -82,7 +81,6 @@ private:
     int m_na, m_nb; //!< determines the number of reciprocal lattice points to use
     mutable double m_qx;
     mutable double m_qy;
-    mutable RealIntegrator m_integrator;
 };
 
 #endif // BORNAGAIN_CORE_AGGREGATE_INTERFERENCEFUNCTION2DLATTICE_H
diff --git a/Core/Aggregate/InterferenceFunction2DParaCrystal.cpp b/Core/Aggregate/InterferenceFunction2DParaCrystal.cpp
index ce1c267879abbdc0b9e6f48a99cf0f129b64bae1..fd46adaea828dc188c27a03476b08a1ba043ee87 100644
--- a/Core/Aggregate/InterferenceFunction2DParaCrystal.cpp
+++ b/Core/Aggregate/InterferenceFunction2DParaCrystal.cpp
@@ -17,6 +17,7 @@
 #include "Core/Basics/MathConstants.h"
 #include "Core/Parametrization/ParameterPool.h"
 #include "Core/Parametrization/RealParameter.h"
+#include "Core/Tools/Integrator.h"
 #include <limits>
 
 InterferenceFunction2DParaCrystal::InterferenceFunction2DParaCrystal(const Lattice2D& lattice,
@@ -100,8 +101,8 @@ double InterferenceFunction2DParaCrystal::iff_without_dw(const kvector_t q) cons
     m_qy = q.y();
     if (!m_integrate_xi)
         return interferenceForXi(m_lattice->rotationAngle());
-    return m_integrator.integrate([&](double xi) -> double { return interferenceForXi(xi); }, 0.0,
-                                  M_TWOPI)
+    return RealIntegrator().integrate([&](double xi) -> double { return interferenceForXi(xi); },
+                                      0.0, M_TWOPI)
            / M_TWOPI;
 }
 
diff --git a/Core/Aggregate/InterferenceFunction2DParaCrystal.h b/Core/Aggregate/InterferenceFunction2DParaCrystal.h
index 5739bb0f7b0574a1e33b41ba4cf288c00f5ae845..ce7cf10f7a6a6a012c3a0281add15901045fb33c 100644
--- a/Core/Aggregate/InterferenceFunction2DParaCrystal.h
+++ b/Core/Aggregate/InterferenceFunction2DParaCrystal.h
@@ -19,7 +19,6 @@
 #include "Core/Basics/Complex.h"
 #include "Core/Correlations/FTDistributions2D.h"
 #include "Core/Lattice/Lattice2D.h"
-#include "Core/Tools/Integrator.h"
 #include <memory>
 
 class IFTDistribution2D;
@@ -91,7 +90,6 @@ private:
     double m_domain_sizes[2]; //!< Coherence domain sizes
     mutable double m_qx;
     mutable double m_qy;
-    mutable RealIntegrator m_integrator;
 };
 
 #endif // BORNAGAIN_CORE_AGGREGATE_INTERFERENCEFUNCTION2DPARACRYSTAL_H
diff --git a/Core/Aggregate/InterferenceFunction2DSuperLattice.cpp b/Core/Aggregate/InterferenceFunction2DSuperLattice.cpp
index 8a94239dc686236ac3dd910fd681bc47f651b365..b925b1753ad06cc17ffef299c41085137e96494b 100644
--- a/Core/Aggregate/InterferenceFunction2DSuperLattice.cpp
+++ b/Core/Aggregate/InterferenceFunction2DSuperLattice.cpp
@@ -17,6 +17,7 @@
 #include "Core/Basics/Exceptions.h"
 #include "Core/Basics/MathConstants.h"
 #include "Core/Parametrization/RealParameter.h"
+#include "Core/Tools/Integrator.h"
 #include "Core/Tools/MathFunctions.h"
 
 #include <limits>
@@ -95,8 +96,8 @@ double InterferenceFunction2DSuperLattice::evaluate(const kvector_t q, double ou
     m_qy = q.y();
     if (!m_integrate_xi)
         return interferenceForXi(mP_lattice->rotationAngle());
-    return m_integrator.integrate([&](double xi) -> double { return interferenceForXi(xi); }, 0.0,
-                                  M_TWOPI)
+    return RealIntegrator().integrate([&](double xi) -> double { return interferenceForXi(xi); },
+                                      0.0, M_TWOPI)
            / M_TWOPI;
 }
 
diff --git a/Core/Aggregate/InterferenceFunction2DSuperLattice.h b/Core/Aggregate/InterferenceFunction2DSuperLattice.h
index 85b4b102c672165007e8f6427056e124e98fa252..aadaceb843ecb22e370ebdef67cdb5d5ff175798 100644
--- a/Core/Aggregate/InterferenceFunction2DSuperLattice.h
+++ b/Core/Aggregate/InterferenceFunction2DSuperLattice.h
@@ -17,7 +17,6 @@
 
 #include "Core/Aggregate/IInterferenceFunction.h"
 #include "Core/Lattice/Lattice2D.h"
-#include "Core/Tools/Integrator.h"
 
 //! Interference function of a 2D superlattice with a configurable interference function for
 //! each lattice site.
@@ -68,7 +67,6 @@ private:
     mutable double m_qx;
     mutable double m_qy;
     mutable double m_xi;
-    mutable RealIntegrator m_integrator;
 };
 
 #endif // BORNAGAIN_CORE_AGGREGATE_INTERFERENCEFUNCTION2DSUPERLATTICE_H
diff --git a/Core/Aggregate/InterferenceFunctionFinite2DLattice.cpp b/Core/Aggregate/InterferenceFunctionFinite2DLattice.cpp
index 095b3e3c832b0971e31ca1a841d118b04dac9243..d68a86f7d7cf4456e5849858fd733b0cb269b7b2 100644
--- a/Core/Aggregate/InterferenceFunctionFinite2DLattice.cpp
+++ b/Core/Aggregate/InterferenceFunctionFinite2DLattice.cpp
@@ -16,6 +16,7 @@
 #include "Core/Basics/Exceptions.h"
 #include "Core/Basics/MathConstants.h"
 #include "Core/Parametrization/RealParameter.h"
+#include "Core/Tools/Integrator.h"
 #include "Core/Tools/MathFunctions.h"
 
 #include <limits>
@@ -114,8 +115,8 @@ double InterferenceFunctionFinite2DLattice::iff_without_dw(const kvector_t q) co
     m_qy = q.y();
     if (!m_integrate_xi)
         return interferenceForXi(mP_lattice->rotationAngle());
-    return m_integrator.integrate([&](double xi) -> double { return interferenceForXi(xi); }, 0.0,
-                                  M_TWOPI)
+    return RealIntegrator().integrate([&](double xi) -> double { return interferenceForXi(xi); },
+                                      0.0, M_TWOPI)
            / M_TWOPI;
 }
 
diff --git a/Core/Aggregate/InterferenceFunctionFinite2DLattice.h b/Core/Aggregate/InterferenceFunctionFinite2DLattice.h
index 2563802c67f63263bba7d41f0a3d0a6491b0788c..dcae005653b816202c29d158b3be07bb3b91b41b 100644
--- a/Core/Aggregate/InterferenceFunctionFinite2DLattice.h
+++ b/Core/Aggregate/InterferenceFunctionFinite2DLattice.h
@@ -17,7 +17,6 @@
 
 #include "Core/Aggregate/IInterferenceFunction.h"
 #include "Core/Lattice/Lattice2D.h"
-#include "Core/Tools/Integrator.h"
 
 //! Interference function of a finite 2D lattice.
 //! @ingroup interference
@@ -63,7 +62,6 @@ private:
     unsigned m_N_1, m_N_2; //!< Size of the finite lattice in lattice units
     mutable double m_qx;
     mutable double m_qy;
-    mutable RealIntegrator m_integrator;
 };
 
 #endif // BORNAGAIN_CORE_AGGREGATE_INTERFERENCEFUNCTIONFINITE2DLATTICE_H
diff --git a/Core/Aggregate/InterferenceFunctionHardDisk.cpp b/Core/Aggregate/InterferenceFunctionHardDisk.cpp
index fc1f769dcde0ab876eeb557032453d0ac95fc475..56948d53d45590707f63f7b486140400b31199d8 100644
--- a/Core/Aggregate/InterferenceFunctionHardDisk.cpp
+++ b/Core/Aggregate/InterferenceFunctionHardDisk.cpp
@@ -14,6 +14,7 @@
 
 #include "Core/Aggregate/InterferenceFunctionHardDisk.h"
 #include "Core/Parametrization/RealParameter.h"
+#include "Core/Tools/Integrator.h"
 #include "Core/Tools/MathFunctions.h"
 #include <cmath>
 
@@ -69,7 +70,7 @@ double InterferenceFunctionHardDisk::iff_without_dw(const kvector_t q) const
     m_s2 = S2(m_packing);
     double c_q =
         2.0 * M_PI
-        * m_integrator.integrate([&](double x) -> double { return integrand(x); }, 0.0, 1.0);
+        * RealIntegrator().integrate([&](double x) -> double { return integrand(x); }, 0.0, 1.0);
     double rho = 4.0 * m_packing / M_PI;
     return 1.0 / (1.0 - rho * c_q);
 }
diff --git a/Core/Aggregate/InterferenceFunctionHardDisk.h b/Core/Aggregate/InterferenceFunctionHardDisk.h
index 3845affa62eaa900c23695297a2b178e177aeda0..18c5b38dddc54ac822adbec1750ebb91237285f7 100644
--- a/Core/Aggregate/InterferenceFunctionHardDisk.h
+++ b/Core/Aggregate/InterferenceFunctionHardDisk.h
@@ -16,7 +16,6 @@
 #define BORNAGAIN_CORE_AGGREGATE_INTERFERENCEFUNCTIONHARDDISK_H
 
 #include "Core/Aggregate/IInterferenceFunction.h"
-#include "Core/Tools/Integrator.h"
 
 //! Percus-Yevick hard disk interference function.
 //!
@@ -49,7 +48,6 @@ private:
     double m_density;
 
     mutable double m_q, m_c_zero, m_s2, m_packing;
-    mutable RealIntegrator m_integrator;
 };
 
 #endif // BORNAGAIN_CORE_AGGREGATE_INTERFERENCEFUNCTIONHARDDISK_H
diff --git a/Core/Correlations/FTDistributions2D.cpp b/Core/Correlations/FTDistributions2D.cpp
index 6c7d334f522a4010b1a3440b276c84c9bd73a53c..8f0b8d942c2e75501cdbd1c7103df766858c9649 100644
--- a/Core/Correlations/FTDistributions2D.cpp
+++ b/Core/Correlations/FTDistributions2D.cpp
@@ -16,6 +16,7 @@
 #include "Core/Basics/Algorithms.h"
 #include "Core/Basics/Exceptions.h"
 #include "Core/Basics/MathConstants.h"
+#include "Core/Tools/Integrator.h"
 #include "Core/Tools/MathFunctions.h"
 #include <limits>
 
@@ -155,7 +156,7 @@ double FTDistribution2DCone::evaluate(double qx, double qy) const
     if (scaled_q < std::numeric_limits<double>::epsilon())
         return 1.0 - 3.0 * scaled_q * scaled_q / 40.0;
     // second part of the integrand: \f$u^2\cdot J_0(u)\f$
-    double integral = m_integrator.integrate(
+    double integral = RealIntegrator().integrate(
         [](double x) -> double { return x * x * MathFunctions::Bessel_J0(x); }, 0.0, scaled_q);
     return 6.0 * (MathFunctions::Bessel_J1c(scaled_q) - integral / scaled_q / scaled_q / scaled_q);
 }
diff --git a/Core/Correlations/FTDistributions2D.h b/Core/Correlations/FTDistributions2D.h
index b9af89f7d0f49405cc48c23e8724134bd576e929..3e782e12f42fcae44e8096c9009e2ab33bfa1b78 100644
--- a/Core/Correlations/FTDistributions2D.h
+++ b/Core/Correlations/FTDistributions2D.h
@@ -19,7 +19,6 @@
 #include "Core/Basics/MathConstants.h"
 #include "Core/Correlations/IDistribution2DSampler.h"
 #include "Core/Parametrization/INode.h"
-#include "Core/Tools/Integrator.h"
 
 //! Interface for two-dimensional distributions in Fourier space.
 //! @ingroup distribution_internal
@@ -129,9 +128,6 @@ public:
 #ifndef SWIG
     std::unique_ptr<IDistribution2DSampler> createSampler() const final;
 #endif
-
-private:
-    mutable RealIntegrator m_integrator;
 };
 
 //! Two-dimensional Voigt distribution in Fourier space;
diff --git a/Core/Correlations/IPeakShape.cpp b/Core/Correlations/IPeakShape.cpp
index 2104274ebd7794a5fc3992f7061326d0655b2bff..af0d83456a2060e3b0cc08581dfb25dddd4e9d72 100644
--- a/Core/Correlations/IPeakShape.cpp
+++ b/Core/Correlations/IPeakShape.cpp
@@ -14,8 +14,8 @@
 
 #include "Core/Correlations/IPeakShape.h"
 #include "Core/Basics/MathConstants.h"
+#include "Core/Tools/Integrator.h"
 #include "Core/Tools/MathFunctions.h"
-
 #include <limits>
 
 namespace
@@ -252,8 +252,8 @@ double MisesFisherGaussPeakShape::evaluate(const kvector_t q, const kvector_t q_
     m_theta = std::acos(q.unit().dot(m_zenith));
     double pre_1 = FisherPrefactor(m_kappa_1);
     double pre_2 = MisesPrefactor(m_kappa_2);
-    double integral =
-        m_integrator.integrate([&](double phi) -> double { return integrand(phi); }, 0.0, M_TWOPI);
+    double integral = RealIntegrator().integrate(
+        [&](double phi) -> double { return integrand(phi); }, 0.0, M_TWOPI);
     return m_max_intensity * radial_part * pre_1 * pre_2 * integral;
 }
 
@@ -300,8 +300,8 @@ double MisesGaussPeakShape::evaluate(const kvector_t q, const kvector_t q_lattic
     m_phi = std::acos(q_ortho.unit().dot(m_ux));
     m_theta = std::acos(q.unit().dot(m_zenith));
     double pre = MisesPrefactor(m_kappa);
-    double integral =
-        m_integrator.integrate([&](double phi) -> double { return integrand(phi); }, 0.0, M_TWOPI);
+    double integral = RealIntegrator().integrate(
+        [&](double phi) -> double { return integrand(phi); }, 0.0, M_TWOPI);
     return m_max_intensity * pre * integral;
 }
 
diff --git a/Core/Correlations/IPeakShape.h b/Core/Correlations/IPeakShape.h
index 0067487823f85b75ea90ab4b63c642ec84595acc..1f7eb7f5c4352b83ed38812f7666ee0120ea4e9e 100644
--- a/Core/Correlations/IPeakShape.h
+++ b/Core/Correlations/IPeakShape.h
@@ -16,7 +16,6 @@
 #define BORNAGAIN_CORE_CORRELATIONS_IPEAKSHAPE_H
 
 #include "Core/Scattering/ISample.h"
-#include "Core/Tools/Integrator.h"
 #include "Core/Vector/Vectors3D.h"
 
 //! Pure virtual interface class that defines the peak shape of a Bragg peak.
@@ -164,7 +163,6 @@ private:
     double m_kappa_1, m_kappa_2;
     mutable double m_theta, m_phi;
     mutable kvector_t m_ux, m_uy, m_up;
-    mutable RealIntegrator m_integrator;
 };
 
 //! Class that implements a peak shape that is a convolution of a Mises-Fisher distribution
@@ -194,7 +192,6 @@ private:
     double m_kappa;
     mutable double m_theta, m_phi, m_qr;
     mutable kvector_t m_ux, m_uy, m_p;
-    mutable RealIntegrator m_integrator;
 };
 
 #endif // BORNAGAIN_CORE_CORRELATIONS_IPEAKSHAPE_H
diff --git a/Core/HardParticle/FormFactorCone.cpp b/Core/HardParticle/FormFactorCone.cpp
index 3959d13a92490ab0788c0a018eace9412e3139bc..77474e7a463a1197233822159870b0d7cb7aa71f 100644
--- a/Core/HardParticle/FormFactorCone.cpp
+++ b/Core/HardParticle/FormFactorCone.cpp
@@ -16,6 +16,7 @@
 #include "Core/Basics/Exceptions.h"
 #include "Core/Basics/MathConstants.h"
 #include "Core/Shapes/DoubleEllipse.h"
+#include "Core/Tools/Integrator.h"
 #include "Core/Tools/MathFunctions.h"
 #include <limits>
 
@@ -73,7 +74,7 @@ complex_t FormFactorCone::evaluate_for_q(cvector_t q) const
         return M_PI / 3. * (R * R * H + (R * R - R2 * R2) * (apex_height - H));
     } else {
         complex_t integral =
-            m_integrator.integrate([&](double Z) { return Integrand(Z); }, 0., m_height);
+            ComplexIntegrator().integrate([&](double Z) { return Integrand(Z); }, 0., m_height);
         return M_TWOPI * integral;
     }
 }
diff --git a/Core/HardParticle/FormFactorCone.h b/Core/HardParticle/FormFactorCone.h
index 73913a44030fafeda81130276c8814823dba70bb..ff8b31ada3d2060c0e5caf6553a38d5b625df64c 100644
--- a/Core/HardParticle/FormFactorCone.h
+++ b/Core/HardParticle/FormFactorCone.h
@@ -16,7 +16,6 @@
 #define BORNAGAIN_CORE_HARDPARTICLE_FORMFACTORCONE_H
 
 #include "Core/Scattering/IFormFactorBorn.h"
-#include "Core/Tools/Integrator.h"
 
 //! A conical frustum (cone truncated parallel to the base) with circular base.
 //! @ingroup hardParticle
@@ -55,7 +54,6 @@ private:
     const double& m_alpha;
     double m_cot_alpha;
     mutable cvector_t m_q;
-    mutable ComplexIntegrator m_integrator;
 };
 
 #endif // BORNAGAIN_CORE_HARDPARTICLE_FORMFACTORCONE_H
diff --git a/Core/HardParticle/FormFactorHemiEllipsoid.cpp b/Core/HardParticle/FormFactorHemiEllipsoid.cpp
index 49875b5501908392aed21b51af8479828c6d2311..ee73ec817ff6627ba6cd075118aa202216fca58f 100644
--- a/Core/HardParticle/FormFactorHemiEllipsoid.cpp
+++ b/Core/HardParticle/FormFactorHemiEllipsoid.cpp
@@ -15,6 +15,7 @@
 #include "Core/HardParticle/FormFactorHemiEllipsoid.h"
 #include "Core/Basics/MathConstants.h"
 #include "Core/Shapes/TruncatedEllipsoid.h"
+#include "Core/Tools/Integrator.h"
 #include "Core/Tools/MathFunctions.h"
 #include <limits>
 
@@ -72,7 +73,7 @@ complex_t FormFactorHemiEllipsoid::evaluate_for_q(cvector_t q) const
 
     if (std::abs(m_q.mag()) <= std::numeric_limits<double>::epsilon())
         return M_TWOPI * R * W * H / 3.;
-    return M_TWOPI * m_integrator.integrate([&](double Z) { return Integrand(Z); }, 0., H);
+    return M_TWOPI * ComplexIntegrator().integrate([&](double Z) { return Integrand(Z); }, 0., H);
 }
 
 void FormFactorHemiEllipsoid::onChange()
diff --git a/Core/HardParticle/FormFactorHemiEllipsoid.h b/Core/HardParticle/FormFactorHemiEllipsoid.h
index e3fdaef86a8ccec9049d6d3d6f2d5174d55df1ef..7453c429dea20fb1cf262cadd74c823b3d6a9ebb 100644
--- a/Core/HardParticle/FormFactorHemiEllipsoid.h
+++ b/Core/HardParticle/FormFactorHemiEllipsoid.h
@@ -16,7 +16,6 @@
 #define BORNAGAIN_CORE_HARDPARTICLE_FORMFACTORHEMIELLIPSOID_H
 
 #include "Core/Scattering/IFormFactorBorn.h"
-#include "Core/Tools/Integrator.h"
 
 //! An hemi ellipsoid,
 //!   obtained by truncating a full ellipsoid in the middle plane spanned by two principal axes.
@@ -53,7 +52,6 @@ private:
     const double& m_radius_y;
     const double& m_height;
     mutable cvector_t m_q;
-    mutable ComplexIntegrator m_integrator;
 };
 
 #endif // BORNAGAIN_CORE_HARDPARTICLE_FORMFACTORHEMIELLIPSOID_H
diff --git a/Core/HardParticle/FormFactorTruncatedSphere.cpp b/Core/HardParticle/FormFactorTruncatedSphere.cpp
index a09bd3b4338bb6a8788c1d8c9970fe05c47aed2d..0bbec494e3501aa2a4e6a248d0aaa7612f7cdac1 100644
--- a/Core/HardParticle/FormFactorTruncatedSphere.cpp
+++ b/Core/HardParticle/FormFactorTruncatedSphere.cpp
@@ -16,6 +16,7 @@
 #include "Core/Basics/Exceptions.h"
 #include "Core/Basics/MathConstants.h"
 #include "Core/Shapes/TruncatedEllipsoid.h"
+#include "Core/Tools/Integrator.h"
 #include "Core/Tools/MathFunctions.h"
 #include "Fit/Tools/RealLimits.h"
 #include <limits>
@@ -76,8 +77,8 @@ complex_t FormFactorTruncatedSphere::evaluate_for_q(cvector_t q) const
                   - m_dh * m_dh * (3. * m_radius - m_dh));
     }
     // else
-    complex_t integral = m_integrator.integrate([&](double Z) { return Integrand(Z); },
-                                                m_radius - m_height, m_radius - m_dh);
+    complex_t integral = ComplexIntegrator().integrate([&](double Z) { return Integrand(Z); },
+                                                       m_radius - m_height, m_radius - m_dh);
     return M_TWOPI * integral * exp_I(q.z() * (m_height - m_radius));
 }
 
diff --git a/Core/HardParticle/FormFactorTruncatedSphere.h b/Core/HardParticle/FormFactorTruncatedSphere.h
index a646b5093474e908499c6447c0b0f5e2100dcf22..aaab6b54cd46e7c74b56be9f32eda6344adacef2 100644
--- a/Core/HardParticle/FormFactorTruncatedSphere.h
+++ b/Core/HardParticle/FormFactorTruncatedSphere.h
@@ -16,7 +16,6 @@
 #define BORNAGAIN_CORE_HARDPARTICLE_FORMFACTORTRUNCATEDSPHERE_H
 
 #include "Core/Scattering/IFormFactorBorn.h"
-#include "Core/Tools/Integrator.h"
 
 //! A truncated Sphere.
 //! @ingroup hardParticle
@@ -55,7 +54,6 @@ private:
     const double& m_height;
     const double& m_dh;
     mutable cvector_t m_q;
-    mutable ComplexIntegrator m_integrator;
 };
 
 #endif // BORNAGAIN_CORE_HARDPARTICLE_FORMFACTORTRUNCATEDSPHERE_H
diff --git a/Core/HardParticle/FormFactorTruncatedSpheroid.cpp b/Core/HardParticle/FormFactorTruncatedSpheroid.cpp
index 561b14d73cbe969407da37cb4c92e11f013a1e52..b7e93ff209271203c2b1c63d88ee074a293e20dc 100644
--- a/Core/HardParticle/FormFactorTruncatedSpheroid.cpp
+++ b/Core/HardParticle/FormFactorTruncatedSpheroid.cpp
@@ -16,6 +16,7 @@
 #include "Core/Basics/Exceptions.h"
 #include "Core/Basics/MathConstants.h"
 #include "Core/Shapes/TruncatedEllipsoid.h"
+#include "Core/Tools/Integrator.h"
 #include "Core/Tools/MathFunctions.h"
 #include <limits>
 
@@ -83,8 +84,8 @@ complex_t FormFactorTruncatedSpheroid::evaluate_for_q(cvector_t q) const
         return M_PI / 3. / fp * (H * H * (3. * R - H / fp) - m_dh * m_dh * (3. * R - m_dh / fp));
     complex_t z_part = std::exp(complex_t(0.0, 1.0) * m_q.z() * (H - fp * R));
     return M_TWOPI * z_part
-           * m_integrator.integrate([&](double Z) { return Integrand(Z); }, fp * R - H,
-                                    fp * R - m_dh);
+           * ComplexIntegrator().integrate([&](double Z) { return Integrand(Z); }, fp * R - H,
+                                           fp * R - m_dh);
 }
 
 IFormFactor* FormFactorTruncatedSpheroid::sliceFormFactor(ZLimits limits, const IRotation& rot,
diff --git a/Core/HardParticle/FormFactorTruncatedSpheroid.h b/Core/HardParticle/FormFactorTruncatedSpheroid.h
index 7b3e449efe0d5a75f4d1bdb4935cdb345c3076f9..887fa7a8ae2b2dcf62ad44f9db191ec3fb66a338 100644
--- a/Core/HardParticle/FormFactorTruncatedSpheroid.h
+++ b/Core/HardParticle/FormFactorTruncatedSpheroid.h
@@ -16,7 +16,6 @@
 #define BORNAGAIN_CORE_HARDPARTICLE_FORMFACTORTRUNCATEDSPHEROID_H
 
 #include "Core/Scattering/IFormFactorBorn.h"
-#include "Core/Tools/Integrator.h"
 
 //! A truncated spheroid.
 //! An ellipsoid with two equal axis, truncated by a plane perpendicular to the third axis.
@@ -58,7 +57,6 @@ private:
     const double& m_height_flattening;
     const double& m_dh;
     mutable cvector_t m_q;
-    mutable ComplexIntegrator m_integrator;
 };
 
 #endif // BORNAGAIN_CORE_HARDPARTICLE_FORMFACTORTRUNCATEDSPHEROID_H
diff --git a/Core/Intensity/SimpleUnitConverters.cpp b/Core/Intensity/SimpleUnitConverters.cpp
index 85de14fef2ad77f9b0c3df04ac995a55db1b37b7..e3238c2d5f0f71d27fb0ccb78af8f9520a3f2620 100644
--- a/Core/Intensity/SimpleUnitConverters.cpp
+++ b/Core/Intensity/SimpleUnitConverters.cpp
@@ -17,9 +17,9 @@
 #include "Core/Basics/Units.h"
 #include "Core/Beam/Beam.h"
 #include "Core/Detector/RectangularDetector.h"
+#include "Core/Detector/RectangularPixel.h"
 #include "Core/Detector/RegionOfInterest.h"
 #include "Core/Detector/SphericalDetector.h"
-#include "Core/Detector/RectangularPixel.h"
 #include "Core/Intensity/AxisNames.h"
 #include <algorithm>
 #include <cmath>
diff --git a/auto/Wrap/libBornAgainCore.py b/auto/Wrap/libBornAgainCore.py
index cd73caeeba071c64e077763666f5980c07078a6d..0cf0c748d5d6fbc77bfc1c7f860a4820cafb6295 100644
--- a/auto/Wrap/libBornAgainCore.py
+++ b/auto/Wrap/libBornAgainCore.py
@@ -12116,7 +12116,6 @@ class IInterferenceFunction(ISample):
         else:
             _self = self
         _libBornAgainCore.IInterferenceFunction_swiginit(self, _libBornAgainCore.new_IInterferenceFunction(_self, *args))
-    __swig_destroy__ = _libBornAgainCore.delete_IInterferenceFunction
 
     def clone(self):
         r"""
@@ -12191,6 +12190,7 @@ class IInterferenceFunction(ISample):
     def iff_without_dw(self, q):
         r"""iff_without_dw(IInterferenceFunction self, kvector_t q) -> double"""
         return _libBornAgainCore.IInterferenceFunction_iff_without_dw(self, q)
+    __swig_destroy__ = _libBornAgainCore.delete_IInterferenceFunction
     def __disown__(self):
         self.this.disown()
         _libBornAgainCore.disown_IInterferenceFunction(self)
diff --git a/auto/Wrap/libBornAgainCore_wrap.cpp b/auto/Wrap/libBornAgainCore_wrap.cpp
index ee57d5406db69109fbe2b5605fe91903dfadc3c9..c65e2a55c3d8bd95f25c9bb03c96ec4ce9ab91b5 100644
--- a/auto/Wrap/libBornAgainCore_wrap.cpp
+++ b/auto/Wrap/libBornAgainCore_wrap.cpp
@@ -82450,28 +82450,6 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_delete_IInterferenceFunction(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
-  PyObject *resultobj = 0;
-  IInterferenceFunction *arg1 = (IInterferenceFunction *) 0 ;
-  void *argp1 = 0 ;
-  int res1 = 0 ;
-  PyObject *swig_obj[1] ;
-  
-  if (!args) SWIG_fail;
-  swig_obj[0] = args;
-  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_IInterferenceFunction, SWIG_POINTER_DISOWN |  0 );
-  if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "delete_IInterferenceFunction" "', argument " "1"" of type '" "IInterferenceFunction *""'"); 
-  }
-  arg1 = reinterpret_cast< IInterferenceFunction * >(argp1);
-  delete arg1;
-  resultobj = SWIG_Py_Void();
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
 SWIGINTERN PyObject *_wrap_IInterferenceFunction_clone(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
   IInterferenceFunction *arg1 = (IInterferenceFunction *) 0 ;
@@ -82888,6 +82866,28 @@ fail:
 }
 
 
+SWIGINTERN PyObject *_wrap_delete_IInterferenceFunction(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+  PyObject *resultobj = 0;
+  IInterferenceFunction *arg1 = (IInterferenceFunction *) 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  PyObject *swig_obj[1] ;
+  
+  if (!args) SWIG_fail;
+  swig_obj[0] = args;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_IInterferenceFunction, SWIG_POINTER_DISOWN |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "delete_IInterferenceFunction" "', argument " "1"" of type '" "IInterferenceFunction *""'"); 
+  }
+  arg1 = reinterpret_cast< IInterferenceFunction * >(argp1);
+  delete arg1;
+  resultobj = SWIG_Py_Void();
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
 SWIGINTERN PyObject *_wrap_disown_IInterferenceFunction(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
   IInterferenceFunction *arg1 = (IInterferenceFunction *) 0 ;
@@ -128385,11 +128385,6 @@ static PyMethodDef SwigMethods[] = {
 		"IInterferenceFunction::IInterferenceFunction(const IInterferenceFunction &other)\n"
 		"\n"
 		""},
-	 { "delete_IInterferenceFunction", _wrap_delete_IInterferenceFunction, METH_O, "\n"
-		"delete_IInterferenceFunction(IInterferenceFunction self)\n"
-		"IInterferenceFunction::~IInterferenceFunction()\n"
-		"\n"
-		""},
 	 { "IInterferenceFunction_clone", _wrap_IInterferenceFunction_clone, METH_O, "\n"
 		"IInterferenceFunction_clone(IInterferenceFunction self) -> IInterferenceFunction\n"
 		"virtual IInterferenceFunction* IInterferenceFunction::clone() const =0\n"
@@ -128440,6 +128435,11 @@ static PyMethodDef SwigMethods[] = {
 		"\n"
 		""},
 	 { "IInterferenceFunction_iff_without_dw", _wrap_IInterferenceFunction_iff_without_dw, METH_VARARGS, "IInterferenceFunction_iff_without_dw(IInterferenceFunction self, kvector_t q) -> double"},
+	 { "delete_IInterferenceFunction", _wrap_delete_IInterferenceFunction, METH_O, "\n"
+		"delete_IInterferenceFunction(IInterferenceFunction self)\n"
+		"IInterferenceFunction::~IInterferenceFunction()\n"
+		"\n"
+		""},
 	 { "disown_IInterferenceFunction", _wrap_disown_IInterferenceFunction, METH_O, NULL},
 	 { "IInterferenceFunction_swigregister", IInterferenceFunction_swigregister, METH_O, NULL},
 	 { "IInterferenceFunction_swiginit", IInterferenceFunction_swiginit, METH_VARARGS, NULL},