diff --git a/Base/Types/Exceptions.cpp b/Base/Types/Exceptions.cpp
index 831f3f8b7b5a060480b9221720736ea1bb643579..03b2ea6b99a85168a5b335d480e0ff5f39a3f3aa 100644
--- a/Base/Types/Exceptions.cpp
+++ b/Base/Types/Exceptions.cpp
@@ -45,18 +45,6 @@ ClassInitializationException::ClassInitializationException(const std::string& me
     LogExceptionMessage(message);
 }
 
-UnknownClassRegistrationException::UnknownClassRegistrationException(const std::string& message)
-    : std::runtime_error(message)
-{
-    LogExceptionMessage(message);
-}
-
-ExistingClassRegistrationException::ExistingClassRegistrationException(const std::string& message)
-    : std::runtime_error(message)
-{
-    LogExceptionMessage(message);
-}
-
 LogicErrorException::LogicErrorException(const std::string& message) : std::logic_error(message)
 {
     LogExceptionMessage(message);
@@ -68,12 +56,6 @@ RuntimeErrorException::RuntimeErrorException(const std::string& message)
     LogExceptionMessage(message);
 }
 
-DivisionByZeroException::DivisionByZeroException(const std::string& message)
-    : std::runtime_error(message)
-{
-    LogExceptionMessage(message);
-}
-
 DomainErrorException::DomainErrorException(const std::string& message) : std::domain_error(message)
 {
     LogExceptionMessage(message);
diff --git a/Base/Types/Exceptions.h b/Base/Types/Exceptions.h
index b9887b74a64b65a48883829c2c403a1683e77ce8..db71f9c1b0463b1d8d23fedf669279c19de74864 100644
--- a/Base/Types/Exceptions.h
+++ b/Base/Types/Exceptions.h
@@ -54,18 +54,6 @@ public:
     ClassInitializationException(const std::string& message);
 };
 
-class UnknownClassRegistrationException : public std::runtime_error
-{
-public:
-    UnknownClassRegistrationException(const std::string& message);
-};
-
-class ExistingClassRegistrationException : public std::runtime_error
-{
-public:
-    ExistingClassRegistrationException(const std::string& message);
-};
-
 class LogicErrorException : public std::logic_error
 {
 public:
@@ -78,12 +66,6 @@ public:
     RuntimeErrorException(const std::string& message);
 };
 
-class DivisionByZeroException : public std::runtime_error
-{
-public:
-    DivisionByZeroException(const std::string& message);
-};
-
 class DomainErrorException : public std::domain_error
 {
 public:
diff --git a/Base/Utils/Bessel.cpp b/Base/Utils/Bessel.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..3dfdf408d293b35fee9894473c0113ed13af0245
--- /dev/null
+++ b/Base/Utils/Bessel.cpp
@@ -0,0 +1,204 @@
+// ************************************************************************** //
+//
+//  BornAgain: simulate and fit scattering at grazing incidence
+//
+//! @file      Base/Utils/Bessel.cpp
+//! @brief     Implements Bessel functions in namespace Math.
+//!
+//! @homepage  http://www.bornagainproject.org
+//! @license   GNU General Public License v3 or higher (see COPYING)
+//! @copyright Forschungszentrum Jülich GmbH 2018
+//! @authors   Scientific Computing Group at MLZ (see CITATION, AUTHORS)
+//
+// ************************************************************************** //
+
+#include "Base/Utils/Bessel.h"
+#include "Base/Const/MathConstants.h"
+#include <gsl/gsl_sf_bessel.h>
+
+namespace
+{
+
+//! Computes the complex Bessel function J0(z), using power series and asymptotic expansion.
+//!
+//! Forked from unoptimized code at http://www.crbond.com/math.htm,
+//! who refers to "Computation of Special Functions", Zhang and Jin, John Wiley and Sons, 1996.
+
+complex_t J0_PowSer(const complex_t z)
+{
+    complex_t cj0;
+    static const double eps = 1e-15;
+    static double a[] = {-7.03125e-2,           0.112152099609375,     -0.5725014209747314,
+                         6.074042001273483,     -1.100171402692467e2,  3.038090510922384e3,
+                         -1.188384262567832e5,  6.252951493434797e6,   -4.259392165047669e8,
+                         3.646840080706556e10,  -3.833534661393944e12, 4.854014686852901e14,
+                         -7.286857349377656e16, 1.279721941975975e19};
+    static double b[] = {7.32421875e-2,         -0.2271080017089844,  1.727727502584457,
+                         -2.438052969955606e1,  5.513358961220206e2,  -1.825775547429318e4,
+                         8.328593040162893e5,   -5.006958953198893e7, 3.836255180230433e9,
+                         -3.649010818849833e11, 4.218971570284096e13, -5.827244631566907e15,
+                         9.476288099260110e17,  -1.792162323051699e20};
+
+    double a0 = std::abs(z);
+    complex_t z1 = z;
+    if (a0 == 0.0)
+        return 1.0;
+    if (std::real(z) < 0.0)
+        z1 = -z;
+    if (a0 <= 12.0) {
+        // standard power series [http://dlmf.nist.gov/10.2 (10.2.2)]
+        complex_t z2 = 0.25 * z * z;
+        cj0 = 1.0;
+        complex_t cr = 1.0;
+        for (size_t k = 1; k <= 40; ++k) {
+            cr *= -z2 / (double)(k * k);
+            cj0 += cr;
+            if (std::abs(cr) < std::abs(cj0) * eps)
+                break;
+        }
+    } else {
+        // Hankel's asymptotic expansion [http://dlmf.nist.gov/10.17 (10.17.3)]
+        size_t kz;
+        if (a0 >= 50.0)
+            kz = 8; // can be changed to 10
+        else if (a0 >= 35.0)
+            kz = 10; //   "      "     "  12
+        else
+            kz = 12; //   "      "     "  14
+        complex_t ct1 = z1 - M_PI_4;
+        complex_t cp0 = 1.0;
+        complex_t cq0 = -0.125;
+        const complex_t z1m2 = 1. / (z1 * z1); // faster than std::pow(z1, -2.0) ??
+        complex_t ptmp = z1m2;
+        for (size_t k = 0; k < kz; ++k) {
+            cp0 += a[k] * ptmp;
+            cq0 += b[k] * ptmp;
+            ptmp *= z1m2;
+        }
+        cj0 = std::sqrt(M_2_PI / z1) * (cp0 * std::cos(ct1) - cq0 / z1 * std::sin(ct1));
+    }
+    return cj0;
+}
+
+//! Computes the complex Bessel function J1(z), using power series and asymptotic expansion.
+//!
+//! Forked from same source as for J0_PowSer
+
+complex_t J1_PowSer(const complex_t z)
+{
+    complex_t cj1;
+    static const double eps = 1e-15;
+
+    static double a1[] = {0.1171875,
+                          -0.1441955566406250,
+                          0.6765925884246826,
+                          -6.883914268109947,
+                          1.215978918765359e2,
+                          -3.302272294480852e3,
+                          1.276412726461746e5,
+                          -6.656367718817688e6,
+                          4.502786003050393e8,
+                          -3.833857520742790e10,
+                          4.011838599133198e12,
+                          -5.060568503314727e14,
+                          7.572616461117958e16,
+                          -1.326257285320556e19};
+    static double b1[] = {-0.1025390625,         0.2775764465332031,    -1.993531733751297,
+                          2.724882731126854e1,   -6.038440767050702e2,  1.971837591223663e4,
+                          -8.902978767070678e5,  5.310411010968522e7,   -4.043620325107754e9,
+                          3.827011346598605e11,  -4.406481417852278e13, 6.065091351222699e15,
+                          -9.833883876590679e17, 1.855045211579828e20};
+
+    double a0 = std::abs(z);
+    if (a0 == 0.0)
+        return 0.0;
+
+    complex_t z1 = z;
+    if (std::real(z) < 0.0)
+        z1 = -z;
+    if (a0 <= 12.0) {
+        // standard power series [http://dlmf.nist.gov/10.2 (10.2.2)]
+        const complex_t z2 = 0.25 * z * z;
+        cj1 = 1.0;
+        complex_t cr = 1.0; // powers will be computed recursively
+        for (int k = 1; k <= 40; ++k) {
+            cr *= -z2 / (double)(k * (k + 1));
+            cj1 += cr;
+            if (std::abs(cr) < std::abs(cj1) * eps)
+                break;
+        }
+        cj1 *= 0.5 * z1;
+    } else {
+        // Hankel's asymptotic expansion [http://dlmf.nist.gov/10.17 (10.17.3)]
+        size_t kz;
+        if (a0 >= 50.0)
+            kz = 8; // can be changed to 10
+        else if (a0 >= 35.0)
+            kz = 10; //   "      "     "  12
+        else
+            kz = 12; //   "      "     "  14
+        complex_t cp1 = 1.0;
+        complex_t cq1 = 0.375;                 // division by z1 postponed to final sum
+        const complex_t z1m2 = 1. / (z1 * z1); // faster than std::pow(z1, -2.0) ??
+        complex_t ptmp = z1m2;                 // powers will be computed recursively
+        for (size_t k = 0; k < kz; ++k) {
+            cp1 += a1[k] * ptmp;
+            cq1 += b1[k] * ptmp; // division by z1 postponed to final sum
+            ptmp *= z1m2;
+        }
+        const complex_t ct2 = z1 - 0.75 * M_PI;
+        cj1 = std::sqrt(M_2_PI / z1) * (cp1 * std::cos(ct2) - cq1 / z1 * std::sin(ct2));
+    }
+    if (std::real(z) < 0.0)
+        cj1 = -cj1;
+    return cj1;
+}
+
+} // namespace
+
+// ************************************************************************** //
+//  Bessel functions
+// ************************************************************************** //
+
+double Math::Bessel::J0(double x)
+{
+    return gsl_sf_bessel_J0(x);
+}
+
+double Math::Bessel::J1(double x)
+{
+    return gsl_sf_bessel_J1(x);
+}
+
+double Math::Bessel::J1c(double x)
+{
+    return x == 0 ? 0.5 : gsl_sf_bessel_J1(x) / x;
+}
+
+double Math::Bessel::I0(double x)
+{
+    return gsl_sf_bessel_I0(x);
+}
+
+complex_t Math::Bessel::J0(const complex_t z)
+{
+    if (std::imag(z) == 0)
+        return gsl_sf_bessel_J0(std::real(z));
+    return J0_PowSer(z);
+}
+
+complex_t Math::Bessel::J1(const complex_t z)
+{
+    if (std::imag(z) == 0)
+        return gsl_sf_bessel_J1(std::real(z));
+    return J1_PowSer(z);
+}
+
+complex_t Math::Bessel::J1c(const complex_t z)
+{
+    if (std::imag(z) == 0) {
+        double xv = std::real(z);
+        return xv == 0 ? 0.5 : gsl_sf_bessel_J1(xv) / xv;
+    }
+    return z == 0. ? 0.5 : J1_PowSer(z) / z;
+}
diff --git a/Base/Utils/Bessel.h b/Base/Utils/Bessel.h
new file mode 100644
index 0000000000000000000000000000000000000000..ddd7ba25a2eafb760df26d85a27c8ac6688f18a5
--- /dev/null
+++ b/Base/Utils/Bessel.h
@@ -0,0 +1,50 @@
+// ************************************************************************** //
+//
+//  BornAgain: simulate and fit scattering at grazing incidence
+//
+//! @file      Base/Utils/Bessel.h
+//! @brief     Defines Bessel functions in namespace Math.
+//!
+//! @homepage  http://www.bornagainproject.org
+//! @license   GNU General Public License v3 or higher (see COPYING)
+//! @copyright Forschungszentrum Jülich GmbH 2018
+//! @authors   Scientific Computing Group at MLZ (see CITATION, AUTHORS)
+//
+// ************************************************************************** //
+
+#ifndef BORNAGAIN_BASE_UTILS_BESSEL_H
+#define BORNAGAIN_BASE_UTILS_BESSEL_H
+
+#include "Base/Types/Complex.h"
+#include <vector>
+
+namespace Math
+{
+namespace Bessel
+{
+
+//! Bessel function of the first kind and order 0
+double J0(double x);
+
+//! Bessel function of the first kind and order 1
+double J1(double x);
+
+//! Bessel function  J1(x)/x
+double J1c(double x);
+
+//! Modified Bessel function of the first kind and order 0
+double I0(double x);
+
+//! Complex Bessel function of the first kind and order 0
+complex_t J0(const complex_t z);
+
+//! Complex Bessel function of the first kind and order 1
+complex_t J1(const complex_t z);
+
+//! Complex Bessel function  J1(x)/x
+complex_t J1c(const complex_t z);
+
+} // namespace Bessel
+} // namespace Math
+
+#endif // BORNAGAIN_BASE_UTILS_BESSEL_H
diff --git a/Base/Utils/MathFunctions.cpp b/Base/Utils/MathFunctions.cpp
index f01f622f305c04302ebeaadbab600ac09f589059..3b497aeacdfff7341c88a28f61e6299ef3e61456 100644
--- a/Base/Utils/MathFunctions.cpp
+++ b/Base/Utils/MathFunctions.cpp
@@ -3,7 +3,7 @@
 //  BornAgain: simulate and fit scattering at grazing incidence
 //
 //! @file      Base/Utils/MathFunctions.cpp
-//! @brief     Implements namespace MathFunctions.
+//! @brief     Implements functions in namespace Math.
 //!
 //! @homepage  http://www.bornagainproject.org
 //! @license   GNU General Public License v3 or higher (see COPYING)
@@ -12,12 +12,11 @@
 //
 // ************************************************************************** //
 
-#include "Base/Utils/MathFunctions.h"
 #include "Base/Const/MathConstants.h"
+#include "Base/Utils/MathFunctions.h"
 #include <chrono>
 #include <cstring>
 #include <fftw3.h>
-#include <gsl/gsl_sf_bessel.h>
 #include <gsl/gsl_sf_erf.h>
 #include <gsl/gsl_sf_expint.h>
 #include <gsl/gsl_sf_trig.h>
@@ -29,39 +28,34 @@
 //  Various functions
 // ************************************************************************** //
 
-double MathFunctions::StandardNormal(double x)
+double Math::StandardNormal(double x)
 {
     return std::exp(-x * x / 2.0) / std::sqrt(M_TWOPI);
 }
 
-double MathFunctions::Gaussian(double x, double average, double std_dev)
+double Math::Gaussian(double x, double average, double std_dev)
 {
     return StandardNormal((x - average) / std_dev) / std_dev;
 }
 
-double MathFunctions::IntegratedGaussian(double x, double average, double std_dev)
+double Math::IntegratedGaussian(double x, double average, double std_dev)
 {
     double normalized_x = (x - average) / std_dev;
     static double root2 = std::sqrt(2.0);
     return (gsl_sf_erf(normalized_x / root2) + 1.0) / 2.0;
 }
 
-double MathFunctions::cot(double x)
+double Math::cot(double x)
 {
     return tan(M_PI_2 - x);
 }
 
-double MathFunctions::Si(double x) // int_0^x du Sin(u)/u
-{
-    return gsl_sf_Si(x);
-}
-
-double MathFunctions::sinc(double x) // Sin(x)/x
+double Math::sinc(double x) // Sin(x)/x
 {
     return gsl_sf_sinc(x / M_PI);
 }
 
-complex_t MathFunctions::sinc(const complex_t z) // Sin(x)/x
+complex_t Math::sinc(const complex_t z) // Sin(x)/x
 {
     // This is an exception from the rule that we must not test floating-point numbers for equality.
     // For small non-zero arguments, sin(z) returns quite accurately z or z-z^3/6.
@@ -72,14 +66,14 @@ complex_t MathFunctions::sinc(const complex_t z) // Sin(x)/x
     return std::sin(z) / z;
 }
 
-complex_t MathFunctions::tanhc(const complex_t z) // tanh(x)/x
+complex_t Math::tanhc(const complex_t z) // tanh(x)/x
 {
     if (std::abs(z) < std::numeric_limits<double>::epsilon())
         return 1.0;
     return std::tanh(z) / z;
 }
 
-double MathFunctions::Laue(const double x, size_t N)
+double Math::Laue(const double x, size_t N)
 {
     static const double SQRT6DOUBLE_EPS = std::sqrt(6.0 * std::numeric_limits<double>::epsilon());
     auto nd = static_cast<double>(N);
@@ -90,306 +84,28 @@ double MathFunctions::Laue(const double x, size_t N)
     return num / den;
 }
 
-double MathFunctions::erf(double arg)
+double Math::erf(double arg)
 {
     if (arg < 0.0)
-        throw std::runtime_error("Error in MathFunctions::erf: negative argument is not allowed");
+        throw std::runtime_error("Error in Math::erf: negative argument is not allowed");
     if (std::isinf(arg))
         return 1.0;
     return gsl_sf_erf(arg);
 }
 
-// ************************************************************************** //
-//  Bessel functions
-// ************************************************************************** //
-
-namespace MathFunctions
-{
-//! Computes complex Bessel function J0(z), using power series and asymptotic expansion
-complex_t Bessel_J0_PowSer(const complex_t z);
-
-//! Computes complex Bessel function J0(z), using power series and asymptotic expansion
-complex_t Bessel_J1_PowSer(const complex_t z);
-} // namespace MathFunctions
-
-double MathFunctions::Bessel_J0(double x)
-{
-    return gsl_sf_bessel_J0(x);
-}
-
-double MathFunctions::Bessel_J1(double x)
-{
-    return gsl_sf_bessel_J1(x);
-}
-
-double MathFunctions::Bessel_J1c(double x)
-{
-    return x == 0 ? 0.5 : gsl_sf_bessel_J1(x) / x;
-}
-
-double MathFunctions::Bessel_I0(double x)
-{
-    return gsl_sf_bessel_I0(x);
-}
-
-complex_t MathFunctions::Bessel_J0(const complex_t z)
-{
-    if (std::imag(z) == 0)
-        return gsl_sf_bessel_J0(std::real(z));
-    return Bessel_J0_PowSer(z);
-}
-
-complex_t MathFunctions::Bessel_J1(const complex_t z)
-{
-    if (std::imag(z) == 0)
-        return gsl_sf_bessel_J1(std::real(z));
-    return Bessel_J1_PowSer(z);
-}
-
-complex_t MathFunctions::Bessel_J1c(const complex_t z)
-{
-    if (std::imag(z) == 0) {
-        double xv = std::real(z);
-        return xv == 0 ? 0.5 : gsl_sf_bessel_J1(xv) / xv;
-    }
-    return z == 0. ? 0.5 : MathFunctions::Bessel_J1_PowSer(z) / z;
-}
-
-//! Computes the complex Bessel function J0(z), using power series and asymptotic expansion.
-//!
-//! Forked from unoptimized code at http://www.crbond.com/math.htm,
-//! who refers to "Computation of Special Functions", Zhang and Jin, John Wiley and Sons, 1996.
-
-complex_t MathFunctions::Bessel_J0_PowSer(const complex_t z)
-{
-    complex_t cj0;
-    static const double eps = 1e-15;
-    static double a[] = {-7.03125e-2,           0.112152099609375,     -0.5725014209747314,
-                         6.074042001273483,     -1.100171402692467e2,  3.038090510922384e3,
-                         -1.188384262567832e5,  6.252951493434797e6,   -4.259392165047669e8,
-                         3.646840080706556e10,  -3.833534661393944e12, 4.854014686852901e14,
-                         -7.286857349377656e16, 1.279721941975975e19};
-    static double b[] = {7.32421875e-2,         -0.2271080017089844,  1.727727502584457,
-                         -2.438052969955606e1,  5.513358961220206e2,  -1.825775547429318e4,
-                         8.328593040162893e5,   -5.006958953198893e7, 3.836255180230433e9,
-                         -3.649010818849833e11, 4.218971570284096e13, -5.827244631566907e15,
-                         9.476288099260110e17,  -1.792162323051699e20};
-
-    double a0 = std::abs(z);
-    complex_t z1 = z;
-    if (a0 == 0.0)
-        return 1.0;
-    if (std::real(z) < 0.0)
-        z1 = -z;
-    if (a0 <= 12.0) {
-        // standard power series [http://dlmf.nist.gov/10.2 (10.2.2)]
-        complex_t z2 = 0.25 * z * z;
-        cj0 = 1.0;
-        complex_t cr = 1.0;
-        for (size_t k = 1; k <= 40; ++k) {
-            cr *= -z2 / (double)(k * k);
-            cj0 += cr;
-            if (std::abs(cr) < std::abs(cj0) * eps)
-                break;
-        }
-    } else {
-        // Hankel's asymptotic expansion [http://dlmf.nist.gov/10.17 (10.17.3)]
-        size_t kz;
-        if (a0 >= 50.0)
-            kz = 8; // can be changed to 10
-        else if (a0 >= 35.0)
-            kz = 10; //   "      "     "  12
-        else
-            kz = 12; //   "      "     "  14
-        complex_t ct1 = z1 - M_PI_4;
-        complex_t cp0 = 1.0;
-        complex_t cq0 = -0.125;
-        const complex_t z1m2 = 1. / (z1 * z1); // faster than std::pow(z1, -2.0) ??
-        complex_t ptmp = z1m2;
-        for (size_t k = 0; k < kz; ++k) {
-            cp0 += a[k] * ptmp;
-            cq0 += b[k] * ptmp;
-            ptmp *= z1m2;
-        }
-        cj0 = std::sqrt(M_2_PI / z1) * (cp0 * std::cos(ct1) - cq0 / z1 * std::sin(ct1));
-    }
-    return cj0;
-}
-
-//! Computes the complex Bessel function J1(z), using power series and asymptotic expansion.
-//!
-//! Forked from same source as for Bessel_J0_PowSer
-
-complex_t MathFunctions::Bessel_J1_PowSer(const complex_t z)
-{
-    complex_t cj1;
-    static const double eps = 1e-15;
-
-    static double a1[] = {0.1171875,
-                          -0.1441955566406250,
-                          0.6765925884246826,
-                          -6.883914268109947,
-                          1.215978918765359e2,
-                          -3.302272294480852e3,
-                          1.276412726461746e5,
-                          -6.656367718817688e6,
-                          4.502786003050393e8,
-                          -3.833857520742790e10,
-                          4.011838599133198e12,
-                          -5.060568503314727e14,
-                          7.572616461117958e16,
-                          -1.326257285320556e19};
-    static double b1[] = {-0.1025390625,         0.2775764465332031,    -1.993531733751297,
-                          2.724882731126854e1,   -6.038440767050702e2,  1.971837591223663e4,
-                          -8.902978767070678e5,  5.310411010968522e7,   -4.043620325107754e9,
-                          3.827011346598605e11,  -4.406481417852278e13, 6.065091351222699e15,
-                          -9.833883876590679e17, 1.855045211579828e20};
-
-    double a0 = std::abs(z);
-    if (a0 == 0.0)
-        return 0.0;
-
-    complex_t z1 = z;
-    if (std::real(z) < 0.0)
-        z1 = -z;
-    if (a0 <= 12.0) {
-        // standard power series [http://dlmf.nist.gov/10.2 (10.2.2)]
-        const complex_t z2 = 0.25 * z * z;
-        cj1 = 1.0;
-        complex_t cr = 1.0; // powers will be computed recursively
-        for (int k = 1; k <= 40; ++k) {
-            cr *= -z2 / (double)(k * (k + 1));
-            cj1 += cr;
-            if (std::abs(cr) < std::abs(cj1) * eps)
-                break;
-        }
-        cj1 *= 0.5 * z1;
-    } else {
-        // Hankel's asymptotic expansion [http://dlmf.nist.gov/10.17 (10.17.3)]
-        size_t kz;
-        if (a0 >= 50.0)
-            kz = 8; // can be changed to 10
-        else if (a0 >= 35.0)
-            kz = 10; //   "      "     "  12
-        else
-            kz = 12; //   "      "     "  14
-        complex_t cp1 = 1.0;
-        complex_t cq1 = 0.375;                 // division by z1 postponed to final sum
-        const complex_t z1m2 = 1. / (z1 * z1); // faster than std::pow(z1, -2.0) ??
-        complex_t ptmp = z1m2;                 // powers will be computed recursively
-        for (size_t k = 0; k < kz; ++k) {
-            cp1 += a1[k] * ptmp;
-            cq1 += b1[k] * ptmp; // division by z1 postponed to final sum
-            ptmp *= z1m2;
-        }
-        const complex_t ct2 = z1 - 0.75 * M_PI;
-        cj1 = std::sqrt(M_2_PI / z1) * (cp1 * std::cos(ct2) - cq1 / z1 * std::sin(ct2));
-    }
-    if (std::real(z) < 0.0)
-        cj1 = -cj1;
-    return cj1;
-}
-
-// ************************************************************************** //
-//  Fourier transform and convolution
-// ************************************************************************** //
-
-//! @brief simple (and unoptimized) wrapper function
-//!   for the discrete fast Fourier transformation library (fftw3)
-
-std::vector<complex_t> MathFunctions::FastFourierTransform(const std::vector<complex_t>& data,
-                                                           MathFunctions::EFFTDirection ftCase)
-{
-    double scale(1.);
-    size_t npx = data.size();
-
-    fftw_complex* ftData = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * npx);
-    fftw_complex* ftResult = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * npx);
-    memset(ftData, 0, sizeof(fftw_complex) * npx);
-    memset(ftResult, 0, sizeof(fftw_complex) * npx);
-
-    for (size_t i = 0; i < npx; i++) {
-        ftData[i][0] = data[i].real();
-        ftData[i][1] = data[i].imag();
-    }
-
-    fftw_plan plan;
-    switch (ftCase) {
-    case MathFunctions::FORWARD_FFT:
-        plan = fftw_plan_dft_1d((int)npx, ftData, ftResult, FFTW_FORWARD, FFTW_ESTIMATE);
-        break;
-    case MathFunctions::BACKWARD_FFT:
-        plan = fftw_plan_dft_1d((int)npx, ftData, ftResult, FFTW_BACKWARD, FFTW_ESTIMATE);
-        scale = 1. / double(npx);
-        break;
-    default:
-        throw std::runtime_error(
-            "MathFunctions::FastFourierTransform -> Panic! Unknown transform case.");
-    }
-
-    fftw_execute(plan);
-
-    // saving data for user
-    std::vector<complex_t> outData;
-    outData.resize(npx);
-    for (size_t i = 0; i < npx; i++)
-        outData[i] = scale * complex_t(ftResult[i][0], ftResult[i][1]);
-
-    fftw_destroy_plan(plan);
-    fftw_free(ftData);
-    fftw_free(ftResult);
-
-    return outData;
-}
-
-//! @brief simple (and unoptimized) wrapper function
-//!   for the discrete fast Fourier transformation library (fftw3);
-//!   transforms real to complex
-
-std::vector<complex_t> MathFunctions::FastFourierTransform(const std::vector<double>& data,
-                                                           MathFunctions::EFFTDirection ftCase)
-{
-    std::vector<complex_t> cdata;
-    cdata.resize(data.size());
-    for (size_t i = 0; i < data.size(); i++)
-        cdata[i] = complex_t(data[i], 0);
-    return MathFunctions::FastFourierTransform(cdata, ftCase);
-}
-
-//! convolution of two real vectors of equal size
-
-std::vector<complex_t> MathFunctions::ConvolveFFT(const std::vector<double>& signal,
-                                                  const std::vector<double>& resfunc)
-{
-    if (signal.size() != resfunc.size())
-        throw std::runtime_error("MathFunctions::ConvolveFFT() -> This convolution works only for "
-                                 "two vectors of equal size. Use Convolve class instead.");
-    std::vector<complex_t> fft_signal =
-        MathFunctions::FastFourierTransform(signal, MathFunctions::FORWARD_FFT);
-    std::vector<complex_t> fft_resfunc =
-        MathFunctions::FastFourierTransform(resfunc, MathFunctions::FORWARD_FFT);
-
-    std::vector<complex_t> fft_prod;
-    fft_prod.resize(fft_signal.size());
-    for (size_t i = 0; i < fft_signal.size(); i++)
-        fft_prod[i] = fft_signal[i] * fft_resfunc[i];
-
-    std::vector<complex_t> result =
-        MathFunctions::FastFourierTransform(fft_prod, MathFunctions::BACKWARD_FFT);
-    return result;
-}
-
 // ************************************************************************** //
 //  Random number generators
 // ************************************************************************** //
 
-double MathFunctions::GenerateUniformRandom()
+/* currently unused
+
+double Math::GenerateUniformRandom()
 {
     int random_int = std::rand();
     return (double)random_int / RAND_MAX;
 }
 
-double MathFunctions::GenerateStandardNormalRandom() // using c++11 standard library
+double Math::GenerateStandardNormalRandom() // using c++11 standard library
 {
     unsigned seed =
         static_cast<unsigned>(std::chrono::system_clock::now().time_since_epoch().count());
@@ -398,12 +114,13 @@ double MathFunctions::GenerateStandardNormalRandom() // using c++11 standard lib
     return distribution(generator);
 }
 
-double MathFunctions::GenerateNormalRandom(double average, double std_dev)
+double Math::GenerateNormalRandom(double average, double std_dev)
 {
     return GenerateStandardNormalRandom() * std_dev + average;
 }
+*/
 
-double MathFunctions::GeneratePoissonRandom(double average) // using c++11 standard library
+double Math::GeneratePoissonRandom(double average) // using c++11 standard library
 {
     unsigned seed =
         static_cast<unsigned>(std::chrono::system_clock::now().time_since_epoch().count());
diff --git a/Base/Utils/MathFunctions.h b/Base/Utils/MathFunctions.h
index 727ae2c390f9d07c99178e81eb5e7e81b17b6802..3f2ab498e311e0f5164f311ccd3476e11fcc1668 100644
--- a/Base/Utils/MathFunctions.h
+++ b/Base/Utils/MathFunctions.h
@@ -3,7 +3,7 @@
 //  BornAgain: simulate and fit scattering at grazing incidence
 //
 //! @file      Base/Utils/MathFunctions.h
-//! @brief     Defines namespace MathFunctions.
+//! @brief     Defines functions in namespace Math.
 //!
 //! @homepage  http://www.bornagainproject.org
 //! @license   GNU General Public License v3 or higher (see COPYING)
@@ -20,7 +20,7 @@
 
 //! Various mathematical functions.
 
-namespace MathFunctions
+namespace Math
 {
 
 // ************************************************************************** //
@@ -34,9 +34,6 @@ double IntegratedGaussian(double x, double average, double std_dev);
 //! cotangent function: \f$cot(x)\equiv1/tan(x)\f$
 double cot(double x);
 
-//! Sine integral function: \f$Si(x)\equiv\int_0^x du \sin(u)/u\f$
-double Si(double x);
-
 //! sinc function: \f$sinc(x)\equiv\sin(x)/x\f$
 double sinc(double x);
 
@@ -52,56 +49,15 @@ double Laue(const double x, size_t N);
 //! Error function of real-valued argument
 double erf(double arg);
 
-// ************************************************************************** //
-//  Bessel functions
-// ************************************************************************** //
-
-//! Bessel function of the first kind and order 0
-double Bessel_J0(double x);
-
-//! Bessel function of the first kind and order 1
-double Bessel_J1(double x);
-
-//! Bessel function  Bessel_J1(x)/x
-double Bessel_J1c(double x);
-
-//! Modified Bessel function of the first kind and order 0
-double Bessel_I0(double x);
-
-//! Complex Bessel function of the first kind and order 0
-complex_t Bessel_J0(const complex_t z);
-
-//! Complex Bessel function of the first kind and order 1
-complex_t Bessel_J1(const complex_t z);
-
-//! Complex Bessel function  Bessel_J1(x)/x
-complex_t Bessel_J1c(const complex_t z);
-
-// ************************************************************************** //
-//  Fourier transform and convolution
-// ************************************************************************** //
-
-// TODO move elsewhere, and rm #include <vector>
-
-enum EFFTDirection { FORWARD_FFT, BACKWARD_FFT };
-
-// TODO: name these two functions differently (SWIG warning 509)
-std::vector<complex_t> FastFourierTransform(const std::vector<complex_t>& data,
-                                            EFFTDirection tcase);
-std::vector<complex_t> FastFourierTransform(const std::vector<double>& data, EFFTDirection tcase);
-
-std::vector<complex_t> ConvolveFFT(const std::vector<double>& signal,
-                                   const std::vector<double>& resfunc);
-
 // ************************************************************************** //
 //  Random number generators
 // ************************************************************************** //
 
-double GenerateUniformRandom();
-double GenerateStandardNormalRandom();
-double GenerateNormalRandom(double average, double std_dev);
+//double GenerateUniformRandom();
+//double GenerateStandardNormalRandom();
+//double GenerateNormalRandom(double average, double std_dev);
 double GeneratePoissonRandom(double average);
 
-} // Namespace MathFunctions
+} // Namespace Math
 
 #endif // BORNAGAIN_BASE_UTILS_MATHFUNCTIONS_H
diff --git a/Base/Vector/BasicVector3D.cpp b/Base/Vector/BasicVector3D.cpp
index 2de37e14714d51829eca13aad665e31e3afc7235..759ad873a9433f5f0e458450416a669bb2fd5a72 100644
--- a/Base/Vector/BasicVector3D.cpp
+++ b/Base/Vector/BasicVector3D.cpp
@@ -14,7 +14,7 @@
 
 #include "Base/Vector/BasicVector3D.h"
 #include "Base/Const/MathConstants.h"
-#include "Base/Types/Exceptions.h"
+#include <stdexcept>
 
 typedef std::complex<double> complex_t;
 
@@ -90,7 +90,7 @@ template <> BasicVector3D<double> BasicVector3D<double>::unit() const
 {
     double len = mag();
     if (len == 0.0)
-        throw Exceptions::DivisionByZeroException("Cannot normalize zero vector");
+        throw std::runtime_error("Cannot normalize zero vector");
     return BasicVector3D<double>(x() / len, y() / len, z() / len);
 }
 
@@ -98,7 +98,7 @@ template <> BasicVector3D<complex_t> BasicVector3D<complex_t>::unit() const
 {
     double len = mag();
     if (len == 0.0)
-        throw Exceptions::DivisionByZeroException("Cannot normalize zero vector");
+        throw std::runtime_error("Cannot normalize zero vector");
     return BasicVector3D<complex_t>(x() / len, y() / len, z() / len);
 }
 
diff --git a/Base/Vector/Transform3D.cpp b/Base/Vector/Transform3D.cpp
index 09bba513cd19bb67cc0cb5ea38878d29c86a7c6b..a52f38962e835f4092f94bf106d9379ce1920836 100644
--- a/Base/Vector/Transform3D.cpp
+++ b/Base/Vector/Transform3D.cpp
@@ -26,11 +26,6 @@ Transform3D::Transform3D(const Eigen::Matrix3d& matrix) : m_matrix(matrix)
     m_inverse_matrix = m_matrix.inverse();
 }
 
-Transform3D Transform3D::createIdentity()
-{
-    return Transform3D();
-}
-
 Transform3D Transform3D::createRotateX(double phi)
 {
     double cosine = std::cos(phi);
diff --git a/Base/Vector/Transform3D.h b/Base/Vector/Transform3D.h
index eac11d2670a1cf54922568b4edc9d98b6245e86d..8eb54d84c095cd959306cdea6396d8a38a45e323 100644
--- a/Base/Vector/Transform3D.h
+++ b/Base/Vector/Transform3D.h
@@ -43,9 +43,6 @@ public:
     //! Clones the transformation
     Transform3D* clone() const;
 
-    //! Creates identity transformation (default)
-    static Transform3D createIdentity();
-
     //! Creates rotation around x-axis
     static Transform3D createRotateX(double phi);
 
diff --git a/Core/Computation/DWBASingleComputation.cpp b/Core/Computation/DWBASingleComputation.cpp
index e3acd2172c4373c3c1631c0678c5abbdede2f924..fc183db8492c4f50c1ae52515e69cae239d25371 100644
--- a/Core/Computation/DWBASingleComputation.cpp
+++ b/Core/Computation/DWBASingleComputation.cpp
@@ -22,8 +22,6 @@ DWBASingleComputation::DWBASingleComputation() = default;
 
 DWBASingleComputation::~DWBASingleComputation() = default;
 
-DWBASingleComputation::DWBASingleComputation(DWBASingleComputation&&) = default;
-
 void DWBASingleComputation::setProgressHandler(ProgressHandler* p_progress)
 {
     m_progress_counter = std::make_unique<DelayedProgressCounter>(p_progress, 100);
@@ -59,8 +57,3 @@ void DWBASingleComputation::compute(SimulationElement& elem) const
     if (m_progress_counter)
         m_progress_counter->stepProgress();
 }
-
-const std::map<size_t, std::vector<HomogeneousRegion>>& DWBASingleComputation::regionMap() const
-{
-    return m_region_map;
-}
diff --git a/Core/Computation/DWBASingleComputation.h b/Core/Computation/DWBASingleComputation.h
index dc00b0ffbce62534a0189b4716277832c1937545..424c755c3e7c543f5ec8a5de68a8f63dab992828 100644
--- a/Core/Computation/DWBASingleComputation.h
+++ b/Core/Computation/DWBASingleComputation.h
@@ -39,7 +39,6 @@ class DWBASingleComputation
 public:
     DWBASingleComputation();
     ~DWBASingleComputation();
-    DWBASingleComputation(DWBASingleComputation&& other);
 
     void setProgressHandler(ProgressHandler* p_progress);
 
@@ -48,9 +47,6 @@ public:
     void setSpecularBinComputation(GISASSpecularComputation* p_spec_comp);
     void compute(SimulationElement& elem) const;
 
-    //! Retrieves a map of regions for the calculation of averaged layers
-    const std::map<size_t, std::vector<HomogeneousRegion>>& regionMap() const;
-
 private:
     std::vector<std::unique_ptr<ParticleLayoutComputation>> m_layout_comps;
     std::unique_ptr<RoughMultiLayerComputation> m_roughness_comp;
diff --git a/Core/Computation/PoissonNoiseBackground.cpp b/Core/Computation/PoissonNoiseBackground.cpp
index 1f2c7637c1205b0371e5f5365558bec35a750087..7fb8924744e36074fbde403637fbdbfa8eb1f86f 100644
--- a/Core/Computation/PoissonNoiseBackground.cpp
+++ b/Core/Computation/PoissonNoiseBackground.cpp
@@ -27,5 +27,5 @@ PoissonNoiseBackground* PoissonNoiseBackground::clone() const
 
 double PoissonNoiseBackground::addBackground(double intensity) const
 {
-    return MathFunctions::GeneratePoissonRandom(intensity);
+    return Math::GeneratePoissonRandom(intensity);
 }
diff --git a/Core/Computation/ProfileHelper.cpp b/Core/Computation/ProfileHelper.cpp
index a402bda3993fa9ac3b37b0b8e785c9e1e7947ad5..298a9370134d6f4df7014be59c5cda9bad2450d2 100644
--- a/Core/Computation/ProfileHelper.cpp
+++ b/Core/Computation/ProfileHelper.cpp
@@ -19,8 +19,16 @@
 namespace
 {
 const double prefactor = std::sqrt(2.0 / M_PI);
-double Transition(double x, double sigma);
-double TransitionTanh(double x);
+double TransitionTanh(double x)
+{
+    return (1.0 - std::tanh(prefactor * x)) / 2.0;
+}
+double Transition(double x, double sigma)
+{
+    if (sigma <= 0.0)
+        return x < 0.0 ? 1.0 : 0.0;
+    return TransitionTanh(x / sigma);
+}
 } // namespace
 
 ProfileHelper::ProfileHelper(const ProcessedSample& sample)
@@ -45,6 +53,8 @@ ProfileHelper::ProfileHelper(const ProcessedSample& sample)
     }
 }
 
+ProfileHelper::~ProfileHelper() = default;
+
 // Note: for refractive index materials, the material interpolation actually happens at the level
 // of n^2. To first order in delta and beta, this implies the same smooth interpolation of delta
 // and beta, as is done here.
@@ -75,19 +85,3 @@ std::pair<double, double> ProfileHelper::defaultLimits() const
     double z_max = m_zlimits.front() + top_margin;
     return {z_min, z_max};
 }
-
-ProfileHelper::~ProfileHelper() = default;
-
-namespace
-{
-double Transition(double x, double sigma)
-{
-    if (sigma <= 0.0)
-        return x < 0.0 ? 1.0 : 0.0;
-    return TransitionTanh(x / sigma);
-}
-double TransitionTanh(double x)
-{
-    return (1.0 - std::tanh(prefactor * x)) / 2.0;
-}
-} // namespace
diff --git a/Core/Simulation/Simulation.cpp b/Core/Simulation/Simulation.cpp
index 02f10d87e5014c09bf32d81edf012d6fd3f60a6d..097d113d6af19571460b9f29bbeaff4262e8b23d 100644
--- a/Core/Simulation/Simulation.cpp
+++ b/Core/Simulation/Simulation.cpp
@@ -158,11 +158,6 @@ void Simulation::setDetectorResolutionFunction(const IResolutionFunction2D& reso
     instrument().setDetectorResolutionFunction(resolution_function);
 }
 
-void Simulation::removeDetectorResolutionFunction()
-{
-    instrument().removeDetectorResolution();
-}
-
 //! Sets the polarization analyzer characteristics of the detector
 void Simulation::setAnalyzerProperties(const kvector_t direction, double efficiency,
                                        double total_transmission)
diff --git a/Core/Simulation/Simulation.h b/Core/Simulation/Simulation.h
index 8dda1d7137ccbe841f16b22260c74466d83c1ca5..e3d6c5e2c26e8d094dfda2451dd7d10b69ff19df 100644
--- a/Core/Simulation/Simulation.h
+++ b/Core/Simulation/Simulation.h
@@ -61,7 +61,6 @@ public:
     void setBeamPolarization(const kvector_t bloch_vector);
 
     void setDetectorResolutionFunction(const IResolutionFunction2D& resolution_function);
-    void removeDetectorResolutionFunction();
 
     void setAnalyzerProperties(const kvector_t direction, double efficiency,
                                double total_transmission);
diff --git a/Device/Beam/FootprintGauss.cpp b/Device/Beam/FootprintGauss.cpp
index 23870246631cfff628c2c77ab3b5ae8d5e25e1e3..8d4ee20787f4412c9040aa578752e40166407bd4 100644
--- a/Device/Beam/FootprintGauss.cpp
+++ b/Device/Beam/FootprintGauss.cpp
@@ -39,7 +39,7 @@ double FootprintGauss::calculate(double alpha) const
     if (widthRatio() == 0.0)
         return 1.0;
     const double arg = std::sin(alpha) * M_SQRT1_2 / widthRatio();
-    return MathFunctions::erf(arg);
+    return Math::erf(arg);
 }
 
 std::string FootprintGauss::print() const
diff --git a/Device/Histo/IHistogram.cpp b/Device/Histo/IHistogram.cpp
index fc37205a696ea33b138a27acd6d656699a980eba..2402d02b4d0294d000c0d2ed144a094bd364ca7d 100644
--- a/Device/Histo/IHistogram.cpp
+++ b/Device/Histo/IHistogram.cpp
@@ -123,6 +123,16 @@ double IHistogram::getYaxisValue(size_t i)
     return m_data.getAxisValue(i, 1);
 }
 
+const OutputData<CumulativeValue>& IHistogram::getData() const
+{
+    return m_data;
+}
+
+OutputData<CumulativeValue>& IHistogram::getData()
+{
+    return m_data;
+}
+
 double IHistogram::getBinContent(size_t i) const
 {
     return m_data[i].getContent();
diff --git a/Device/Histo/IHistogram.h b/Device/Histo/IHistogram.h
index f27aeef472ca6ebcae599a105ce3f8fc8483cf53..614aeb55df9c37bca32dbcc7e34acf6bbbfdb37c 100644
--- a/Device/Histo/IHistogram.h
+++ b/Device/Histo/IHistogram.h
@@ -30,7 +30,7 @@ public:
 
     IHistogram();
     IHistogram(const IHistogram& other);
-    virtual ~IHistogram() {}
+    virtual ~IHistogram() = default;
 
     IHistogram(const IAxis& axis_x);
     IHistogram(const IAxis& axis_x, const IAxis& axis_y);
@@ -86,12 +86,12 @@ public:
     //! Returns the center of bin i of the y axis.
     double getYaxisValue(size_t i);
 
+    const OutputData<CumulativeValue>& getData() const;
+    OutputData<CumulativeValue>& getData();
+
     //! Returns content (accumulated value) of bin i.
     double getBinContent(size_t i) const;
 
-    const OutputData<CumulativeValue>& getData() const { return m_data; }
-    OutputData<CumulativeValue>& getData() { return m_data; }
-
     //! Returns content (accumulated value) of the 2D histogram bin.
     double getBinContent(size_t binx, size_t biny) const;
 
diff --git a/Device/Histo/SimulationResult.cpp b/Device/Histo/SimulationResult.cpp
index 5108a31b2ef35f8d2d706b778b7d1d844e49e23c..5baf590c6e04d22735f675e0f779b28d6f3ee685 100644
--- a/Device/Histo/SimulationResult.cpp
+++ b/Device/Histo/SimulationResult.cpp
@@ -114,8 +114,8 @@ double SimulationResult::max() const
 {
     ASSERT(m_data);
     double result = 0;
-    for (size_t i=0; i<size(); ++i)
-        if ((*m_data)[i]>result)
+    for (size_t i = 0; i < size(); ++i)
+        if ((*m_data)[i] > result)
             result = (*m_data)[i];
     return result;
 }
diff --git a/Device/InputOutput/OutputDataReadStrategy.h b/Device/InputOutput/OutputDataReadStrategy.h
index 07a57fc273e53dc84c80b4c8ebba7dcb47bc97fc..39d7153af76172ec9463f73b6b5f20d41efd4728 100644
--- a/Device/InputOutput/OutputDataReadStrategy.h
+++ b/Device/InputOutput/OutputDataReadStrategy.h
@@ -25,7 +25,7 @@ template <class T> class OutputData;
 class IOutputDataReadStrategy
 {
 public:
-    virtual ~IOutputDataReadStrategy() {}
+    virtual ~IOutputDataReadStrategy() = default;
     virtual OutputData<double>* readOutputData(std::istream& input_stream) = 0;
 };
 
diff --git a/Device/InputOutput/OutputDataWriteStrategy.h b/Device/InputOutput/OutputDataWriteStrategy.h
index e38b0d4e9c223b5ab5f2954c14d671fb2fba21aa..81168c5f583a8a02d0cd128d5707905ff068c8c2 100644
--- a/Device/InputOutput/OutputDataWriteStrategy.h
+++ b/Device/InputOutput/OutputDataWriteStrategy.h
@@ -25,8 +25,7 @@ template <class T> class OutputData;
 class IOutputDataWriteStrategy
 {
 public:
-    IOutputDataWriteStrategy() {}
-    virtual ~IOutputDataWriteStrategy() {}
+    virtual ~IOutputDataWriteStrategy() = default;
 
     virtual void writeOutputData(const OutputData<double>& data, std::ostream& output_stream) = 0;
 };
diff --git a/Device/Instrument/FourierTransform.h b/Device/Instrument/FourierTransform.h
index 31a9c152e411a32dcaeb325fd7fd07afb4dc64d0..b7693e78d0b44cf980bca1514341d282f9d0e0a3 100644
--- a/Device/Instrument/FourierTransform.h
+++ b/Device/Instrument/FourierTransform.h
@@ -3,7 +3,7 @@
 //  BornAgain: simulate and fit scattering at grazing incidence
 //
 //! @file      Device/Instrument/FourierTransform.h
-//! @brief     Defines class MathFunctions::FourierTransform.
+//! @brief     Defines class Math::FourierTransform.
 //!
 //! @homepage  http://www.bornagainproject.org
 //! @license   GNU General Public License v3 or higher (see COPYING)
diff --git a/Device/Resolution/Convolve.h b/Device/Resolution/Convolve.h
index 93937d9917a124f7e5b421df46bd2c6fe1267602..795d29b122c1b0efbe80bb1e79a809f9f3a195f9 100644
--- a/Device/Resolution/Convolve.h
+++ b/Device/Resolution/Convolve.h
@@ -3,7 +3,7 @@
 //  BornAgain: simulate and fit scattering at grazing incidence
 //
 //! @file      Device/Resolution/Convolve.h
-//! @brief     Defines class MathFunctions::Convolve.
+//! @brief     Defines class Math::Convolve.
 //!
 //! @homepage  http://www.bornagainproject.org
 //! @license   GNU General Public License v3 or higher (see COPYING)
diff --git a/Device/Resolution/ResolutionFunction2DGaussian.cpp b/Device/Resolution/ResolutionFunction2DGaussian.cpp
index 4fc44e516da5366ffa050b940192f1a3fe1caf33..08cfbd3cadd7fa82321bef636f11f5bdcfe1cc2d 100644
--- a/Device/Resolution/ResolutionFunction2DGaussian.cpp
+++ b/Device/Resolution/ResolutionFunction2DGaussian.cpp
@@ -26,6 +26,6 @@ ResolutionFunction2DGaussian::ResolutionFunction2DGaussian(double sigma_x, doubl
 
 double ResolutionFunction2DGaussian::evaluateCDF(double x, double y) const
 {
-    return MathFunctions::IntegratedGaussian(x, 0.0, m_sigma_x)
-           * MathFunctions::IntegratedGaussian(y, 0.0, m_sigma_y);
+    return Math::IntegratedGaussian(x, 0.0, m_sigma_x)
+           * Math::IntegratedGaussian(y, 0.0, m_sigma_y);
 }
diff --git a/Param/Base/ParameterPool.h b/Param/Base/ParameterPool.h
index a77d7ef7d70279fe37f6076bb282bacb16d2a2ac..dc06f8f7bf3f3d384bfd81f9a2ce59c983873b27 100644
--- a/Param/Base/ParameterPool.h
+++ b/Param/Base/ParameterPool.h
@@ -69,8 +69,8 @@ public:
 
     void removeParameter(const std::string& name);
 
-    const RealParameter* operator[](size_t index) const;
-    RealParameter* operator[](size_t index);
+    RealParameter* operator[](size_t index); // used by libBornAgainParam.i
+    const RealParameter* operator[](size_t index) const; // needed by the above
 
 private:
     virtual void print(std::ostream& ostr) const;
diff --git a/Param/Node/INodeVisitor.cpp b/Param/Node/INodeVisitor.cpp
index a4149896f17f12332a1c6d0d18dba89021cb20c5..288d0930e8d060ab9d94cb8112a02c2d6f3a3575 100644
--- a/Param/Node/INodeVisitor.cpp
+++ b/Param/Node/INodeVisitor.cpp
@@ -26,15 +26,3 @@ void VisitNodesPreorder(const INode& node, INodeVisitor& visitor)
         it.next();
     }
 }
-
-void VisitNodesPostorder(const INode& node, INodeVisitor& visitor)
-{
-    NodeIterator<PostorderStrategy> it(&node);
-    it.first();
-    while (!it.isDone()) {
-        visitor.setDepth(it.depth());
-        const INode* child = it.getCurrent();
-        child->accept(&visitor);
-        it.next();
-    }
-}
diff --git a/Param/Node/INodeVisitor.h b/Param/Node/INodeVisitor.h
index 8358ace6a2fbda94214e1c3aa502eac9dbc87fab..ea14ac42d3788c08a3555da5befbf19f6b04952e 100644
--- a/Param/Node/INodeVisitor.h
+++ b/Param/Node/INodeVisitor.h
@@ -287,6 +287,5 @@ private:
 };
 
 void VisitNodesPreorder(const INode& node, INodeVisitor& visitor);
-void VisitNodesPostorder(const INode& node, INodeVisitor& visitor);
 
 #endif // BORNAGAIN_PARAM_NODE_INODEVISITOR_H
diff --git a/Param/Node/IterationStrategy.cpp b/Param/Node/IterationStrategy.cpp
index 12aceb9a3fe61fd79bb05c2abeeb5aeef9404372..91ef77c0444943516be01334f550778a1d0914a5 100644
--- a/Param/Node/IterationStrategy.cpp
+++ b/Param/Node/IterationStrategy.cpp
@@ -51,41 +51,3 @@ bool PreorderStrategy::isDone(IteratorMemento& iterator_stack) const
 {
     return iterator_stack.empty();
 }
-
-PostorderStrategy::PostorderStrategy() = default;
-
-PostorderStrategy* PostorderStrategy::clone() const
-{
-    return new PostorderStrategy();
-}
-
-IteratorMemento PostorderStrategy::first(const INode* p_root)
-{
-    IteratorMemento iterator_stack;
-    iterator_stack.push_state(IteratorState(p_root));
-    std::vector<const INode*> children = p_root->getChildren();
-    while (children.size() > 0) {
-        iterator_stack.push_state(IteratorState(children));
-        children = iterator_stack.getCurrent()->getChildren();
-    }
-    return iterator_stack;
-}
-
-void PostorderStrategy::next(IteratorMemento& iterator_stack) const
-{
-    iterator_stack.next();
-    if (iterator_stack.get_state().isEnd()) {
-        iterator_stack.pop_state();
-        return;
-    }
-    std::vector<const INode*> children = iterator_stack.getCurrent()->getChildren();
-    while (children.size() > 0) {
-        iterator_stack.push_state(IteratorState(children));
-        children = iterator_stack.getCurrent()->getChildren();
-    }
-}
-
-bool PostorderStrategy::isDone(IteratorMemento& iterator_stack) const
-{
-    return iterator_stack.empty();
-}
diff --git a/Param/Node/IterationStrategy.h b/Param/Node/IterationStrategy.h
index af52ca26817444738d0a1a86b7890c2045027adf..b66eaa0c216f1ea3e23c4a3ea21979dcbffd683d 100644
--- a/Param/Node/IterationStrategy.h
+++ b/Param/Node/IterationStrategy.h
@@ -45,17 +45,4 @@ public:
     virtual bool isDone(IteratorMemento& iterator_stack) const;
 };
 
-//! Traverse tree; visit children before their parents.
-class PostorderStrategy : public IterationStrategy
-{
-public:
-    PostorderStrategy();
-
-    virtual PostorderStrategy* clone() const;
-
-    virtual IteratorMemento first(const INode* p_root);
-    virtual void next(IteratorMemento& iterator_stack) const;
-    virtual bool isDone(IteratorMemento& iterator_stack) const;
-};
-
 #endif // BORNAGAIN_PARAM_NODE_ITERATIONSTRATEGY_H
diff --git a/Sample/Aggregate/InterferenceFunction2DSuperLattice.cpp b/Sample/Aggregate/InterferenceFunction2DSuperLattice.cpp
index 37851e1ad1a64ec9dce407926ad6b330af7156f0..4749ff677e35050f8103ed82f92739b6791e40bc 100644
--- a/Sample/Aggregate/InterferenceFunction2DSuperLattice.cpp
+++ b/Sample/Aggregate/InterferenceFunction2DSuperLattice.cpp
@@ -124,7 +124,7 @@ std::vector<const INode*> InterferenceFunction2DSuperLattice::getChildren() cons
 
 double InterferenceFunction2DSuperLattice::iff_without_dw(const kvector_t q) const
 {
-    using MathFunctions::Laue;
+    using Math::Laue;
 
     const double a = m_lattice->length1();
     const double b = m_lattice->length2();
diff --git a/Sample/Aggregate/InterferenceFunctionFinite2DLattice.cpp b/Sample/Aggregate/InterferenceFunctionFinite2DLattice.cpp
index 3e63690f75b9f62615496bc152d5f8f348a6fa34..ef2f13ab6be1094db9e7a6942d1a6949373e0a1b 100644
--- a/Sample/Aggregate/InterferenceFunctionFinite2DLattice.cpp
+++ b/Sample/Aggregate/InterferenceFunctionFinite2DLattice.cpp
@@ -21,7 +21,7 @@
 
 #include <limits>
 
-using MathFunctions::Laue;
+using Math::Laue;
 
 //! Constructor of two-dimensional finite lattice interference function.
 //! @param lattice: object specifying a 2d lattice structure
diff --git a/Sample/Aggregate/InterferenceFunctionFinite3DLattice.cpp b/Sample/Aggregate/InterferenceFunctionFinite3DLattice.cpp
index b6021ac32cdae6e784bf05c4a04d2adb4e49b0df..569eeecb300e353403403d73a5c3456c840a17c1 100644
--- a/Sample/Aggregate/InterferenceFunctionFinite3DLattice.cpp
+++ b/Sample/Aggregate/InterferenceFunctionFinite3DLattice.cpp
@@ -53,7 +53,7 @@ std::vector<const INode*> InterferenceFunctionFinite3DLattice::getChildren() con
 
 double InterferenceFunctionFinite3DLattice::iff_without_dw(const kvector_t q) const
 {
-    using MathFunctions::Laue;
+    using Math::Laue;
     const double qadiv2 = q.dot(m_lattice->getBasisVectorA()) / 2.0;
     const double qbdiv2 = q.dot(m_lattice->getBasisVectorB()) / 2.0;
     const double qcdiv2 = q.dot(m_lattice->getBasisVectorC()) / 2.0;
diff --git a/Sample/Aggregate/InterferenceFunctionHardDisk.cpp b/Sample/Aggregate/InterferenceFunctionHardDisk.cpp
index b679b3028850ceb23571d6cd6dcf9f6e39da9358..9e94a765d48cebbd32e0347ebc31f941eee31811 100644
--- a/Sample/Aggregate/InterferenceFunctionHardDisk.cpp
+++ b/Sample/Aggregate/InterferenceFunctionHardDisk.cpp
@@ -13,8 +13,8 @@
 // ************************************************************************** //
 
 #include "Sample/Aggregate/InterferenceFunctionHardDisk.h"
+#include "Base/Utils/Bessel.h"
 #include "Base/Utils/Integrator.h"
-#include "Base/Utils/MathFunctions.h"
 #include "Param/Base/RealParameter.h"
 #include <cmath>
 
@@ -83,7 +83,7 @@ double InterferenceFunctionHardDisk::packingRatio() const
 double InterferenceFunctionHardDisk::integrand(double x) const
 {
     double cx = m_c_zero * (1.0 + 4.0 * m_packing * (W2(x / 2.0) - 1.0) + m_s2 * x);
-    return x * cx * MathFunctions::Bessel_J0(m_q * x);
+    return x * cx * Math::Bessel::J0(m_q * x);
 }
 
 namespace
diff --git a/Sample/Correlations/FTDecay1D.cpp b/Sample/Correlations/FTDecay1D.cpp
index bf24738878fcc5afb882fff7c61ef8a1ce78be39..88a5ca16f3ee48711539e9030ef7e5e23327e52c 100644
--- a/Sample/Correlations/FTDecay1D.cpp
+++ b/Sample/Correlations/FTDecay1D.cpp
@@ -97,7 +97,7 @@ FTDecayFunction1DTriangle* FTDecayFunction1DTriangle::clone() const
 
 double FTDecayFunction1DTriangle::evaluate(double q) const
 {
-    double sincqw2 = MathFunctions::sinc(q * m_decay_length / 2.0);
+    double sincqw2 = Math::sinc(q * m_decay_length / 2.0);
     return m_decay_length * sincqw2 * sincqw2;
 }
 
diff --git a/Sample/Correlations/FTDistributions1D.cpp b/Sample/Correlations/FTDistributions1D.cpp
index 18f5ca2042b2370e45429d3253eaebb5942188ea..835cb691205b0b8c7b18ed103a09ca3b394b6b2e 100644
--- a/Sample/Correlations/FTDistributions1D.cpp
+++ b/Sample/Correlations/FTDistributions1D.cpp
@@ -124,7 +124,7 @@ FTDistribution1DGate* FTDistribution1DGate::clone() const
 
 double FTDistribution1DGate::evaluate(double q) const
 {
-    return MathFunctions::sinc(q * m_omega);
+    return Math::sinc(q * m_omega);
 }
 
 double FTDistribution1DGate::qSecondDerivative() const
@@ -158,7 +158,7 @@ FTDistribution1DTriangle* FTDistribution1DTriangle::clone() const
 
 double FTDistribution1DTriangle::evaluate(double q) const
 {
-    double sincqw2 = MathFunctions::sinc(q * m_omega / 2.0);
+    double sincqw2 = Math::sinc(q * m_omega / 2.0);
     return sincqw2 * sincqw2;
 }
 
@@ -196,7 +196,7 @@ double FTDistribution1DCosine::evaluate(double q) const
     double qw = std::abs(q * m_omega);
     if (std::abs(1.0 - qw * qw / M_PI / M_PI) < std::numeric_limits<double>::epsilon())
         return 0.5;
-    return MathFunctions::sinc(qw) / (1.0 - qw * qw / M_PI / M_PI);
+    return Math::sinc(qw) / (1.0 - qw * qw / M_PI / M_PI);
 }
 
 double FTDistribution1DCosine::qSecondDerivative() const
diff --git a/Sample/Correlations/FTDistributions2D.cpp b/Sample/Correlations/FTDistributions2D.cpp
index 624165fdf9512650ba9853bf4f71b84d14484556..8e1f54a1b40faa046072bf8b2c9e651f13df6aac 100644
--- a/Sample/Correlations/FTDistributions2D.cpp
+++ b/Sample/Correlations/FTDistributions2D.cpp
@@ -14,8 +14,8 @@
 
 #include "Sample/Correlations/FTDistributions2D.h"
 #include "Base/Types/Exceptions.h"
+#include "Base/Utils/Bessel.h"
 #include "Base/Utils/Integrator.h"
-#include "Base/Utils/MathFunctions.h"
 #include <limits>
 
 // ************************************************************************** //
@@ -119,7 +119,7 @@ FTDistribution2DGate* FTDistribution2DGate::clone() const
 double FTDistribution2DGate::evaluate(double qx, double qy) const
 {
     double scaled_q = std::sqrt(sumsq(qx, qy));
-    return MathFunctions::Bessel_J1c(scaled_q) * 2.0;
+    return Math::Bessel::J1c(scaled_q) * 2.0;
 }
 
 std::unique_ptr<IDistribution2DSampler> FTDistribution2DGate::createSampler() const
@@ -153,8 +153,8 @@ double FTDistribution2DCone::evaluate(double qx, double qy) const
         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 = 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);
+        [](double x) -> double { return x * x * Math::Bessel::J0(x); }, 0.0, scaled_q);
+    return 6.0 * (Math::Bessel::J1c(scaled_q) - integral / scaled_q / scaled_q / scaled_q);
 }
 
 std::unique_ptr<IDistribution2DSampler> FTDistribution2DCone::createSampler() const
diff --git a/Sample/Correlations/IPeakShape.cpp b/Sample/Correlations/IPeakShape.cpp
index 495d4e2484dc3ff076cbb4846cb42c8c301d2283..f4965089bac08f4736c4de2cbcbfd91eff664e27 100644
--- a/Sample/Correlations/IPeakShape.cpp
+++ b/Sample/Correlations/IPeakShape.cpp
@@ -14,8 +14,8 @@
 
 #include "Sample/Correlations/IPeakShape.h"
 #include "Base/Const/MathConstants.h"
+#include "Base/Utils/Bessel.h"
 #include "Base/Utils/Integrator.h"
-#include "Base/Utils/MathFunctions.h"
 #include <limits>
 
 namespace
@@ -56,7 +56,7 @@ double MisesPrefactor(double kappa)
     if (kappa > maxkappa2) {
         return std::sqrt(kappa / 2.0 / M_PI) / (1.0 + 1.0 / (8.0 * kappa));
     } else {
-        return std::exp(kappa) / (2.0 * M_PI * MathFunctions::Bessel_I0(kappa));
+        return std::exp(kappa) / (2.0 * M_PI * Math::Bessel::I0(kappa));
     }
 }
 
diff --git a/Sample/HardParticle/FormFactorAnisoPyramid.cpp b/Sample/HardParticle/FormFactorAnisoPyramid.cpp
index 3edfdd5e3122f07cd76a20909551f0c05979aa8a..c0b7c992a1d067f58d3be1ea6456b392bdfe591f 100644
--- a/Sample/HardParticle/FormFactorAnisoPyramid.cpp
+++ b/Sample/HardParticle/FormFactorAnisoPyramid.cpp
@@ -49,7 +49,7 @@ IFormFactor* FormFactorAnisoPyramid::sliceFormFactor(ZLimits limits, const IRota
                                                      kvector_t translation) const
 {
     auto effects = computeSlicingEffects(limits, translation, m_height);
-    double dbase_edge = 2 * effects.dz_bottom * MathFunctions::cot(m_alpha);
+    double dbase_edge = 2 * effects.dz_bottom * Math::cot(m_alpha);
     FormFactorAnisoPyramid slicedff(m_length - dbase_edge, m_width - dbase_edge,
                                     m_height - effects.dz_bottom - effects.dz_top, m_alpha);
     return createTransformedFormFactor(slicedff, rot, effects.position);
@@ -57,7 +57,7 @@ IFormFactor* FormFactorAnisoPyramid::sliceFormFactor(ZLimits limits, const IRota
 
 void FormFactorAnisoPyramid::onChange()
 {
-    double cot_alpha = MathFunctions::cot(m_alpha);
+    double cot_alpha = Math::cot(m_alpha);
     if (!std::isfinite(cot_alpha) || cot_alpha < 0)
         throw Exceptions::OutOfBoundsException("AnisoPyramid: angle alpha out of bounds");
     double r = cot_alpha * 2 * m_height / m_length;
diff --git a/Sample/HardParticle/FormFactorBox.cpp b/Sample/HardParticle/FormFactorBox.cpp
index 749060f42f4a3d140fc23ae3e9b1ad010df2a580..65faec6ffb3fd79d1d4f201e4873700806dc19fc 100644
--- a/Sample/HardParticle/FormFactorBox.cpp
+++ b/Sample/HardParticle/FormFactorBox.cpp
@@ -35,9 +35,8 @@ FormFactorBox::FormFactorBox(double length, double width, double height)
 complex_t FormFactorBox::evaluate_for_q(cvector_t q) const
 {
     complex_t qzHdiv2 = m_height / 2 * q.z();
-    return m_length * m_width * m_height * MathFunctions::sinc(m_length / 2 * q.x())
-           * MathFunctions::sinc(m_width / 2 * q.y()) * MathFunctions::sinc(qzHdiv2)
-           * exp_I(qzHdiv2);
+    return m_length * m_width * m_height * Math::sinc(m_length / 2 * q.x())
+           * Math::sinc(m_width / 2 * q.y()) * Math::sinc(qzHdiv2) * exp_I(qzHdiv2);
 }
 
 IFormFactor* FormFactorBox::sliceFormFactor(ZLimits limits, const IRotation& rot,
diff --git a/Sample/HardParticle/FormFactorCone.cpp b/Sample/HardParticle/FormFactorCone.cpp
index a938ad329f20de42ca64f6de1c3fab1045cd0a03..0495251272c917bb1d5bbf37b90725c106ba733e 100644
--- a/Sample/HardParticle/FormFactorCone.cpp
+++ b/Sample/HardParticle/FormFactorCone.cpp
@@ -15,6 +15,7 @@
 #include "Sample/HardParticle/FormFactorCone.h"
 #include "Base/Const/MathConstants.h"
 #include "Base/Types/Exceptions.h"
+#include "Base/Utils/Bessel.h"
 #include "Base/Utils/Integrator.h"
 #include "Base/Utils/MathFunctions.h"
 #include "Sample/Shapes/DoubleEllipse.h"
@@ -29,7 +30,7 @@ FormFactorCone::FormFactorCone(const std::vector<double> P)
                       P),
       m_radius(m_P[0]), m_height(m_P[1]), m_alpha(m_P[2])
 {
-    m_cot_alpha = MathFunctions::cot(m_alpha);
+    m_cot_alpha = Math::cot(m_alpha);
     if (!std::isfinite(m_cot_alpha) || m_cot_alpha < 0)
         throw Exceptions::OutOfBoundsException("pyramid angle alpha out of bounds");
     if (m_cot_alpha * m_height > m_radius) {
@@ -54,7 +55,7 @@ complex_t FormFactorCone::Integrand(double Z) const
 {
     double Rz = m_radius - Z * m_cot_alpha;
     complex_t q_p = std::sqrt(m_q.x() * m_q.x() + m_q.y() * m_q.y()); // sqrt(x*x + y*y)
-    return Rz * Rz * MathFunctions::Bessel_J1c(q_p * Rz) * exp_I(m_q.z() * Z);
+    return Rz * Rz * Math::Bessel::J1c(q_p * Rz) * exp_I(m_q.z() * Z);
 }
 
 complex_t FormFactorCone::evaluate_for_q(cvector_t q) const
@@ -87,7 +88,7 @@ IFormFactor* FormFactorCone::sliceFormFactor(ZLimits limits, const IRotation& ro
 
 void FormFactorCone::onChange()
 {
-    m_cot_alpha = MathFunctions::cot(m_alpha);
+    m_cot_alpha = Math::cot(m_alpha);
     double radius2 = m_radius - m_height * m_cot_alpha;
     m_shape = std::make_unique<DoubleEllipse>(m_radius, m_radius, m_height, radius2, radius2);
 }
diff --git a/Sample/HardParticle/FormFactorCone6.cpp b/Sample/HardParticle/FormFactorCone6.cpp
index 83e37e02d2cf24fb26be9f1ee6f745b4201008a9..85366b8ecd6d90587af437e3b2f88b57cf30aec6 100644
--- a/Sample/HardParticle/FormFactorCone6.cpp
+++ b/Sample/HardParticle/FormFactorCone6.cpp
@@ -49,7 +49,7 @@ IFormFactor* FormFactorCone6::sliceFormFactor(ZLimits limits, const IRotation& r
                                               kvector_t translation) const
 {
     auto effects = computeSlicingEffects(limits, translation, m_height);
-    double dbase_edge = effects.dz_bottom * MathFunctions::cot(m_alpha);
+    double dbase_edge = effects.dz_bottom * Math::cot(m_alpha);
     FormFactorCone6 slicedff(m_base_edge - dbase_edge,
                              m_height - effects.dz_bottom - effects.dz_top, m_alpha);
     return createTransformedFormFactor(slicedff, rot, effects.position);
@@ -57,7 +57,7 @@ IFormFactor* FormFactorCone6::sliceFormFactor(ZLimits limits, const IRotation& r
 
 void FormFactorCone6::onChange()
 {
-    double cot_alpha = MathFunctions::cot(m_alpha);
+    double cot_alpha = Math::cot(m_alpha);
     if (!std::isfinite(cot_alpha) || cot_alpha < 0)
         throw Exceptions::OutOfBoundsException("pyramid angle alpha out of bounds");
     double r = cot_alpha * 2 / sqrt(3) * m_height / m_base_edge; // L(top)/L(base)
diff --git a/Sample/HardParticle/FormFactorCuboctahedron.cpp b/Sample/HardParticle/FormFactorCuboctahedron.cpp
index 347bcc083c860be1ef8c9d1deda9084ce93516b3..1dbdca18aba3d7326ee02e4dd1830e4c0f26c91c 100644
--- a/Sample/HardParticle/FormFactorCuboctahedron.cpp
+++ b/Sample/HardParticle/FormFactorCuboctahedron.cpp
@@ -55,13 +55,13 @@ IFormFactor* FormFactorCuboctahedron::sliceFormFactor(ZLimits limits, const IRot
 {
     auto effects = computeSlicingEffects(limits, translation, m_height * (1 + m_height_ratio));
     if (effects.dz_bottom > m_height) {
-        double dbase_edge = 2 * (effects.dz_bottom - m_height) * MathFunctions::cot(m_alpha);
+        double dbase_edge = 2 * (effects.dz_bottom - m_height) * Math::cot(m_alpha);
         FormFactorPyramid slicedff(
             m_length - dbase_edge,
             m_height * (1 + m_height_ratio) - effects.dz_bottom - effects.dz_top, m_alpha);
         return createTransformedFormFactor(slicedff, rot, effects.position);
     } else if (effects.dz_top > m_height_ratio * m_height) {
-        double dbase_edge = 2 * (m_height - effects.dz_bottom) * MathFunctions::cot(m_alpha);
+        double dbase_edge = 2 * (m_height - effects.dz_bottom) * Math::cot(m_alpha);
         FormFactorPyramid slicedff(
             m_length - dbase_edge,
             m_height * (1 + m_height_ratio) - effects.dz_bottom - effects.dz_top, M_PI - m_alpha);
@@ -75,7 +75,7 @@ IFormFactor* FormFactorCuboctahedron::sliceFormFactor(ZLimits limits, const IRot
 
 void FormFactorCuboctahedron::onChange()
 {
-    double cot_alpha = MathFunctions::cot(m_alpha);
+    double cot_alpha = Math::cot(m_alpha);
     if (!std::isfinite(cot_alpha) || cot_alpha < 0)
         throw Exceptions::OutOfBoundsException("pyramid angle alpha out of bounds");
     double x = m_height_ratio;
diff --git a/Sample/HardParticle/FormFactorCylinder.cpp b/Sample/HardParticle/FormFactorCylinder.cpp
index 68acfd06720f1c8022f13c85bda1767cccd20349..4220f819ec42a83a74620a2c615f993ccbe1eb02 100644
--- a/Sample/HardParticle/FormFactorCylinder.cpp
+++ b/Sample/HardParticle/FormFactorCylinder.cpp
@@ -14,6 +14,7 @@
 
 #include "Sample/HardParticle/FormFactorCylinder.h"
 #include "Base/Const/MathConstants.h"
+#include "Base/Utils/Bessel.h"
 #include "Base/Utils/MathFunctions.h"
 #include "Sample/Shapes/DoubleEllipse.h"
 
@@ -39,9 +40,9 @@ complex_t FormFactorCylinder::evaluate_for_q(cvector_t q) const
     double H = m_height;
 
     complex_t qzH_half = q.z() * H / 2.0;
-    complex_t z_part = H * MathFunctions::sinc(qzH_half) * exp_I(qzH_half);
+    complex_t z_part = H * Math::sinc(qzH_half) * exp_I(qzH_half);
     complex_t qxy = std::sqrt(q.x() * q.x() + q.y() * q.y());
-    complex_t radial_part = M_TWOPI * R * R * MathFunctions::Bessel_J1c(qxy * R);
+    complex_t radial_part = M_TWOPI * R * R * Math::Bessel::J1c(qxy * R);
     complex_t result = radial_part * z_part;
 
     return result;
diff --git a/Sample/HardParticle/FormFactorEllipsoidalCylinder.cpp b/Sample/HardParticle/FormFactorEllipsoidalCylinder.cpp
index fbf4e354fdec35b27ea79f359f97a8f78fa2e6ba..4a10ec0b592aff8ef21cbb2fffd8b1178473f601 100644
--- a/Sample/HardParticle/FormFactorEllipsoidalCylinder.cpp
+++ b/Sample/HardParticle/FormFactorEllipsoidalCylinder.cpp
@@ -14,6 +14,7 @@
 
 #include "Sample/HardParticle/FormFactorEllipsoidalCylinder.h"
 #include "Base/Const/MathConstants.h"
+#include "Base/Utils/Bessel.h"
 #include "Base/Utils/MathFunctions.h"
 #include "Sample/Shapes/DoubleEllipse.h"
 
@@ -46,9 +47,9 @@ complex_t FormFactorEllipsoidalCylinder::evaluate_for_q(cvector_t q) const
     complex_t qyRb = q.y() * m_radius_y;
     complex_t qzHdiv2 = m_height / 2 * q.z();
 
-    complex_t Fz = exp_I(qzHdiv2) * MathFunctions::sinc(qzHdiv2);
+    complex_t Fz = exp_I(qzHdiv2) * Math::sinc(qzHdiv2);
     complex_t gamma = std::sqrt((qxRa) * (qxRa) + (qyRb) * (qyRb));
-    complex_t J1_gamma_div_gamma = MathFunctions::Bessel_J1c(gamma);
+    complex_t J1_gamma_div_gamma = Math::Bessel::J1c(gamma);
 
     return M_TWOPI * m_radius_x * m_radius_y * m_height * Fz * J1_gamma_div_gamma;
 }
diff --git a/Sample/HardParticle/FormFactorFullSphere.cpp b/Sample/HardParticle/FormFactorFullSphere.cpp
index 425906d335d48de002d68c25785c0ad87ab49bf0..11996457af40be63f59feb4a115bc90927a87e17 100644
--- a/Sample/HardParticle/FormFactorFullSphere.cpp
+++ b/Sample/HardParticle/FormFactorFullSphere.cpp
@@ -63,11 +63,10 @@ IFormFactor* FormFactorFullSphere::sliceFormFactor(ZLimits limits, const IRotati
     kvector_t rotation_offset =
         m_position_at_center ? kvector_t(0.0, 0.0, 0.0) : rot.transformed(center) - center;
     kvector_t new_translation = translation + rotation_offset;
-    std::unique_ptr<IRotation> P_identity(IRotation::createIdentity());
     double height = 2.0 * m_radius;
     auto effects = computeSlicingEffects(limits, new_translation, height);
     FormFactorTruncatedSphere slicedff(m_radius, height - effects.dz_bottom, effects.dz_top);
-    return createTransformedFormFactor(slicedff, *P_identity, effects.position);
+    return createTransformedFormFactor(slicedff, IdentityRotation(), effects.position);
 }
 
 void FormFactorFullSphere::onChange() {}
diff --git a/Sample/HardParticle/FormFactorHemiEllipsoid.cpp b/Sample/HardParticle/FormFactorHemiEllipsoid.cpp
index 4fc7f6337e1e13940bdf2677903b1385c2822b5c..4add3bed3683c0672468ba26780cfb0ecae03fb0 100644
--- a/Sample/HardParticle/FormFactorHemiEllipsoid.cpp
+++ b/Sample/HardParticle/FormFactorHemiEllipsoid.cpp
@@ -14,8 +14,8 @@
 
 #include "Sample/HardParticle/FormFactorHemiEllipsoid.h"
 #include "Base/Const/MathConstants.h"
+#include "Base/Utils/Bessel.h"
 #include "Base/Utils/Integrator.h"
-#include "Base/Utils/MathFunctions.h"
 #include "Sample/Shapes/TruncatedEllipsoid.h"
 #include <limits>
 
@@ -55,7 +55,7 @@ complex_t FormFactorHemiEllipsoid::Integrand(double Z) const
     complex_t qyWz = m_q.y() * Wz;
 
     complex_t gamma = std::sqrt(qxRz * qxRz + qyWz * qyWz);
-    complex_t J1_gamma_div_gamma = MathFunctions::Bessel_J1c(gamma);
+    complex_t J1_gamma_div_gamma = Math::Bessel::J1c(gamma);
 
     return Rz * Wz * J1_gamma_div_gamma * exp_I(m_q.z() * Z);
 }
diff --git a/Sample/HardParticle/FormFactorLongBoxGauss.cpp b/Sample/HardParticle/FormFactorLongBoxGauss.cpp
index c0e411680af00e482b14738d15c4905861c34349..00940fc5f55e4e1e22b3642cd1a83dfb0404ba66 100644
--- a/Sample/HardParticle/FormFactorLongBoxGauss.cpp
+++ b/Sample/HardParticle/FormFactorLongBoxGauss.cpp
@@ -39,8 +39,8 @@ complex_t FormFactorLongBoxGauss::evaluate_for_q(cvector_t q) const
     complex_t qyWdiv2 = m_width * q.y() / 2.0;
     complex_t qzHdiv2 = m_height * q.z() / 2.0;
 
-    return m_height * m_length * m_width * exp_I(qzHdiv2) * std::exp(-qxL2)
-           * MathFunctions::sinc(qyWdiv2) * MathFunctions::sinc(qzHdiv2);
+    return m_height * m_length * m_width * exp_I(qzHdiv2) * std::exp(-qxL2) * Math::sinc(qyWdiv2)
+           * Math::sinc(qzHdiv2);
 }
 
 IFormFactor* FormFactorLongBoxGauss::sliceFormFactor(ZLimits limits, const IRotation& rot,
diff --git a/Sample/HardParticle/FormFactorLongBoxLorentz.cpp b/Sample/HardParticle/FormFactorLongBoxLorentz.cpp
index ce8c005d1cc746e042e8c25382403a138b8d52d1..74a9676d00836ecea9d946fa6968ff96aee37b25 100644
--- a/Sample/HardParticle/FormFactorLongBoxLorentz.cpp
+++ b/Sample/HardParticle/FormFactorLongBoxLorentz.cpp
@@ -40,7 +40,7 @@ complex_t FormFactorLongBoxLorentz::evaluate_for_q(cvector_t q) const
     complex_t qzHdiv2 = m_height * q.z() / 2.0;
 
     return m_height * m_length * m_width * std::exp(complex_t(0., 1.) * qzHdiv2) / (1.0 + qxL2)
-           * MathFunctions::sinc(qyWdiv2) * MathFunctions::sinc(qzHdiv2);
+           * Math::sinc(qyWdiv2) * Math::sinc(qzHdiv2);
 }
 
 IFormFactor* FormFactorLongBoxLorentz::sliceFormFactor(ZLimits limits, const IRotation& rot,
diff --git a/Sample/HardParticle/FormFactorPyramid.cpp b/Sample/HardParticle/FormFactorPyramid.cpp
index bf7abf6e0d80c7832012139f6c12373db5ab0de0..863b6d569213abdb8d898459fc3a36a665fdbc75 100644
--- a/Sample/HardParticle/FormFactorPyramid.cpp
+++ b/Sample/HardParticle/FormFactorPyramid.cpp
@@ -48,7 +48,7 @@ IFormFactor* FormFactorPyramid::sliceFormFactor(ZLimits limits, const IRotation&
                                                 kvector_t translation) const
 {
     auto effects = computeSlicingEffects(limits, translation, m_height);
-    double dbase_edge = 2 * effects.dz_bottom * MathFunctions::cot(m_alpha);
+    double dbase_edge = 2 * effects.dz_bottom * Math::cot(m_alpha);
     FormFactorPyramid slicedff(m_base_edge - dbase_edge,
                                m_height - effects.dz_bottom - effects.dz_top, m_alpha);
     return createTransformedFormFactor(slicedff, rot, effects.position);
@@ -56,7 +56,7 @@ IFormFactor* FormFactorPyramid::sliceFormFactor(ZLimits limits, const IRotation&
 
 void FormFactorPyramid::onChange()
 {
-    double cot_alpha = MathFunctions::cot(m_alpha);
+    double cot_alpha = Math::cot(m_alpha);
     if (!std::isfinite(cot_alpha))
         throw Exceptions::OutOfBoundsException("pyramid angle alpha out of bounds");
     double r = cot_alpha * 2 * m_height / m_base_edge; // [L(base)-L(top)]/L(base)
diff --git a/Sample/HardParticle/FormFactorTetrahedron.cpp b/Sample/HardParticle/FormFactorTetrahedron.cpp
index cda831b1085f2acd75f0a40baa0ef9af486c1645..030f6c9b517e77375e7de40b765da2266f023583 100644
--- a/Sample/HardParticle/FormFactorTetrahedron.cpp
+++ b/Sample/HardParticle/FormFactorTetrahedron.cpp
@@ -46,7 +46,7 @@ IFormFactor* FormFactorTetrahedron::sliceFormFactor(ZLimits limits, const IRotat
                                                     kvector_t translation) const
 {
     auto effects = computeSlicingEffects(limits, translation, m_height);
-    double dbase_edge = 2 * sqrt(3) * effects.dz_bottom * MathFunctions::cot(m_alpha);
+    double dbase_edge = 2 * sqrt(3) * effects.dz_bottom * Math::cot(m_alpha);
     FormFactorTetrahedron slicedff(m_base_edge - dbase_edge,
                                    m_height - effects.dz_bottom - effects.dz_top, m_alpha);
     return createTransformedFormFactor(slicedff, rot, effects.position);
@@ -54,7 +54,7 @@ IFormFactor* FormFactorTetrahedron::sliceFormFactor(ZLimits limits, const IRotat
 
 void FormFactorTetrahedron::onChange()
 {
-    double cot_alpha = MathFunctions::cot(m_alpha);
+    double cot_alpha = Math::cot(m_alpha);
     if (!std::isfinite(cot_alpha) || cot_alpha < 0)
         throw Exceptions::OutOfBoundsException("pyramid angle alpha out of bounds");
     double r = cot_alpha * 2 * std::sqrt(3.) * m_height / m_base_edge; // L(top)/L(base)
diff --git a/Sample/HardParticle/FormFactorTruncatedSphere.cpp b/Sample/HardParticle/FormFactorTruncatedSphere.cpp
index 7191783b624c949957cc887787d163a884a0a34e..ebd59fd86244ca6a0094559fc644ef36a021e806 100644
--- a/Sample/HardParticle/FormFactorTruncatedSphere.cpp
+++ b/Sample/HardParticle/FormFactorTruncatedSphere.cpp
@@ -15,8 +15,8 @@
 #include "Sample/HardParticle/FormFactorTruncatedSphere.h"
 #include "Base/Const/MathConstants.h"
 #include "Base/Types/Exceptions.h"
+#include "Base/Utils/Bessel.h"
 #include "Base/Utils/Integrator.h"
-#include "Base/Utils/MathFunctions.h"
 #include "Fit/Tools/RealLimits.h"
 #include "Sample/Shapes/TruncatedEllipsoid.h"
 #include <limits>
@@ -60,7 +60,7 @@ complex_t FormFactorTruncatedSphere::Integrand(double Z) const
     complex_t qx = m_q.x();
     complex_t qy = m_q.y();
     complex_t q_p = std::sqrt(qx * qx + qy * qy); // NOT the modulus!
-    return Rz * Rz * MathFunctions::Bessel_J1c(q_p * Rz) * exp_I(m_q.z() * Z);
+    return Rz * Rz * Math::Bessel::J1c(q_p * Rz) * exp_I(m_q.z() * Z);
 }
 
 //! Complex form factor.
diff --git a/Sample/HardParticle/FormFactorTruncatedSpheroid.cpp b/Sample/HardParticle/FormFactorTruncatedSpheroid.cpp
index 9e0606ad3c85234a79bf11c3e4085561dbb3a8a1..ac5d0775a38fe1a673b9b752c8b524f22d624bff 100644
--- a/Sample/HardParticle/FormFactorTruncatedSpheroid.cpp
+++ b/Sample/HardParticle/FormFactorTruncatedSpheroid.cpp
@@ -15,8 +15,8 @@
 #include "Sample/HardParticle/FormFactorTruncatedSpheroid.h"
 #include "Base/Const/MathConstants.h"
 #include "Base/Types/Exceptions.h"
+#include "Base/Utils/Bessel.h"
 #include "Base/Utils/Integrator.h"
-#include "Base/Utils/MathFunctions.h"
 #include "Sample/Shapes/TruncatedEllipsoid.h"
 #include <limits>
 
@@ -64,7 +64,7 @@ complex_t FormFactorTruncatedSpheroid::Integrand(double Z) const
 
     double Rz = std::sqrt(R * R - Z * Z / (fp * fp));
     complex_t qrRz = std::sqrt(m_q.x() * m_q.x() + m_q.y() * m_q.y()) * Rz;
-    complex_t J1_qrRz_div_qrRz = MathFunctions::Bessel_J1c(qrRz);
+    complex_t J1_qrRz_div_qrRz = Math::Bessel::J1c(qrRz);
 
     return Rz * Rz * J1_qrRz_div_qrRz * exp_I(m_q.z() * Z);
 }
diff --git a/Sample/HardParticle/PolyhedralComponents.cpp b/Sample/HardParticle/PolyhedralComponents.cpp
index 77c6e3c7ee4277aa61a32962dfe9ab982875bd1d..e6265a6bb9e8cf5039a1ba2a28bf7317952f6d51 100644
--- a/Sample/HardParticle/PolyhedralComponents.cpp
+++ b/Sample/HardParticle/PolyhedralComponents.cpp
@@ -310,14 +310,13 @@ complex_t PolyhedralFace::edge_sum_ff(cvector_t q, cvector_t qpa, bool sym_Ci) c
         } else {
             vfac = -vfacsum; // to improve numeric accuracy: qcE_J = - sum_{j=0}^{J-1} qcE_j
         }
-        complex_t term = vfac * MathFunctions::sinc(qE) * Rfac;
+        complex_t term = vfac * Math::sinc(qE) * Rfac;
         sum += term;
 #ifdef POLYHEDRAL_DIAGNOSTIC
         if (diagnosis.debmsg >= 2)
             std::cout << std::scientific << std::showpos << std::setprecision(16)
                       << "    sum=" << sum << " term=" << term << " vf=" << vfac << " qE=" << qE
-                      << " qR=" << qR << " sinc=" << MathFunctions::sinc(qE) << " Rfac=" << Rfac
-                      << "\n";
+                      << " qR=" << qR << " sinc=" << Math::sinc(qE) << " Rfac=" << Rfac << "\n";
 #endif
     }
     return sum;
diff --git a/Sample/HardParticle/Prism.cpp b/Sample/HardParticle/Prism.cpp
index 48f3a072e4895498f5d0b0cea44a46f56672ab30..8cd96e507577a4181755588bd73199d109a47264 100644
--- a/Sample/HardParticle/Prism.cpp
+++ b/Sample/HardParticle/Prism.cpp
@@ -60,7 +60,7 @@ complex_t Prism::evaluate_for_q(const cvector_t& q) const
         diagnosis.nExpandedFaces = 0;
 #endif
         cvector_t qxy(q.x(), q.y(), 0.);
-        return m_height * exp_I(m_height / 2 * q.z()) * MathFunctions::sinc(m_height / 2 * q.z())
+        return m_height * exp_I(m_height / 2 * q.z()) * Math::sinc(m_height / 2 * q.z())
                * m_base->ff_2D(qxy);
     } catch (std::logic_error& e) {
         throw std::logic_error(std::string("Bug in Prism: ") + e.what()
diff --git a/Sample/HardParticle/Ripples.cpp b/Sample/HardParticle/Ripples.cpp
index 011ecac1550e810a904b5825bcc724853bae885a..5b9dd4d836e8fe2c39c0e73a56bf796bc3cf0427 100644
--- a/Sample/HardParticle/Ripples.cpp
+++ b/Sample/HardParticle/Ripples.cpp
@@ -19,7 +19,7 @@
 
 complex_t ripples::factor_x_box(complex_t q, double r)
 {
-    return r * MathFunctions::sinc(q * r / 2.0);
+    return r * Math::sinc(q * r / 2.0);
 }
 
 complex_t ripples::factor_x_Gauss(complex_t q, double r)
@@ -38,8 +38,7 @@ complex_t ripples::profile_yz_bar(complex_t qy, complex_t qz, double width, doub
     const complex_t qyWdiv2 = width * qy / 2.0;
     const complex_t qzHdiv2 = height * qz / 2.0;
 
-    return height * width * exp_I(qzHdiv2) * MathFunctions::sinc(qyWdiv2)
-           * MathFunctions::sinc(qzHdiv2);
+    return height * width * exp_I(qzHdiv2) * Math::sinc(qyWdiv2) * Math::sinc(qzHdiv2);
 }
 
 //! Complex form factor of triangular ripple.
@@ -55,7 +54,7 @@ complex_t ripples::profile_yz_cosine(complex_t qy, complex_t qz, double width, d
         complex_t aaa2 = aaa * aaa;
         if (aaa2 == 1.)
             return factor * M_PI_4 * height;
-        return factor * M_PI_2 * height * MathFunctions::sinc(qy * width * 0.5) / (1.0 - aaa2);
+        return factor * M_PI_2 * height * Math::sinc(qy * width * 0.5) / (1.0 - aaa2);
     }
 
     // numerical integration otherwise
@@ -94,8 +93,7 @@ complex_t ripples::profile_yz_triangular(complex_t qy, complex_t qz, double widt
     } else {
         const complex_t gamma_p = (a + qyW2) * 0.5;
         const complex_t gamma_m = (a - qyW2) * 0.5;
-        result = exp_I(gamma_m) * MathFunctions::sinc(gamma_p)
-                 - exp_I(gamma_p) * MathFunctions::sinc(gamma_m);
+        result = exp_I(gamma_m) * Math::sinc(gamma_p) - exp_I(gamma_p) * Math::sinc(gamma_m);
         result = mul_I(exp_I(-qyd) * result / (qyW2 * 2.));
     }
     return factor * result;
diff --git a/Sample/Particle/IParticle.cpp b/Sample/Particle/IParticle.cpp
index 58bb692971dd59e7c5b1f155effec18b8d221ed6..cac3fcf579543f0eb0ab5ff9d3d49001394c18ea 100644
--- a/Sample/Particle/IParticle.cpp
+++ b/Sample/Particle/IParticle.cpp
@@ -94,7 +94,7 @@ SafePointerVector<IParticle> IParticle::decompose() const
 ParticleLimits IParticle::bottomTopZ() const
 {
     std::unique_ptr<IFormFactor> P_ff(createFormFactor());
-    std::unique_ptr<IRotation> P_rot(IRotation::createIdentity());
+    std::unique_ptr<IRotation> P_rot(new IdentityRotation);
     return {P_ff->bottomZ(*P_rot), P_ff->topZ(*P_rot)};
 }
 
diff --git a/Sample/Particle/MesoCrystal.cpp b/Sample/Particle/MesoCrystal.cpp
index 4ddb436ef968f1a87b16e1eb5181f60192eb75e3..ead55e97686eab54757ddfc15bf6f62702b9b8c0 100644
--- a/Sample/Particle/MesoCrystal.cpp
+++ b/Sample/Particle/MesoCrystal.cpp
@@ -47,7 +47,7 @@ SlicedParticle MesoCrystal::createSlicedParticle(ZLimits limits) const
 {
     if (!m_particle_structure || !m_meso_form_factor)
         return {};
-    std::unique_ptr<IRotation> P_rotation(IRotation::createIdentity());
+    std::unique_ptr<IRotation> P_rotation(new IdentityRotation);
     if (m_rotation)
         P_rotation.reset(m_rotation->clone());
     std::unique_ptr<IFormFactor> P_tem_ff(
diff --git a/Sample/Particle/Particle.cpp b/Sample/Particle/Particle.cpp
index 487299279a713ae447eab64550a8cda0555f43a5..8510998e5cf10f66c78481c841158ea3b0574389 100644
--- a/Sample/Particle/Particle.cpp
+++ b/Sample/Particle/Particle.cpp
@@ -59,7 +59,7 @@ SlicedParticle Particle::createSlicedParticle(ZLimits limits) const
 {
     if (!m_form_factor)
         return {};
-    std::unique_ptr<IRotation> P_rotation(IRotation::createIdentity());
+    std::unique_ptr<IRotation> P_rotation(new IdentityRotation);
     if (m_rotation)
         P_rotation.reset(m_rotation->clone());
     std::unique_ptr<IFormFactor> P_tem_ff(
diff --git a/Sample/Particle/ParticleCoreShell.cpp b/Sample/Particle/ParticleCoreShell.cpp
index 051b092d902578c7894ebd8fdac5679a79a33fc5..98f1bbe46260bd0b3e593733eec0a0622d693ed4 100644
--- a/Sample/Particle/ParticleCoreShell.cpp
+++ b/Sample/Particle/ParticleCoreShell.cpp
@@ -42,7 +42,7 @@ SlicedParticle ParticleCoreShell::createSlicedParticle(ZLimits limits) const
 {
     if (!m_core || !m_shell)
         return {};
-    std::unique_ptr<IRotation> P_rotation(IRotation::createIdentity());
+    std::unique_ptr<IRotation> P_rotation(new IdentityRotation);
     if (m_rotation)
         P_rotation.reset(m_rotation->clone());
 
diff --git a/Sample/SampleBuilderEngine/IRegistry.h b/Sample/SampleBuilderEngine/IRegistry.h
index 8660aff3179bedba7a495ae7eb78c1341e3e5e33..0d33299a298769eddab85fab3bf5c32defc9e64f 100644
--- a/Sample/SampleBuilderEngine/IRegistry.h
+++ b/Sample/SampleBuilderEngine/IRegistry.h
@@ -15,11 +15,11 @@
 #ifndef BORNAGAIN_SAMPLE_SAMPLEBUILDERENGINE_IREGISTRY_H
 #define BORNAGAIN_SAMPLE_SAMPLEBUILDERENGINE_IREGISTRY_H
 
-#include "Base/Types/Exceptions.h"
 #include <map>
 #include <memory>
 #include <string>
 #include <vector>
+#include <stdexcept>
 
 //! @class IRegistry
 //! @ingroup tools_internal
@@ -32,8 +32,7 @@ public:
     {
         auto it = m_data.find(key);
         if (it == m_data.end())
-            throw Exceptions::UnknownClassRegistrationException(
-                "IRegistry::createItem() -> Error. Not existing item key '" + key + "'");
+            throw std::runtime_error("Key '" + key + "' not found in registry");
         return it->second.get();
     }
 
@@ -51,8 +50,7 @@ protected:
     void add(const std::string& key, ValueType* item)
     {
         if (m_data.find(key) != m_data.end())
-            throw Exceptions::ExistingClassRegistrationException(
-                "IRegistry::createItem() -> Error. Already existing item with key '" + key + "'");
+            throw std::runtime_error("Key '" + key + "' already in registry");
         m_data[key] = std::unique_ptr<ValueType>(item);
     }
 
diff --git a/Sample/Scattering/Rotations.cpp b/Sample/Scattering/Rotations.cpp
index cf9836c6317fa192e704dff8aaddb96a05a94c57..448421f27a474bb66d39d616327104f5ec0bd4a9 100644
--- a/Sample/Scattering/Rotations.cpp
+++ b/Sample/Scattering/Rotations.cpp
@@ -50,11 +50,6 @@ IRotation* IRotation::createRotation(const Transform3D& transform)
     ASSERT(0); // impossible case
 }
 
-IRotation* IRotation::createIdentity()
-{
-    return new RotationZ(0.0);
-}
-
 kvector_t IRotation::transformed(const kvector_t& v) const
 {
     return getTransform3D().transformed(v);
@@ -91,7 +86,7 @@ IdentityRotation::IdentityRotation()
 
 Transform3D IdentityRotation::getTransform3D() const
 {
-    return Transform3D::createIdentity();
+    return {};
 }
 
 // ************************************************************************** //
diff --git a/Sample/Scattering/Rotations.h b/Sample/Scattering/Rotations.h
index 1b4b94a6737a91657f61c58d78a7d1527982a5d9..d9babbe11ecc17e124e8876aaf2b33aa095b35db 100644
--- a/Sample/Scattering/Rotations.h
+++ b/Sample/Scattering/Rotations.h
@@ -27,7 +27,6 @@ class IRotation : public ICloneable, public INode
 {
 public:
     static IRotation* createRotation(const Transform3D& transform);
-    static IRotation* createIdentity();
 
     IRotation(const NodeMeta& meta, const std::vector<double>& PValues);
 
@@ -51,7 +50,7 @@ IRotation* createProduct(const IRotation& left, const IRotation& right);
 
 //! The identity rotation, which leaves everything in place.
 
-class IdentityRotation : public IRotation // TODO get rid of this class
+class IdentityRotation : public IRotation // TODO RECONSIDER: merge with class IRotation
 {
 public:
     IdentityRotation();
diff --git a/Sample/Specular/SpecularMagneticNewTanhStrategy.cpp b/Sample/Specular/SpecularMagneticNewTanhStrategy.cpp
index 1b6677960c85aff80248e9a5c0dd3cb1d1d36ab2..7bf2b77c5f78491c26060a0873d7d15071a6d5ad 100644
--- a/Sample/Specular/SpecularMagneticNewTanhStrategy.cpp
+++ b/Sample/Specular/SpecularMagneticNewTanhStrategy.cpp
@@ -36,8 +36,8 @@ SpecularMagneticNewTanhStrategy::computeRoughnessMatrix(const MatrixRTCoefficien
         const double factor1 = 2. * (1. + b.z());
         Q << (1. + b.z()), (I * b.y() - b.x()), (b.x() + I * b.y()), (b.z() + 1.);
 
-        complex_t l1 = std::sqrt(MathFunctions::tanhc(sigeff * coeff.m_lambda(1)));
-        complex_t l2 = std::sqrt(MathFunctions::tanhc(sigeff * coeff.m_lambda(0)));
+        complex_t l1 = std::sqrt(Math::tanhc(sigeff * coeff.m_lambda(1)));
+        complex_t l2 = std::sqrt(Math::tanhc(sigeff * coeff.m_lambda(0)));
 
         if (inverse) {
             l1 = 1. / l1;
@@ -50,7 +50,7 @@ SpecularMagneticNewTanhStrategy::computeRoughnessMatrix(const MatrixRTCoefficien
 
     } else if (b.mag() < 10 * std::numeric_limits<double>::epsilon()) {
         complex_t alpha =
-            std::sqrt(MathFunctions::tanhc(0.5 * sigeff * (coeff.m_lambda(1) + coeff.m_lambda(0))));
+            std::sqrt(Math::tanhc(0.5 * sigeff * (coeff.m_lambda(1) + coeff.m_lambda(0))));
         if (inverse)
             alpha = 1. / alpha;
         const Eigen::Matrix2cd lambda = Eigen::DiagonalMatrix<complex_t, 2>({alpha, alpha});
diff --git a/Sample/Specular/SpecularScalarTanhStrategy.cpp b/Sample/Specular/SpecularScalarTanhStrategy.cpp
index 85d3229304612c124ec7bde141d3db0b066da4a1..5c13a941e79b2d5cf4107c3494d0e605b14d672f 100644
--- a/Sample/Specular/SpecularScalarTanhStrategy.cpp
+++ b/Sample/Specular/SpecularScalarTanhStrategy.cpp
@@ -28,8 +28,7 @@ SpecularScalarTanhStrategy::transition(complex_t kzi, complex_t kzi1, double sig
     complex_t roughness = 1;
     if (sigma > 0.0) {
         const double sigeff = pi2_15 * sigma;
-        roughness =
-            std::sqrt(MathFunctions::tanhc(sigeff * kzi1) / MathFunctions::tanhc(sigeff * kzi));
+        roughness = std::sqrt(Math::tanhc(sigeff * kzi1) / Math::tanhc(sigeff * kzi));
     }
     const complex_t inv_roughness = 1.0 / roughness;
     const complex_t kz_ratio = kzi1 / kzi * roughness;
diff --git a/Tests/UnitTests/Core/Detector/BesselTest.cpp b/Tests/UnitTests/Core/Detector/BesselTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..bc787064c77576272c31f5dfc7fc8f9e58829dff
--- /dev/null
+++ b/Tests/UnitTests/Core/Detector/BesselTest.cpp
@@ -0,0 +1,37 @@
+#include "Base/Utils/Bessel.h"
+#include "Tests/GTestWrapper/google_test.h"
+
+class BesselTest : public ::testing::Test
+{
+};
+
+// Test complex Bessel function J1
+TEST_F(BesselTest, BesselJ1)
+{
+    const double eps = 4.7e-16; // more than twice the machine precision
+    complex_t res;
+
+    // Test four arbitrary function values.
+    // Reference values are computed using dev-tools/math/cbesselJ01.c,
+    //   from which BesselJ1 has been derived. So this is _not_ an independent
+    //   for numeric accuracy, but only for invariance.
+    // However, the four specific results below have also been checked against the
+    //   online calculator http://keisan.casio.com/exec/system/1180573474.
+    //   Agreement is excellent, except for the dominantly imaginary argument .01+100i.
+    //   In conclusion, Bessel::J1 is clearly good enough for our purpose.
+    res = Math::Bessel::J1(complex_t(0.8, 1.5));
+    EXPECT_NEAR(res.real(), 0.72837687825769404, eps * std::abs(res)); // Keisan ..69398...
+    EXPECT_NEAR(res.imag(), 0.75030568686427268, eps * std::abs(res)); // Keisan ..27264...
+
+    res = Math::Bessel::J1(complex_t(1e-2, 1e2));
+    EXPECT_NEAR(res.real(), 1.0630504683139779e+40, eps * std::abs(res)); // Keisan 1.063015...
+    EXPECT_NEAR(res.imag(), 1.0683164984973165e+42, eps * std::abs(res)); // Keisan ..73162...
+
+    res = Math::Bessel::J1(complex_t(-1e2, 1e-2));
+    EXPECT_NEAR(res.real(), 0.077149198549289394, eps * std::abs(res)); // Keisan ..89252...
+    EXPECT_NEAR(res.imag(), 2.075766253119904e-4, eps * std::abs(res)); // Keisan ..19951...
+
+    res = Math::Bessel::J1(complex_t(7, 9));
+    EXPECT_NEAR(res.real(), 370.00180888861155, eps * std::abs(res)); // Keisan ..61107...
+    EXPECT_NEAR(res.imag(), 856.00300811057934, eps * std::abs(res)); // Keisan ..57940...
+}
diff --git a/Tests/UnitTests/Core/Detector/SpecialFunctionsTest.cpp b/Tests/UnitTests/Core/Detector/SpecialFunctionsTest.cpp
index 80fb7aa08733867e29c1c3e14d0add5c86099ddd..e2406da58a8e3779c8470f7f267bdc53f42998f4 100644
--- a/Tests/UnitTests/Core/Detector/SpecialFunctionsTest.cpp
+++ b/Tests/UnitTests/Core/Detector/SpecialFunctionsTest.cpp
@@ -10,37 +10,6 @@ class SpecialFunctionsTest : public ::testing::Test
 {
 };
 
-// Test complex Bessel function J1
-TEST_F(SpecialFunctionsTest, BesselJ1)
-{
-    const double eps = 4.7e-16; // more than twice the machine precision
-    complex_t res;
-
-    // Test four arbitrary function values.
-    // Reference values are computed using dev-tools/math/cbesselJ01.c,
-    //   from which BesselJ1 has been derived. So this is _not_ an independent
-    //   for numeric accuracy, but only for invariance.
-    // However, the four specific results below have also been checked against the
-    //   online calculator http://keisan.casio.com/exec/system/1180573474.
-    //   Agreement is excellent, except for the dominantly imaginary argument .01+100i.
-    //   In conclusion, Bessel_J1 is clearly good enough for our purpose.
-    res = MathFunctions::Bessel_J1(complex_t(0.8, 1.5));
-    EXPECT_NEAR(res.real(), 0.72837687825769404, eps * std::abs(res)); // Keisan ..69398...
-    EXPECT_NEAR(res.imag(), 0.75030568686427268, eps * std::abs(res)); // Keisan ..27264...
-
-    res = MathFunctions::Bessel_J1(complex_t(1e-2, 1e2));
-    EXPECT_NEAR(res.real(), 1.0630504683139779e+40, eps * std::abs(res)); // Keisan 1.063015...
-    EXPECT_NEAR(res.imag(), 1.0683164984973165e+42, eps * std::abs(res)); // Keisan ..73162...
-
-    res = MathFunctions::Bessel_J1(complex_t(-1e2, 1e-2));
-    EXPECT_NEAR(res.real(), 0.077149198549289394, eps * std::abs(res)); // Keisan ..89252...
-    EXPECT_NEAR(res.imag(), 2.075766253119904e-4, eps * std::abs(res)); // Keisan ..19951...
-
-    res = MathFunctions::Bessel_J1(complex_t(7, 9));
-    EXPECT_NEAR(res.real(), 370.00180888861155, eps * std::abs(res)); // Keisan ..61107...
-    EXPECT_NEAR(res.imag(), 856.00300811057934, eps * std::abs(res)); // Keisan ..57940...
-}
-
 // Test accuracy of complex function sinc(z) near the removable singularity at z=0
 
 TEST_F(SpecialFunctionsTest, csinc)
@@ -51,47 +20,47 @@ TEST_F(SpecialFunctionsTest, csinc)
         double ph = M_TWOPI * i / 24;
         // std::cout << "---------------------------------------------------------------------\n";
         // std::cout << "phase = " << ph << "\n";
-        EXPECT_EQ(MathFunctions::sinc(complex_t(0, 0)), complex_t(1., 0.));
+        EXPECT_EQ(Math::sinc(complex_t(0, 0)), complex_t(1., 0.));
         complex_t z;
         z = std::polar(1e-17, ph);
-        EXPECT_CNEAR(MathFunctions::sinc(z), complex_t(1., 0.), eps);
+        EXPECT_CNEAR(Math::sinc(z), complex_t(1., 0.), eps);
         z = std::polar(2e-17, ph);
-        EXPECT_CNEAR(MathFunctions::sinc(z), complex_t(1., 0.), eps);
+        EXPECT_CNEAR(Math::sinc(z), complex_t(1., 0.), eps);
         z = std::polar(5e-17, ph);
-        EXPECT_CNEAR(MathFunctions::sinc(z), complex_t(1., 0.), eps);
+        EXPECT_CNEAR(Math::sinc(z), complex_t(1., 0.), eps);
         z = std::polar(1e-16, ph);
-        EXPECT_CNEAR(MathFunctions::sinc(z), complex_t(1., 0.), eps);
+        EXPECT_CNEAR(Math::sinc(z), complex_t(1., 0.), eps);
         z = std::polar(2e-16, ph);
-        EXPECT_CNEAR(MathFunctions::sinc(z), complex_t(1., 0.), eps);
+        EXPECT_CNEAR(Math::sinc(z), complex_t(1., 0.), eps);
         z = std::polar(5e-16, ph);
-        EXPECT_CNEAR(MathFunctions::sinc(z), complex_t(1., 0.), eps);
+        EXPECT_CNEAR(Math::sinc(z), complex_t(1., 0.), eps);
         z = std::polar(1e-15, ph);
-        EXPECT_CNEAR(MathFunctions::sinc(z), complex_t(1., 0.), eps);
+        EXPECT_CNEAR(Math::sinc(z), complex_t(1., 0.), eps);
         z = std::polar(1e-13, ph);
-        EXPECT_CNEAR(MathFunctions::sinc(z), complex_t(1., 0.), eps);
+        EXPECT_CNEAR(Math::sinc(z), complex_t(1., 0.), eps);
         z = std::polar(1e-11, ph);
-        EXPECT_CNEAR(MathFunctions::sinc(z), complex_t(1., 0.), eps);
+        EXPECT_CNEAR(Math::sinc(z), complex_t(1., 0.), eps);
         z = std::polar(1e-9, ph);
-        EXPECT_CNEAR(MathFunctions::sinc(z), complex_t(1., 0.), eps);
+        EXPECT_CNEAR(Math::sinc(z), complex_t(1., 0.), eps);
         z = std::polar(5e-8, ph);
-        EXPECT_CNEAR(MathFunctions::sinc(z), 1. - z * z / 6., eps);
+        EXPECT_CNEAR(Math::sinc(z), 1. - z * z / 6., eps);
         z = std::polar(2e-8, ph);
-        EXPECT_CNEAR(MathFunctions::sinc(z), 1. - z * z / 6., eps);
+        EXPECT_CNEAR(Math::sinc(z), 1. - z * z / 6., eps);
         z = std::polar(1e-8, ph);
-        EXPECT_CNEAR(MathFunctions::sinc(z), 1. - z * z / 6., eps);
+        EXPECT_CNEAR(Math::sinc(z), 1. - z * z / 6., eps);
         z = std::polar(5e-7, ph);
-        EXPECT_CNEAR(MathFunctions::sinc(z), 1. - z * z / 6., eps);
+        EXPECT_CNEAR(Math::sinc(z), 1. - z * z / 6., eps);
         z = std::polar(2e-7, ph);
-        EXPECT_CNEAR(MathFunctions::sinc(z), 1. - z * z / 6., eps);
+        EXPECT_CNEAR(Math::sinc(z), 1. - z * z / 6., eps);
         z = std::polar(1e-7, ph);
-        EXPECT_CNEAR(MathFunctions::sinc(z), 1. - z * z / 6., eps);
+        EXPECT_CNEAR(Math::sinc(z), 1. - z * z / 6., eps);
         z = std::polar(1e-6, ph);
-        EXPECT_CNEAR(MathFunctions::sinc(z), 1. - z * z / 6., eps);
+        EXPECT_CNEAR(Math::sinc(z), 1. - z * z / 6., eps);
         z = std::polar(1e-5, ph);
-        EXPECT_CNEAR(MathFunctions::sinc(z), 1. - z * z / 6., eps);
+        EXPECT_CNEAR(Math::sinc(z), 1. - z * z / 6., eps);
         z = std::polar(1e-4, ph);
-        EXPECT_CNEAR(MathFunctions::sinc(z), 1. - z * z / 6. * (1. - z * z / 20.), eps);
+        EXPECT_CNEAR(Math::sinc(z), 1. - z * z / 6. * (1. - z * z / 20.), eps);
         z = std::polar(1e-3, ph);
-        EXPECT_CNEAR(MathFunctions::sinc(z), 1. - z * z / 6. * (1. - z * z / 20.), eps);
+        EXPECT_CNEAR(Math::sinc(z), 1. - z * z / 6. * (1. - z * z / 20.), eps);
     }
 }
diff --git a/auto/Wrap/doxygenBase.i b/auto/Wrap/doxygenBase.i
index 04d4c0a481982acf192f85d509d43f123b5dc7a1..65b76c17333b35efdebfe75d29b0c6f9dca7249d 100644
--- a/auto/Wrap/doxygenBase.i
+++ b/auto/Wrap/doxygenBase.i
@@ -376,13 +376,6 @@ Creates a new clipped axis.
 ";
 
 
-// File: classExceptions_1_1DivisionByZeroException.xml
-%feature("docstring") Exceptions::DivisionByZeroException "";
-
-%feature("docstring")  Exceptions::DivisionByZeroException::DivisionByZeroException "Exceptions::DivisionByZeroException::DivisionByZeroException(const std::string &message)
-";
-
-
 // File: classExceptions_1_1DomainErrorException.xml
 %feature("docstring") Exceptions::DomainErrorException "";
 
@@ -390,13 +383,6 @@ Creates a new clipped axis.
 ";
 
 
-// File: classExceptions_1_1ExistingClassRegistrationException.xml
-%feature("docstring") Exceptions::ExistingClassRegistrationException "";
-
-%feature("docstring")  Exceptions::ExistingClassRegistrationException::ExistingClassRegistrationException "Exceptions::ExistingClassRegistrationException::ExistingClassRegistrationException(const std::string &message)
-";
-
-
 // File: classExceptions_1_1FileIsBadException.xml
 %feature("docstring") Exceptions::FileIsBadException "";
 
@@ -1070,13 +1056,6 @@ Determine if the transformation is trivial (identity)
 ";
 
 
-// File: classExceptions_1_1UnknownClassRegistrationException.xml
-%feature("docstring") Exceptions::UnknownClassRegistrationException "";
-
-%feature("docstring")  Exceptions::UnknownClassRegistrationException::UnknownClassRegistrationException "Exceptions::UnknownClassRegistrationException::UnknownClassRegistrationException(const std::string &message)
-";
-
-
 // File: classVariableBinAxis.xml
 %feature("docstring") VariableBinAxis "
 
@@ -1153,6 +1132,9 @@ Creates a new clipped axis.
 ";
 
 
+// File: namespace_0d30.xml
+
+
 // File: namespacealgo.xml
 %feature("docstring")  algo::almostEqual "bool algo::almostEqual(double a, double b)
 
@@ -1242,129 +1224,84 @@ Returns true if file with given name exists on disk.
 ";
 
 
-// File: namespaceMathFunctions.xml
-%feature("docstring")  MathFunctions::Bessel_J0_PowSer "complex_t MathFunctions::Bessel_J0_PowSer(const complex_t z)
-
-Computes complex Bessel function J0(z), using power series and asymptotic expansion.
-
-Computes the complex Bessel function J0(z), using power series and asymptotic expansion.
-
-Forked from unoptimized code at http://www.crbond.com/math.htm, who refers to \"Computation of Special Functions\", Zhang and Jin, John Wiley and Sons, 1996. 
+// File: namespaceMath.xml
+%feature("docstring")  Math::Bessel::StandardNormal "double Math::StandardNormal(double x)
 ";
 
-%feature("docstring")  MathFunctions::Bessel_J1_PowSer "complex_t MathFunctions::Bessel_J1_PowSer(const complex_t z)
-
-Computes complex Bessel function J0(z), using power series and asymptotic expansion.
-
-Computes the complex Bessel function J1(z), using power series and asymptotic expansion.
-
-Forked from same source as for Bessel_J0_PowSer 
+%feature("docstring")  Math::Bessel::Gaussian "double Math::Gaussian(double x, double average, double std_dev)
 ";
 
-%feature("docstring")  MathFunctions::StandardNormal "double MathFunctions::StandardNormal(double x)
+%feature("docstring")  Math::Bessel::IntegratedGaussian "double Math::IntegratedGaussian(double x, double average, double std_dev)
 ";
 
-%feature("docstring")  MathFunctions::Gaussian "double MathFunctions::Gaussian(double x, double average, double std_dev)
-";
-
-%feature("docstring")  MathFunctions::IntegratedGaussian "double MathFunctions::IntegratedGaussian(double x, double average, double std_dev)
-";
-
-%feature("docstring")  MathFunctions::cot "double MathFunctions::cot(double x)
+%feature("docstring")  Math::Bessel::cot "double Math::cot(double x)
 
 cotangent function:  $cot(x)\\\\equiv1/tan(x)$
 ";
 
-%feature("docstring")  MathFunctions::Si "double MathFunctions::Si(double x)
-
-Sine integral function:  $Si(x)\\\\equiv\\\\int_0^x du \\\\sin(u)/u$. 
-";
-
-%feature("docstring")  MathFunctions::sinc "double MathFunctions::sinc(double x)
+%feature("docstring")  Math::Bessel::sinc "double Math::sinc(double x)
 
 sinc function:  $sinc(x)\\\\equiv\\\\sin(x)/x$
 ";
 
-%feature("docstring")  MathFunctions::sinc "complex_t MathFunctions::sinc(const complex_t z)
+%feature("docstring")  Math::Bessel::sinc "complex_t Math::sinc(const complex_t z)
 
 Complex sinc function:  $sinc(x)\\\\equiv\\\\sin(x)/x$. 
 ";
 
-%feature("docstring")  MathFunctions::tanhc "complex_t MathFunctions::tanhc(const complex_t z)
+%feature("docstring")  Math::Bessel::tanhc "complex_t Math::tanhc(const complex_t z)
 
 Complex tanhc function:  $tanhc(x)\\\\equiv\\\\tanh(x)/x$. 
 ";
 
-%feature("docstring")  MathFunctions::Laue "double MathFunctions::Laue(const double x, size_t N)
+%feature("docstring")  Math::Bessel::Laue "double Math::Laue(const double x, size_t N)
 
 Real Laue function:  $Laue(x,N)\\\\equiv\\\\sin(Nx)/sin(x)$. 
 ";
 
-%feature("docstring")  MathFunctions::erf "double MathFunctions::erf(double arg)
+%feature("docstring")  Math::Bessel::erf "double Math::erf(double arg)
 
 Error function of real-valued argument. 
 ";
 
-%feature("docstring")  MathFunctions::Bessel_J0 "double MathFunctions::Bessel_J0(double x)
+%feature("docstring")  Math::Bessel::GeneratePoissonRandom "double Math::GeneratePoissonRandom(double average)
+";
+
+
+// File: namespaceMath_1_1Bessel.xml
+%feature("docstring")  Math::Bessel::J0 "double Math::Bessel::J0(double x)
 
 Bessel function of the first kind and order 0. 
 ";
 
-%feature("docstring")  MathFunctions::Bessel_J1 "double MathFunctions::Bessel_J1(double x)
+%feature("docstring")  Math::Bessel::J1 "double Math::Bessel::J1(double x)
 
 Bessel function of the first kind and order 1. 
 ";
 
-%feature("docstring")  MathFunctions::Bessel_J1c "double MathFunctions::Bessel_J1c(double x)
+%feature("docstring")  Math::Bessel::J1c "double Math::Bessel::J1c(double x)
 
-Bessel function Bessel_J1(x)/x. 
+Bessel function J1(x)/x. 
 ";
 
-%feature("docstring")  MathFunctions::Bessel_I0 "double MathFunctions::Bessel_I0(double x)
+%feature("docstring")  Math::Bessel::I0 "double Math::Bessel::I0(double x)
 
 Modified Bessel function of the first kind and order 0. 
 ";
 
-%feature("docstring")  MathFunctions::Bessel_J0 "complex_t MathFunctions::Bessel_J0(const complex_t z)
+%feature("docstring")  Math::Bessel::J0 "complex_t Math::Bessel::J0(const complex_t z)
 
 Complex Bessel function of the first kind and order 0. 
 ";
 
-%feature("docstring")  MathFunctions::Bessel_J1 "complex_t MathFunctions::Bessel_J1(const complex_t z)
+%feature("docstring")  Math::Bessel::J1 "complex_t Math::Bessel::J1(const complex_t z)
 
 Complex Bessel function of the first kind and order 1. 
 ";
 
-%feature("docstring")  MathFunctions::Bessel_J1c "complex_t MathFunctions::Bessel_J1c(const complex_t z)
-
-Complex Bessel function Bessel_J1(x)/x. 
-";
-
-%feature("docstring")  MathFunctions::FastFourierTransform "std::vector< complex_t > MathFunctions::FastFourierTransform(const std::vector< complex_t > &data, EFFTDirection tcase)
+%feature("docstring")  Math::Bessel::J1c "complex_t Math::Bessel::J1c(const complex_t z)
 
-simple (and unoptimized) wrapper function for the discrete fast Fourier transformation library (fftw3) 
-";
-
-%feature("docstring")  MathFunctions::FastFourierTransform "std::vector< complex_t > MathFunctions::FastFourierTransform(const std::vector< double > &data, EFFTDirection tcase)
-
-simple (and unoptimized) wrapper function for the discrete fast Fourier transformation library (fftw3); transforms real to complex 
-";
-
-%feature("docstring")  MathFunctions::ConvolveFFT "std::vector< complex_t > MathFunctions::ConvolveFFT(const std::vector< double > &signal, const std::vector< double > &resfunc)
-
-convolution of two real vectors of equal size 
-";
-
-%feature("docstring")  MathFunctions::GenerateUniformRandom "double MathFunctions::GenerateUniformRandom()
-";
-
-%feature("docstring")  MathFunctions::GenerateStandardNormalRandom "double MathFunctions::GenerateStandardNormalRandom()
-";
-
-%feature("docstring")  MathFunctions::GenerateNormalRandom "double MathFunctions::GenerateNormalRandom(double average, double std_dev)
-";
-
-%feature("docstring")  MathFunctions::GeneratePoissonRandom "double MathFunctions::GeneratePoissonRandom(double average)
+Complex Bessel function J1(x)/x. 
 ";
 
 
@@ -1604,6 +1541,12 @@ Returns exp(I*z), where I is the imaginary unit.
 // File: Assert_8h.xml
 
 
+// File: Bessel_8cpp.xml
+
+
+// File: Bessel_8h.xml
+
+
 // File: FileSystemUtils_8cpp.xml
 
 
diff --git a/auto/Wrap/doxygenCore.i b/auto/Wrap/doxygenCore.i
index 5cbdf2f416319422e172f50ac1dff396929e020d..1a91485096268f543732d7a2b7a5fcc4c4eaf3cf 100644
--- a/auto/Wrap/doxygenCore.i
+++ b/auto/Wrap/doxygenCore.i
@@ -420,9 +420,6 @@ C++ includes: DWBASingleComputation.h
 %feature("docstring")  DWBASingleComputation::~DWBASingleComputation "DWBASingleComputation::~DWBASingleComputation()
 ";
 
-%feature("docstring")  DWBASingleComputation::DWBASingleComputation "DWBASingleComputation::DWBASingleComputation(DWBASingleComputation &&other)
-";
-
 %feature("docstring")  DWBASingleComputation::setProgressHandler "void DWBASingleComputation::setProgressHandler(ProgressHandler *p_progress)
 ";
 
@@ -438,11 +435,6 @@ C++ includes: DWBASingleComputation.h
 %feature("docstring")  DWBASingleComputation::compute "void DWBASingleComputation::compute(SimulationElement &elem) const
 ";
 
-%feature("docstring")  DWBASingleComputation::regionMap "const std::map< size_t, std::vector< HomogeneousRegion > > & DWBASingleComputation::regionMap() const
-
-Retrieves a map of regions for the calculation of averaged layers. 
-";
-
 
 // File: classFitObjective.xml
 %feature("docstring") FitObjective "
@@ -2021,9 +2013,6 @@ Sets the beam polarization according to the given Bloch vector.
 %feature("docstring")  Simulation::setDetectorResolutionFunction "void Simulation::setDetectorResolutionFunction(const IResolutionFunction2D &resolution_function)
 ";
 
-%feature("docstring")  Simulation::removeDetectorResolutionFunction "void Simulation::removeDetectorResolutionFunction()
-";
-
 %feature("docstring")  Simulation::setAnalyzerProperties "void Simulation::setAnalyzerProperties(const kvector_t direction, double efficiency, double total_transmission)
 
 Sets the polarization analyzer characteristics of the detector. 
diff --git a/auto/Wrap/doxygenDevice.i b/auto/Wrap/doxygenDevice.i
index 0896fd73285617fc92988420c74bf8191db50fbc..b72cfdd6285377a8a403c9746d6dff371db1ec29 100644
--- a/auto/Wrap/doxygenDevice.i
+++ b/auto/Wrap/doxygenDevice.i
@@ -1184,7 +1184,7 @@ C++ includes: IHistogram.h
 %feature("docstring")  IHistogram::IHistogram "IHistogram::IHistogram(const IHistogram &other)
 ";
 
-%feature("docstring")  IHistogram::~IHistogram "virtual IHistogram::~IHistogram()
+%feature("docstring")  IHistogram::~IHistogram "virtual IHistogram::~IHistogram()=default
 ";
 
 %feature("docstring")  IHistogram::IHistogram "IHistogram::IHistogram(const IAxis &axis_x)
@@ -1276,15 +1276,15 @@ Returns the center of bin i of the x axis.
 Returns the center of bin i of the y axis. 
 ";
 
-%feature("docstring")  IHistogram::getBinContent "double IHistogram::getBinContent(size_t i) const
-
-Returns content (accumulated value) of bin i. 
+%feature("docstring")  IHistogram::getData "const OutputData< CumulativeValue > & IHistogram::getData() const
 ";
 
-%feature("docstring")  IHistogram::getData "const OutputData<CumulativeValue>& IHistogram::getData() const
+%feature("docstring")  IHistogram::getData "OutputData< CumulativeValue > & IHistogram::getData()
 ";
 
-%feature("docstring")  IHistogram::getData "OutputData<CumulativeValue>& IHistogram::getData()
+%feature("docstring")  IHistogram::getBinContent "double IHistogram::getBinContent(size_t i) const
+
+Returns content (accumulated value) of bin i. 
 ";
 
 %feature("docstring")  IHistogram::getBinContent "double IHistogram::getBinContent(size_t binx, size_t biny) const
@@ -1606,7 +1606,7 @@ Interface for reading strategy of  OutputData from file.
 C++ includes: OutputDataReadStrategy.h
 ";
 
-%feature("docstring")  IOutputDataReadStrategy::~IOutputDataReadStrategy "virtual IOutputDataReadStrategy::~IOutputDataReadStrategy()
+%feature("docstring")  IOutputDataReadStrategy::~IOutputDataReadStrategy "virtual IOutputDataReadStrategy::~IOutputDataReadStrategy()=default
 ";
 
 %feature("docstring")  IOutputDataReadStrategy::readOutputData "virtual OutputData<double>* IOutputDataReadStrategy::readOutputData(std::istream &input_stream)=0
@@ -1621,10 +1621,7 @@ Strategy interface to write OututData in file
 C++ includes: OutputDataWriteStrategy.h
 ";
 
-%feature("docstring")  IOutputDataWriteStrategy::IOutputDataWriteStrategy "IOutputDataWriteStrategy::IOutputDataWriteStrategy()
-";
-
-%feature("docstring")  IOutputDataWriteStrategy::~IOutputDataWriteStrategy "virtual IOutputDataWriteStrategy::~IOutputDataWriteStrategy()
+%feature("docstring")  IOutputDataWriteStrategy::~IOutputDataWriteStrategy "virtual IOutputDataWriteStrategy::~IOutputDataWriteStrategy()=default
 ";
 
 %feature("docstring")  IOutputDataWriteStrategy::writeOutputData "virtual void IOutputDataWriteStrategy::writeOutputData(const OutputData< double > &data, std::ostream &output_stream)=0
@@ -2870,7 +2867,7 @@ Returns underlying unit converter.
 %feature("docstring")  SimulationResult::size "size_t SimulationResult::size() const
 ";
 
-%feature("docstring")  SimulationResult::max "double & SimulationResult::max() const
+%feature("docstring")  SimulationResult::max "double SimulationResult::max() const
 ";
 
 %feature("docstring")  SimulationResult::empty "bool SimulationResult::empty() const
diff --git a/auto/Wrap/doxygenParam.i b/auto/Wrap/doxygenParam.i
index fdc889a7f8714366d92ae291da9b1c8d8ce0612e..ba1ccc8a8d389c1982f204d04e1bffaa30f281f1 100644
--- a/auto/Wrap/doxygenParam.i
+++ b/auto/Wrap/doxygenParam.i
@@ -1308,30 +1308,6 @@ C++ includes: ParameterSample.h
 ";
 
 
-// File: classPostorderStrategy.xml
-%feature("docstring") PostorderStrategy "
-
-Traverse tree; visit children before their parents.
-
-C++ includes: IterationStrategy.h
-";
-
-%feature("docstring")  PostorderStrategy::PostorderStrategy "PostorderStrategy::PostorderStrategy()
-";
-
-%feature("docstring")  PostorderStrategy::clone "PostorderStrategy * PostorderStrategy::clone() const
-";
-
-%feature("docstring")  PostorderStrategy::first "IteratorMemento PostorderStrategy::first(const INode *p_root)
-";
-
-%feature("docstring")  PostorderStrategy::next "void PostorderStrategy::next(IteratorMemento &iterator_stack) const
-";
-
-%feature("docstring")  PostorderStrategy::isDone "bool PostorderStrategy::isDone(IteratorMemento &iterator_stack) const
-";
-
-
 // File: classPreorderStrategy.xml
 %feature("docstring") PreorderStrategy "
 
@@ -1723,17 +1699,11 @@ Prints RealLimits in the form of argument (in the context of  ParameterDistribut
 %feature("docstring")  VisitNodesPreorder "void VisitNodesPreorder(const INode &node, INodeVisitor &visitor)
 ";
 
-%feature("docstring")  VisitNodesPostorder "void VisitNodesPostorder(const INode &node, INodeVisitor &visitor)
-";
-
 
 // File: INodeVisitor_8h.xml
 %feature("docstring")  VisitNodesPreorder "void VisitNodesPreorder(const INode &node, INodeVisitor &visitor)
 ";
 
-%feature("docstring")  VisitNodesPostorder "void VisitNodesPostorder(const INode &node, INodeVisitor &visitor)
-";
-
 
 // File: IterationStrategy_8cpp.xml
 
diff --git a/auto/Wrap/libBornAgainBase.py b/auto/Wrap/libBornAgainBase.py
index c56b4dfa91bb2dd6d74547cdb1a04cb0856d83de..917855e1cde810a9f16e9f16374947e7a26b6463 100644
--- a/auto/Wrap/libBornAgainBase.py
+++ b/auto/Wrap/libBornAgainBase.py
@@ -1756,188 +1756,42 @@ def deg2rad(angle):
     return _libBornAgainBase.deg2rad(angle)
 
 def StandardNormal(x):
-    r"""
-    StandardNormal(double x) -> double
-    double MathFunctions::StandardNormal(double x)
-
-    """
+    r"""StandardNormal(double x) -> double"""
     return _libBornAgainBase.StandardNormal(x)
 
 def Gaussian(x, average, std_dev):
-    r"""
-    Gaussian(double x, double average, double std_dev) -> double
-    double MathFunctions::Gaussian(double x, double average, double std_dev)
-
-    """
+    r"""Gaussian(double x, double average, double std_dev) -> double"""
     return _libBornAgainBase.Gaussian(x, average, std_dev)
 
 def IntegratedGaussian(x, average, std_dev):
-    r"""
-    IntegratedGaussian(double x, double average, double std_dev) -> double
-    double MathFunctions::IntegratedGaussian(double x, double average, double std_dev)
-
-    """
+    r"""IntegratedGaussian(double x, double average, double std_dev) -> double"""
     return _libBornAgainBase.IntegratedGaussian(x, average, std_dev)
 
 def cot(x):
-    r"""
-    cot(double x) -> double
-    double MathFunctions::cot(double x)
-
-    cotangent function:  $cot(x)\\equiv1/tan(x)$
-
-    """
+    r"""cot(double x) -> double"""
     return _libBornAgainBase.cot(x)
 
-def Si(x):
-    r"""
-    Si(double x) -> double
-    double MathFunctions::Si(double x)
-
-    Sine integral function:  $Si(x)\\equiv\\int_0^x du \\sin(u)/u$. 
-
-    """
-    return _libBornAgainBase.Si(x)
-
 def sinc(*args):
     r"""
     sinc(double x) -> double
     sinc(complex_t const z) -> complex_t
-    complex_t MathFunctions::sinc(const complex_t z)
-
-    Complex sinc function:  $sinc(x)\\equiv\\sin(x)/x$. 
-
     """
     return _libBornAgainBase.sinc(*args)
 
 def tanhc(z):
-    r"""
-    tanhc(complex_t const z) -> complex_t
-    complex_t MathFunctions::tanhc(const complex_t z)
-
-    Complex tanhc function:  $tanhc(x)\\equiv\\tanh(x)/x$. 
-
-    """
+    r"""tanhc(complex_t const z) -> complex_t"""
     return _libBornAgainBase.tanhc(z)
 
 def Laue(x, N):
-    r"""
-    Laue(double const x, size_t N) -> double
-    double MathFunctions::Laue(const double x, size_t N)
-
-    Real Laue function:  $Laue(x,N)\\equiv\\sin(Nx)/sin(x)$. 
-
-    """
+    r"""Laue(double const x, size_t N) -> double"""
     return _libBornAgainBase.Laue(x, N)
 
 def erf(arg):
-    r"""
-    erf(double arg) -> double
-    double MathFunctions::erf(double arg)
-
-    Error function of real-valued argument. 
-
-    """
+    r"""erf(double arg) -> double"""
     return _libBornAgainBase.erf(arg)
 
-def Bessel_I0(x):
-    r"""
-    Bessel_I0(double x) -> double
-    double MathFunctions::Bessel_I0(double x)
-
-    Modified Bessel function of the first kind and order 0. 
-
-    """
-    return _libBornAgainBase.Bessel_I0(x)
-
-def Bessel_J0(*args):
-    r"""
-    Bessel_J0(double x) -> double
-    Bessel_J0(complex_t const z) -> complex_t
-    complex_t MathFunctions::Bessel_J0(const complex_t z)
-
-    Complex Bessel function of the first kind and order 0. 
-
-    """
-    return _libBornAgainBase.Bessel_J0(*args)
-
-def Bessel_J1(*args):
-    r"""
-    Bessel_J1(double x) -> double
-    Bessel_J1(complex_t const z) -> complex_t
-    complex_t MathFunctions::Bessel_J1(const complex_t z)
-
-    Complex Bessel function of the first kind and order 1. 
-
-    """
-    return _libBornAgainBase.Bessel_J1(*args)
-
-def Bessel_J1c(*args):
-    r"""
-    Bessel_J1c(double x) -> double
-    Bessel_J1c(complex_t const z) -> complex_t
-    complex_t MathFunctions::Bessel_J1c(const complex_t z)
-
-    Complex Bessel function Bessel_J1(x)/x. 
-
-    """
-    return _libBornAgainBase.Bessel_J1c(*args)
-FORWARD_FFT = _libBornAgainBase.FORWARD_FFT
-
-BACKWARD_FFT = _libBornAgainBase.BACKWARD_FFT
-
-
-def FastFourierTransform(*args):
-    r"""
-    FastFourierTransform(vector_complex_t data, MathFunctions::EFFTDirection tcase) -> vector_complex_t
-    FastFourierTransform(vdouble1d_t data, MathFunctions::EFFTDirection tcase) -> vector_complex_t
-    std::vector< complex_t > MathFunctions::FastFourierTransform(const std::vector< double > &data, EFFTDirection tcase)
-
-    simple (and unoptimized) wrapper function for the discrete fast Fourier transformation library (fftw3); transforms real to complex 
-
-    """
-    return _libBornAgainBase.FastFourierTransform(*args)
-
-def ConvolveFFT(signal, resfunc):
-    r"""
-    ConvolveFFT(vdouble1d_t signal, vdouble1d_t resfunc) -> vector_complex_t
-    std::vector< complex_t > MathFunctions::ConvolveFFT(const std::vector< double > &signal, const std::vector< double > &resfunc)
-
-    convolution of two real vectors of equal size 
-
-    """
-    return _libBornAgainBase.ConvolveFFT(signal, resfunc)
-
-def GenerateUniformRandom():
-    r"""
-    GenerateUniformRandom() -> double
-    double MathFunctions::GenerateUniformRandom()
-
-    """
-    return _libBornAgainBase.GenerateUniformRandom()
-
-def GenerateStandardNormalRandom():
-    r"""
-    GenerateStandardNormalRandom() -> double
-    double MathFunctions::GenerateStandardNormalRandom()
-
-    """
-    return _libBornAgainBase.GenerateStandardNormalRandom()
-
-def GenerateNormalRandom(average, std_dev):
-    r"""
-    GenerateNormalRandom(double average, double std_dev) -> double
-    double MathFunctions::GenerateNormalRandom(double average, double std_dev)
-
-    """
-    return _libBornAgainBase.GenerateNormalRandom(average, std_dev)
-
 def GeneratePoissonRandom(average):
-    r"""
-    GeneratePoissonRandom(double average) -> double
-    double MathFunctions::GeneratePoissonRandom(double average)
-
-    """
+    r"""GeneratePoissonRandom(double average) -> double"""
     return _libBornAgainBase.GeneratePoissonRandom(average)
 class ThreadInfo(object):
     r"""
diff --git a/auto/Wrap/libBornAgainBase_wrap.cpp b/auto/Wrap/libBornAgainBase_wrap.cpp
index 559625f8d083a1aa30879a61d43e3ecf05410043..4942a7b52da378f19316f6e10fb006a5ffd17d69 100644
--- a/auto/Wrap/libBornAgainBase_wrap.cpp
+++ b/auto/Wrap/libBornAgainBase_wrap.cpp
@@ -24350,7 +24350,7 @@ SWIGINTERN PyObject *_wrap_StandardNormal(PyObject *SWIGUNUSEDPARM(self), PyObje
     SWIG_exception_fail(SWIG_ArgError(ecode1), "in method '" "StandardNormal" "', argument " "1"" of type '" "double""'");
   } 
   arg1 = static_cast< double >(val1);
-  result = (double)MathFunctions::StandardNormal(arg1);
+  result = (double)Math::StandardNormal(arg1);
   resultobj = SWIG_From_double(static_cast< double >(result));
   return resultobj;
 fail:
@@ -24388,7 +24388,7 @@ SWIGINTERN PyObject *_wrap_Gaussian(PyObject *SWIGUNUSEDPARM(self), PyObject *ar
     SWIG_exception_fail(SWIG_ArgError(ecode3), "in method '" "Gaussian" "', argument " "3"" of type '" "double""'");
   } 
   arg3 = static_cast< double >(val3);
-  result = (double)MathFunctions::Gaussian(arg1,arg2,arg3);
+  result = (double)Math::Gaussian(arg1,arg2,arg3);
   resultobj = SWIG_From_double(static_cast< double >(result));
   return resultobj;
 fail:
@@ -24426,7 +24426,7 @@ SWIGINTERN PyObject *_wrap_IntegratedGaussian(PyObject *SWIGUNUSEDPARM(self), Py
     SWIG_exception_fail(SWIG_ArgError(ecode3), "in method '" "IntegratedGaussian" "', argument " "3"" of type '" "double""'");
   } 
   arg3 = static_cast< double >(val3);
-  result = (double)MathFunctions::IntegratedGaussian(arg1,arg2,arg3);
+  result = (double)Math::IntegratedGaussian(arg1,arg2,arg3);
   resultobj = SWIG_From_double(static_cast< double >(result));
   return resultobj;
 fail:
@@ -24449,30 +24449,7 @@ SWIGINTERN PyObject *_wrap_cot(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
     SWIG_exception_fail(SWIG_ArgError(ecode1), "in method '" "cot" "', argument " "1"" of type '" "double""'");
   } 
   arg1 = static_cast< double >(val1);
-  result = (double)MathFunctions::cot(arg1);
-  resultobj = SWIG_From_double(static_cast< double >(result));
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_Si(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
-  PyObject *resultobj = 0;
-  double arg1 ;
-  double val1 ;
-  int ecode1 = 0 ;
-  PyObject *swig_obj[1] ;
-  double result;
-  
-  if (!args) SWIG_fail;
-  swig_obj[0] = args;
-  ecode1 = SWIG_AsVal_double(swig_obj[0], &val1);
-  if (!SWIG_IsOK(ecode1)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode1), "in method '" "Si" "', argument " "1"" of type '" "double""'");
-  } 
-  arg1 = static_cast< double >(val1);
-  result = (double)MathFunctions::Si(arg1);
+  result = (double)Math::cot(arg1);
   resultobj = SWIG_From_double(static_cast< double >(result));
   return resultobj;
 fail:
@@ -24493,7 +24470,7 @@ SWIGINTERN PyObject *_wrap_sinc__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize
     SWIG_exception_fail(SWIG_ArgError(ecode1), "in method '" "sinc" "', argument " "1"" of type '" "double""'");
   } 
   arg1 = static_cast< double >(val1);
-  result = (double)MathFunctions::sinc(arg1);
+  result = (double)Math::sinc(arg1);
   resultobj = SWIG_From_double(static_cast< double >(result));
   return resultobj;
 fail:
@@ -24514,7 +24491,7 @@ SWIGINTERN PyObject *_wrap_sinc__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize
     SWIG_exception_fail(SWIG_ArgError(ecode1), "in method '" "sinc" "', argument " "1"" of type '" "complex_t""'");
   } 
   arg1 = static_cast< complex_t >(val1);
-  result = MathFunctions::sinc(arg1);
+  result = Math::sinc(arg1);
   resultobj = SWIG_From_std_complex_Sl_double_Sg_(static_cast< std::complex<double> >(result));
   return resultobj;
 fail:
@@ -24554,8 +24531,8 @@ SWIGINTERN PyObject *_wrap_sinc(PyObject *self, PyObject *args) {
 fail:
   SWIG_Python_RaiseOrModifyTypeError("Wrong number or type of arguments for overloaded function 'sinc'.\n"
     "  Possible C/C++ prototypes are:\n"
-    "    MathFunctions::sinc(double)\n"
-    "    MathFunctions::sinc(complex_t const)\n");
+    "    Math::sinc(double)\n"
+    "    Math::sinc(complex_t const)\n");
   return 0;
 }
 
@@ -24575,7 +24552,7 @@ SWIGINTERN PyObject *_wrap_tanhc(PyObject *SWIGUNUSEDPARM(self), PyObject *args)
     SWIG_exception_fail(SWIG_ArgError(ecode1), "in method '" "tanhc" "', argument " "1"" of type '" "complex_t""'");
   } 
   arg1 = static_cast< complex_t >(val1);
-  result = MathFunctions::tanhc(arg1);
+  result = Math::tanhc(arg1);
   resultobj = SWIG_From_std_complex_Sl_double_Sg_(static_cast< std::complex<double> >(result));
   return resultobj;
 fail:
@@ -24605,7 +24582,7 @@ SWIGINTERN PyObject *_wrap_Laue(PyObject *SWIGUNUSEDPARM(self), PyObject *args)
     SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "Laue" "', argument " "2"" of type '" "size_t""'");
   } 
   arg2 = static_cast< size_t >(val2);
-  result = (double)MathFunctions::Laue(arg1,arg2);
+  result = (double)Math::Laue(arg1,arg2);
   resultobj = SWIG_From_double(static_cast< double >(result));
   return resultobj;
 fail:
@@ -24628,488 +24605,7 @@ SWIGINTERN PyObject *_wrap_erf(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
     SWIG_exception_fail(SWIG_ArgError(ecode1), "in method '" "erf" "', argument " "1"" of type '" "double""'");
   } 
   arg1 = static_cast< double >(val1);
-  result = (double)MathFunctions::erf(arg1);
-  resultobj = SWIG_From_double(static_cast< double >(result));
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_Bessel_J0__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
-  PyObject *resultobj = 0;
-  double arg1 ;
-  double val1 ;
-  int ecode1 = 0 ;
-  double result;
-  
-  if ((nobjs < 1) || (nobjs > 1)) SWIG_fail;
-  ecode1 = SWIG_AsVal_double(swig_obj[0], &val1);
-  if (!SWIG_IsOK(ecode1)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode1), "in method '" "Bessel_J0" "', argument " "1"" of type '" "double""'");
-  } 
-  arg1 = static_cast< double >(val1);
-  result = (double)MathFunctions::Bessel_J0(arg1);
-  resultobj = SWIG_From_double(static_cast< double >(result));
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_Bessel_J1__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
-  PyObject *resultobj = 0;
-  double arg1 ;
-  double val1 ;
-  int ecode1 = 0 ;
-  double result;
-  
-  if ((nobjs < 1) || (nobjs > 1)) SWIG_fail;
-  ecode1 = SWIG_AsVal_double(swig_obj[0], &val1);
-  if (!SWIG_IsOK(ecode1)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode1), "in method '" "Bessel_J1" "', argument " "1"" of type '" "double""'");
-  } 
-  arg1 = static_cast< double >(val1);
-  result = (double)MathFunctions::Bessel_J1(arg1);
-  resultobj = SWIG_From_double(static_cast< double >(result));
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_Bessel_J1c__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
-  PyObject *resultobj = 0;
-  double arg1 ;
-  double val1 ;
-  int ecode1 = 0 ;
-  double result;
-  
-  if ((nobjs < 1) || (nobjs > 1)) SWIG_fail;
-  ecode1 = SWIG_AsVal_double(swig_obj[0], &val1);
-  if (!SWIG_IsOK(ecode1)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode1), "in method '" "Bessel_J1c" "', argument " "1"" of type '" "double""'");
-  } 
-  arg1 = static_cast< double >(val1);
-  result = (double)MathFunctions::Bessel_J1c(arg1);
-  resultobj = SWIG_From_double(static_cast< double >(result));
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_Bessel_I0(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
-  PyObject *resultobj = 0;
-  double arg1 ;
-  double val1 ;
-  int ecode1 = 0 ;
-  PyObject *swig_obj[1] ;
-  double result;
-  
-  if (!args) SWIG_fail;
-  swig_obj[0] = args;
-  ecode1 = SWIG_AsVal_double(swig_obj[0], &val1);
-  if (!SWIG_IsOK(ecode1)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode1), "in method '" "Bessel_I0" "', argument " "1"" of type '" "double""'");
-  } 
-  arg1 = static_cast< double >(val1);
-  result = (double)MathFunctions::Bessel_I0(arg1);
-  resultobj = SWIG_From_double(static_cast< double >(result));
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_Bessel_J0__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
-  PyObject *resultobj = 0;
-  complex_t arg1 ;
-  std::complex< double > val1 ;
-  int ecode1 = 0 ;
-  complex_t result;
-  
-  if ((nobjs < 1) || (nobjs > 1)) SWIG_fail;
-  ecode1 = SWIG_AsVal_std_complex_Sl_double_Sg_(swig_obj[0], &val1);
-  if (!SWIG_IsOK(ecode1)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode1), "in method '" "Bessel_J0" "', argument " "1"" of type '" "complex_t""'");
-  } 
-  arg1 = static_cast< complex_t >(val1);
-  result = MathFunctions::Bessel_J0(arg1);
-  resultobj = SWIG_From_std_complex_Sl_double_Sg_(static_cast< std::complex<double> >(result));
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_Bessel_J0(PyObject *self, PyObject *args) {
-  Py_ssize_t argc;
-  PyObject *argv[2] = {
-    0
-  };
-  
-  if (!(argc = SWIG_Python_UnpackTuple(args, "Bessel_J0", 0, 1, argv))) SWIG_fail;
-  --argc;
-  if (argc == 1) {
-    int _v;
-    {
-      int res = SWIG_AsVal_double(argv[0], NULL);
-      _v = SWIG_CheckState(res);
-    }
-    if (_v) {
-      return _wrap_Bessel_J0__SWIG_0(self, argc, argv);
-    }
-  }
-  if (argc == 1) {
-    int _v;
-    {
-      int res = SWIG_AsVal_std_complex_Sl_double_Sg_(argv[0], NULL);
-      _v = SWIG_CheckState(res);
-    }
-    if (_v) {
-      return _wrap_Bessel_J0__SWIG_1(self, argc, argv);
-    }
-  }
-  
-fail:
-  SWIG_Python_RaiseOrModifyTypeError("Wrong number or type of arguments for overloaded function 'Bessel_J0'.\n"
-    "  Possible C/C++ prototypes are:\n"
-    "    MathFunctions::Bessel_J0(double)\n"
-    "    MathFunctions::Bessel_J0(complex_t const)\n");
-  return 0;
-}
-
-
-SWIGINTERN PyObject *_wrap_Bessel_J1__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
-  PyObject *resultobj = 0;
-  complex_t arg1 ;
-  std::complex< double > val1 ;
-  int ecode1 = 0 ;
-  complex_t result;
-  
-  if ((nobjs < 1) || (nobjs > 1)) SWIG_fail;
-  ecode1 = SWIG_AsVal_std_complex_Sl_double_Sg_(swig_obj[0], &val1);
-  if (!SWIG_IsOK(ecode1)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode1), "in method '" "Bessel_J1" "', argument " "1"" of type '" "complex_t""'");
-  } 
-  arg1 = static_cast< complex_t >(val1);
-  result = MathFunctions::Bessel_J1(arg1);
-  resultobj = SWIG_From_std_complex_Sl_double_Sg_(static_cast< std::complex<double> >(result));
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_Bessel_J1(PyObject *self, PyObject *args) {
-  Py_ssize_t argc;
-  PyObject *argv[2] = {
-    0
-  };
-  
-  if (!(argc = SWIG_Python_UnpackTuple(args, "Bessel_J1", 0, 1, argv))) SWIG_fail;
-  --argc;
-  if (argc == 1) {
-    int _v;
-    {
-      int res = SWIG_AsVal_double(argv[0], NULL);
-      _v = SWIG_CheckState(res);
-    }
-    if (_v) {
-      return _wrap_Bessel_J1__SWIG_0(self, argc, argv);
-    }
-  }
-  if (argc == 1) {
-    int _v;
-    {
-      int res = SWIG_AsVal_std_complex_Sl_double_Sg_(argv[0], NULL);
-      _v = SWIG_CheckState(res);
-    }
-    if (_v) {
-      return _wrap_Bessel_J1__SWIG_1(self, argc, argv);
-    }
-  }
-  
-fail:
-  SWIG_Python_RaiseOrModifyTypeError("Wrong number or type of arguments for overloaded function 'Bessel_J1'.\n"
-    "  Possible C/C++ prototypes are:\n"
-    "    MathFunctions::Bessel_J1(double)\n"
-    "    MathFunctions::Bessel_J1(complex_t const)\n");
-  return 0;
-}
-
-
-SWIGINTERN PyObject *_wrap_Bessel_J1c__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
-  PyObject *resultobj = 0;
-  complex_t arg1 ;
-  std::complex< double > val1 ;
-  int ecode1 = 0 ;
-  complex_t result;
-  
-  if ((nobjs < 1) || (nobjs > 1)) SWIG_fail;
-  ecode1 = SWIG_AsVal_std_complex_Sl_double_Sg_(swig_obj[0], &val1);
-  if (!SWIG_IsOK(ecode1)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode1), "in method '" "Bessel_J1c" "', argument " "1"" of type '" "complex_t""'");
-  } 
-  arg1 = static_cast< complex_t >(val1);
-  result = MathFunctions::Bessel_J1c(arg1);
-  resultobj = SWIG_From_std_complex_Sl_double_Sg_(static_cast< std::complex<double> >(result));
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_Bessel_J1c(PyObject *self, PyObject *args) {
-  Py_ssize_t argc;
-  PyObject *argv[2] = {
-    0
-  };
-  
-  if (!(argc = SWIG_Python_UnpackTuple(args, "Bessel_J1c", 0, 1, argv))) SWIG_fail;
-  --argc;
-  if (argc == 1) {
-    int _v;
-    {
-      int res = SWIG_AsVal_double(argv[0], NULL);
-      _v = SWIG_CheckState(res);
-    }
-    if (_v) {
-      return _wrap_Bessel_J1c__SWIG_0(self, argc, argv);
-    }
-  }
-  if (argc == 1) {
-    int _v;
-    {
-      int res = SWIG_AsVal_std_complex_Sl_double_Sg_(argv[0], NULL);
-      _v = SWIG_CheckState(res);
-    }
-    if (_v) {
-      return _wrap_Bessel_J1c__SWIG_1(self, argc, argv);
-    }
-  }
-  
-fail:
-  SWIG_Python_RaiseOrModifyTypeError("Wrong number or type of arguments for overloaded function 'Bessel_J1c'.\n"
-    "  Possible C/C++ prototypes are:\n"
-    "    MathFunctions::Bessel_J1c(double)\n"
-    "    MathFunctions::Bessel_J1c(complex_t const)\n");
-  return 0;
-}
-
-
-SWIGINTERN PyObject *_wrap_FastFourierTransform__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
-  PyObject *resultobj = 0;
-  std::vector< complex_t,std::allocator< complex_t > > *arg1 = 0 ;
-  MathFunctions::EFFTDirection arg2 ;
-  int res1 = SWIG_OLDOBJ ;
-  int val2 ;
-  int ecode2 = 0 ;
-  std::vector< complex_t,std::allocator< complex_t > > result;
-  
-  if ((nobjs < 2) || (nobjs > 2)) SWIG_fail;
-  {
-    std::vector< std::complex< double >,std::allocator< std::complex< double > > > *ptr = (std::vector< std::complex< double >,std::allocator< std::complex< double > > > *)0;
-    res1 = swig::asptr(swig_obj[0], &ptr);
-    if (!SWIG_IsOK(res1)) {
-      SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "FastFourierTransform" "', argument " "1"" of type '" "std::vector< complex_t,std::allocator< complex_t > > const &""'"); 
-    }
-    if (!ptr) {
-      SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "FastFourierTransform" "', argument " "1"" of type '" "std::vector< complex_t,std::allocator< complex_t > > const &""'"); 
-    }
-    arg1 = ptr;
-  }
-  ecode2 = SWIG_AsVal_int(swig_obj[1], &val2);
-  if (!SWIG_IsOK(ecode2)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "FastFourierTransform" "', argument " "2"" of type '" "MathFunctions::EFFTDirection""'");
-  } 
-  arg2 = static_cast< MathFunctions::EFFTDirection >(val2);
-  result = MathFunctions::FastFourierTransform((std::vector< std::complex< double >,std::allocator< std::complex< double > > > const &)*arg1,arg2);
-  resultobj = swig::from(static_cast< std::vector< std::complex< double >,std::allocator< std::complex< double > > > >(result));
-  if (SWIG_IsNewObj(res1)) delete arg1;
-  return resultobj;
-fail:
-  if (SWIG_IsNewObj(res1)) delete arg1;
-  return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_FastFourierTransform__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
-  PyObject *resultobj = 0;
-  std::vector< double,std::allocator< double > > *arg1 = 0 ;
-  MathFunctions::EFFTDirection arg2 ;
-  int res1 = SWIG_OLDOBJ ;
-  int val2 ;
-  int ecode2 = 0 ;
-  std::vector< complex_t,std::allocator< complex_t > > result;
-  
-  if ((nobjs < 2) || (nobjs > 2)) SWIG_fail;
-  {
-    std::vector< double,std::allocator< double > > *ptr = (std::vector< double,std::allocator< double > > *)0;
-    res1 = swig::asptr(swig_obj[0], &ptr);
-    if (!SWIG_IsOK(res1)) {
-      SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "FastFourierTransform" "', argument " "1"" of type '" "std::vector< double,std::allocator< double > > const &""'"); 
-    }
-    if (!ptr) {
-      SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "FastFourierTransform" "', argument " "1"" of type '" "std::vector< double,std::allocator< double > > const &""'"); 
-    }
-    arg1 = ptr;
-  }
-  ecode2 = SWIG_AsVal_int(swig_obj[1], &val2);
-  if (!SWIG_IsOK(ecode2)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "FastFourierTransform" "', argument " "2"" of type '" "MathFunctions::EFFTDirection""'");
-  } 
-  arg2 = static_cast< MathFunctions::EFFTDirection >(val2);
-  result = MathFunctions::FastFourierTransform((std::vector< double,std::allocator< double > > const &)*arg1,arg2);
-  resultobj = swig::from(static_cast< std::vector< std::complex< double >,std::allocator< std::complex< double > > > >(result));
-  if (SWIG_IsNewObj(res1)) delete arg1;
-  return resultobj;
-fail:
-  if (SWIG_IsNewObj(res1)) delete arg1;
-  return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_FastFourierTransform(PyObject *self, PyObject *args) {
-  Py_ssize_t argc;
-  PyObject *argv[3] = {
-    0
-  };
-  
-  if (!(argc = SWIG_Python_UnpackTuple(args, "FastFourierTransform", 0, 2, argv))) SWIG_fail;
-  --argc;
-  if (argc == 2) {
-    int _v;
-    int res = swig::asptr(argv[0], (std::vector< std::complex< double >,std::allocator< std::complex< double > > >**)(0));
-    _v = SWIG_CheckState(res);
-    if (_v) {
-      {
-        int res = SWIG_AsVal_int(argv[1], NULL);
-        _v = SWIG_CheckState(res);
-      }
-      if (_v) {
-        return _wrap_FastFourierTransform__SWIG_0(self, argc, argv);
-      }
-    }
-  }
-  if (argc == 2) {
-    int _v;
-    int res = swig::asptr(argv[0], (std::vector< double,std::allocator< double > >**)(0));
-    _v = SWIG_CheckState(res);
-    if (_v) {
-      {
-        int res = SWIG_AsVal_int(argv[1], NULL);
-        _v = SWIG_CheckState(res);
-      }
-      if (_v) {
-        return _wrap_FastFourierTransform__SWIG_1(self, argc, argv);
-      }
-    }
-  }
-  
-fail:
-  SWIG_Python_RaiseOrModifyTypeError("Wrong number or type of arguments for overloaded function 'FastFourierTransform'.\n"
-    "  Possible C/C++ prototypes are:\n"
-    "    MathFunctions::FastFourierTransform(std::vector< complex_t,std::allocator< complex_t > > const &,MathFunctions::EFFTDirection)\n"
-    "    MathFunctions::FastFourierTransform(std::vector< double,std::allocator< double > > const &,MathFunctions::EFFTDirection)\n");
-  return 0;
-}
-
-
-SWIGINTERN PyObject *_wrap_ConvolveFFT(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
-  PyObject *resultobj = 0;
-  std::vector< double,std::allocator< double > > *arg1 = 0 ;
-  std::vector< double,std::allocator< double > > *arg2 = 0 ;
-  int res1 = SWIG_OLDOBJ ;
-  int res2 = SWIG_OLDOBJ ;
-  PyObject *swig_obj[2] ;
-  std::vector< complex_t,std::allocator< complex_t > > result;
-  
-  if (!SWIG_Python_UnpackTuple(args, "ConvolveFFT", 2, 2, swig_obj)) SWIG_fail;
-  {
-    std::vector< double,std::allocator< double > > *ptr = (std::vector< double,std::allocator< double > > *)0;
-    res1 = swig::asptr(swig_obj[0], &ptr);
-    if (!SWIG_IsOK(res1)) {
-      SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "ConvolveFFT" "', argument " "1"" of type '" "std::vector< double,std::allocator< double > > const &""'"); 
-    }
-    if (!ptr) {
-      SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "ConvolveFFT" "', argument " "1"" of type '" "std::vector< double,std::allocator< double > > const &""'"); 
-    }
-    arg1 = ptr;
-  }
-  {
-    std::vector< double,std::allocator< double > > *ptr = (std::vector< double,std::allocator< double > > *)0;
-    res2 = swig::asptr(swig_obj[1], &ptr);
-    if (!SWIG_IsOK(res2)) {
-      SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "ConvolveFFT" "', argument " "2"" of type '" "std::vector< double,std::allocator< double > > const &""'"); 
-    }
-    if (!ptr) {
-      SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "ConvolveFFT" "', argument " "2"" of type '" "std::vector< double,std::allocator< double > > const &""'"); 
-    }
-    arg2 = ptr;
-  }
-  result = MathFunctions::ConvolveFFT((std::vector< double,std::allocator< double > > const &)*arg1,(std::vector< double,std::allocator< double > > const &)*arg2);
-  resultobj = swig::from(static_cast< std::vector< std::complex< double >,std::allocator< std::complex< double > > > >(result));
-  if (SWIG_IsNewObj(res1)) delete arg1;
-  if (SWIG_IsNewObj(res2)) delete arg2;
-  return resultobj;
-fail:
-  if (SWIG_IsNewObj(res1)) delete arg1;
-  if (SWIG_IsNewObj(res2)) delete arg2;
-  return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_GenerateUniformRandom(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
-  PyObject *resultobj = 0;
-  double result;
-  
-  if (!SWIG_Python_UnpackTuple(args, "GenerateUniformRandom", 0, 0, 0)) SWIG_fail;
-  result = (double)MathFunctions::GenerateUniformRandom();
-  resultobj = SWIG_From_double(static_cast< double >(result));
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_GenerateStandardNormalRandom(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
-  PyObject *resultobj = 0;
-  double result;
-  
-  if (!SWIG_Python_UnpackTuple(args, "GenerateStandardNormalRandom", 0, 0, 0)) SWIG_fail;
-  result = (double)MathFunctions::GenerateStandardNormalRandom();
-  resultobj = SWIG_From_double(static_cast< double >(result));
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_GenerateNormalRandom(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
-  PyObject *resultobj = 0;
-  double arg1 ;
-  double arg2 ;
-  double val1 ;
-  int ecode1 = 0 ;
-  double val2 ;
-  int ecode2 = 0 ;
-  PyObject *swig_obj[2] ;
-  double result;
-  
-  if (!SWIG_Python_UnpackTuple(args, "GenerateNormalRandom", 2, 2, swig_obj)) SWIG_fail;
-  ecode1 = SWIG_AsVal_double(swig_obj[0], &val1);
-  if (!SWIG_IsOK(ecode1)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode1), "in method '" "GenerateNormalRandom" "', argument " "1"" of type '" "double""'");
-  } 
-  arg1 = static_cast< double >(val1);
-  ecode2 = SWIG_AsVal_double(swig_obj[1], &val2);
-  if (!SWIG_IsOK(ecode2)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "GenerateNormalRandom" "', argument " "2"" of type '" "double""'");
-  } 
-  arg2 = static_cast< double >(val2);
-  result = (double)MathFunctions::GenerateNormalRandom(arg1,arg2);
+  result = (double)Math::erf(arg1);
   resultobj = SWIG_From_double(static_cast< double >(result));
   return resultobj;
 fail:
@@ -25132,7 +24628,7 @@ SWIGINTERN PyObject *_wrap_GeneratePoissonRandom(PyObject *SWIGUNUSEDPARM(self),
     SWIG_exception_fail(SWIG_ArgError(ecode1), "in method '" "GeneratePoissonRandom" "', argument " "1"" of type '" "double""'");
   } 
   arg1 = static_cast< double >(val1);
-  result = (double)MathFunctions::GeneratePoissonRandom(arg1);
+  result = (double)Math::GeneratePoissonRandom(arg1);
   resultobj = SWIG_From_double(static_cast< double >(result));
   return resultobj;
 fail:
@@ -33779,130 +33275,18 @@ static PyMethodDef SwigMethods[] = {
 		"double Units::deg2rad(double angle)\n"
 		"\n"
 		""},
-	 { "StandardNormal", _wrap_StandardNormal, METH_O, "\n"
-		"StandardNormal(double x) -> double\n"
-		"double MathFunctions::StandardNormal(double x)\n"
-		"\n"
-		""},
-	 { "Gaussian", _wrap_Gaussian, METH_VARARGS, "\n"
-		"Gaussian(double x, double average, double std_dev) -> double\n"
-		"double MathFunctions::Gaussian(double x, double average, double std_dev)\n"
-		"\n"
-		""},
-	 { "IntegratedGaussian", _wrap_IntegratedGaussian, METH_VARARGS, "\n"
-		"IntegratedGaussian(double x, double average, double std_dev) -> double\n"
-		"double MathFunctions::IntegratedGaussian(double x, double average, double std_dev)\n"
-		"\n"
-		""},
-	 { "cot", _wrap_cot, METH_O, "\n"
-		"cot(double x) -> double\n"
-		"double MathFunctions::cot(double x)\n"
-		"\n"
-		"cotangent function:  $cot(x)\\\\equiv1/tan(x)$\n"
-		"\n"
-		""},
-	 { "Si", _wrap_Si, METH_O, "\n"
-		"Si(double x) -> double\n"
-		"double MathFunctions::Si(double x)\n"
-		"\n"
-		"Sine integral function:  $Si(x)\\\\equiv\\\\int_0^x du \\\\sin(u)/u$. \n"
-		"\n"
-		""},
+	 { "StandardNormal", _wrap_StandardNormal, METH_O, "StandardNormal(double x) -> double"},
+	 { "Gaussian", _wrap_Gaussian, METH_VARARGS, "Gaussian(double x, double average, double std_dev) -> double"},
+	 { "IntegratedGaussian", _wrap_IntegratedGaussian, METH_VARARGS, "IntegratedGaussian(double x, double average, double std_dev) -> double"},
+	 { "cot", _wrap_cot, METH_O, "cot(double x) -> double"},
 	 { "sinc", _wrap_sinc, METH_VARARGS, "\n"
 		"sinc(double x) -> double\n"
 		"sinc(complex_t const z) -> complex_t\n"
-		"complex_t MathFunctions::sinc(const complex_t z)\n"
-		"\n"
-		"Complex sinc function:  $sinc(x)\\\\equiv\\\\sin(x)/x$. \n"
-		"\n"
-		""},
-	 { "tanhc", _wrap_tanhc, METH_O, "\n"
-		"tanhc(complex_t const z) -> complex_t\n"
-		"complex_t MathFunctions::tanhc(const complex_t z)\n"
-		"\n"
-		"Complex tanhc function:  $tanhc(x)\\\\equiv\\\\tanh(x)/x$. \n"
-		"\n"
-		""},
-	 { "Laue", _wrap_Laue, METH_VARARGS, "\n"
-		"Laue(double const x, size_t N) -> double\n"
-		"double MathFunctions::Laue(const double x, size_t N)\n"
-		"\n"
-		"Real Laue function:  $Laue(x,N)\\\\equiv\\\\sin(Nx)/sin(x)$. \n"
-		"\n"
-		""},
-	 { "erf", _wrap_erf, METH_O, "\n"
-		"erf(double arg) -> double\n"
-		"double MathFunctions::erf(double arg)\n"
-		"\n"
-		"Error function of real-valued argument. \n"
-		"\n"
-		""},
-	 { "Bessel_I0", _wrap_Bessel_I0, METH_O, "\n"
-		"Bessel_I0(double x) -> double\n"
-		"double MathFunctions::Bessel_I0(double x)\n"
-		"\n"
-		"Modified Bessel function of the first kind and order 0. \n"
-		"\n"
-		""},
-	 { "Bessel_J0", _wrap_Bessel_J0, METH_VARARGS, "\n"
-		"Bessel_J0(double x) -> double\n"
-		"Bessel_J0(complex_t const z) -> complex_t\n"
-		"complex_t MathFunctions::Bessel_J0(const complex_t z)\n"
-		"\n"
-		"Complex Bessel function of the first kind and order 0. \n"
-		"\n"
-		""},
-	 { "Bessel_J1", _wrap_Bessel_J1, METH_VARARGS, "\n"
-		"Bessel_J1(double x) -> double\n"
-		"Bessel_J1(complex_t const z) -> complex_t\n"
-		"complex_t MathFunctions::Bessel_J1(const complex_t z)\n"
-		"\n"
-		"Complex Bessel function of the first kind and order 1. \n"
-		"\n"
-		""},
-	 { "Bessel_J1c", _wrap_Bessel_J1c, METH_VARARGS, "\n"
-		"Bessel_J1c(double x) -> double\n"
-		"Bessel_J1c(complex_t const z) -> complex_t\n"
-		"complex_t MathFunctions::Bessel_J1c(const complex_t z)\n"
-		"\n"
-		"Complex Bessel function Bessel_J1(x)/x. \n"
-		"\n"
-		""},
-	 { "FastFourierTransform", _wrap_FastFourierTransform, METH_VARARGS, "\n"
-		"FastFourierTransform(vector_complex_t data, MathFunctions::EFFTDirection tcase) -> vector_complex_t\n"
-		"FastFourierTransform(vdouble1d_t data, MathFunctions::EFFTDirection tcase) -> vector_complex_t\n"
-		"std::vector< complex_t > MathFunctions::FastFourierTransform(const std::vector< double > &data, EFFTDirection tcase)\n"
-		"\n"
-		"simple (and unoptimized) wrapper function for the discrete fast Fourier transformation library (fftw3); transforms real to complex \n"
-		"\n"
-		""},
-	 { "ConvolveFFT", _wrap_ConvolveFFT, METH_VARARGS, "\n"
-		"ConvolveFFT(vdouble1d_t signal, vdouble1d_t resfunc) -> vector_complex_t\n"
-		"std::vector< complex_t > MathFunctions::ConvolveFFT(const std::vector< double > &signal, const std::vector< double > &resfunc)\n"
-		"\n"
-		"convolution of two real vectors of equal size \n"
-		"\n"
-		""},
-	 { "GenerateUniformRandom", _wrap_GenerateUniformRandom, METH_NOARGS, "\n"
-		"GenerateUniformRandom() -> double\n"
-		"double MathFunctions::GenerateUniformRandom()\n"
-		"\n"
-		""},
-	 { "GenerateStandardNormalRandom", _wrap_GenerateStandardNormalRandom, METH_NOARGS, "\n"
-		"GenerateStandardNormalRandom() -> double\n"
-		"double MathFunctions::GenerateStandardNormalRandom()\n"
-		"\n"
-		""},
-	 { "GenerateNormalRandom", _wrap_GenerateNormalRandom, METH_VARARGS, "\n"
-		"GenerateNormalRandom(double average, double std_dev) -> double\n"
-		"double MathFunctions::GenerateNormalRandom(double average, double std_dev)\n"
-		"\n"
-		""},
-	 { "GeneratePoissonRandom", _wrap_GeneratePoissonRandom, METH_O, "\n"
-		"GeneratePoissonRandom(double average) -> double\n"
-		"double MathFunctions::GeneratePoissonRandom(double average)\n"
-		"\n"
 		""},
+	 { "tanhc", _wrap_tanhc, METH_O, "tanhc(complex_t const z) -> complex_t"},
+	 { "Laue", _wrap_Laue, METH_VARARGS, "Laue(double const x, size_t N) -> double"},
+	 { "erf", _wrap_erf, METH_O, "erf(double arg) -> double"},
+	 { "GeneratePoissonRandom", _wrap_GeneratePoissonRandom, METH_O, "GeneratePoissonRandom(double average) -> double"},
 	 { "new_ThreadInfo", _wrap_new_ThreadInfo, METH_NOARGS, "\n"
 		"new_ThreadInfo() -> ThreadInfo\n"
 		"ThreadInfo::ThreadInfo()\n"
@@ -35856,8 +35240,6 @@ SWIG_init(void) {
   SWIG_addvarlink(globals, "deg", Swig_var_deg_get, Swig_var_deg_set);
   SWIG_addvarlink(globals, "tesla", Swig_var_tesla_get, Swig_var_tesla_set);
   SWIG_addvarlink(globals, "gauss", Swig_var_gauss_get, Swig_var_gauss_set);
-  SWIG_Python_SetConstant(d, "FORWARD_FFT",SWIG_From_int(static_cast< int >(MathFunctions::FORWARD_FFT)));
-  SWIG_Python_SetConstant(d, "BACKWARD_FFT",SWIG_From_int(static_cast< int >(MathFunctions::BACKWARD_FFT)));
 #if PY_VERSION_HEX >= 0x03000000
   return m;
 #else
diff --git a/auto/Wrap/libBornAgainCore.py b/auto/Wrap/libBornAgainCore.py
index a51214000b7d909b02bf839cbde5f6b96ec50e44..4f87c61aca4bc9766b065cdb2393c4ab0134db5e 100644
--- a/auto/Wrap/libBornAgainCore.py
+++ b/auto/Wrap/libBornAgainCore.py
@@ -3689,14 +3689,6 @@ class Simulation(libBornAgainBase.ICloneable, libBornAgainParam.INode):
         """
         return _libBornAgainCore.Simulation_setDetectorResolutionFunction(self, resolution_function)
 
-    def removeDetectorResolutionFunction(self):
-        r"""
-        removeDetectorResolutionFunction(Simulation self)
-        void Simulation::removeDetectorResolutionFunction()
-
-        """
-        return _libBornAgainCore.Simulation_removeDetectorResolutionFunction(self)
-
     def setAnalyzerProperties(self, direction, efficiency, total_transmission):
         r"""
         setAnalyzerProperties(Simulation self, kvector_t direction, double efficiency, double total_transmission)
diff --git a/auto/Wrap/libBornAgainCore_wrap.cpp b/auto/Wrap/libBornAgainCore_wrap.cpp
index c474f1ac296e522c58784b6efabdceb99fa6f5aa..43e8de257f395921f5a83cc9fb2ba9ff4a3eea27 100644
--- a/auto/Wrap/libBornAgainCore_wrap.cpp
+++ b/auto/Wrap/libBornAgainCore_wrap.cpp
@@ -39532,28 +39532,6 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_Simulation_removeDetectorResolutionFunction(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
-  PyObject *resultobj = 0;
-  Simulation *arg1 = (Simulation *) 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_Simulation, 0 |  0 );
-  if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "Simulation_removeDetectorResolutionFunction" "', argument " "1"" of type '" "Simulation *""'"); 
-  }
-  arg1 = reinterpret_cast< Simulation * >(argp1);
-  (arg1)->removeDetectorResolutionFunction();
-  resultobj = SWIG_Py_Void();
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
 SWIGINTERN PyObject *_wrap_Simulation_setAnalyzerProperties(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
   Simulation *arg1 = (Simulation *) 0 ;
@@ -44027,11 +44005,6 @@ static PyMethodDef SwigMethods[] = {
 		"void Simulation::setDetectorResolutionFunction(const IResolutionFunction2D &resolution_function)\n"
 		"\n"
 		""},
-	 { "Simulation_removeDetectorResolutionFunction", _wrap_Simulation_removeDetectorResolutionFunction, METH_O, "\n"
-		"Simulation_removeDetectorResolutionFunction(Simulation self)\n"
-		"void Simulation::removeDetectorResolutionFunction()\n"
-		"\n"
-		""},
 	 { "Simulation_setAnalyzerProperties", _wrap_Simulation_setAnalyzerProperties, METH_VARARGS, "\n"
 		"Simulation_setAnalyzerProperties(Simulation self, kvector_t direction, double efficiency, double total_transmission)\n"
 		"void Simulation::setAnalyzerProperties(const kvector_t direction, double efficiency, double total_transmission)\n"
diff --git a/auto/Wrap/libBornAgainDevice.py b/auto/Wrap/libBornAgainDevice.py
index 6a4529adb44be29552a1703d0a1bb676eaaf2a74..d099567a6641d090c9a8e13e3077cf5fe50a0144 100644
--- a/auto/Wrap/libBornAgainDevice.py
+++ b/auto/Wrap/libBornAgainDevice.py
@@ -5430,7 +5430,7 @@ class IHistogram(object):
         r"""
         getData(IHistogram self) -> OutputData< CumulativeValue > const
         getData(IHistogram self) -> OutputData< CumulativeValue > &
-        OutputData<CumulativeValue>& IHistogram::getData()
+        OutputData< CumulativeValue > & IHistogram::getData()
 
         """
         return _libBornAgainDevice.IHistogram_getData(self, *args)
@@ -6127,7 +6127,7 @@ class SimulationResult(object):
     def max(self):
         r"""
         max(SimulationResult self) -> double
-        double & SimulationResult::max() const
+        double SimulationResult::max() const
 
         """
         return _libBornAgainDevice.SimulationResult_max(self)
diff --git a/auto/Wrap/libBornAgainDevice_wrap.cpp b/auto/Wrap/libBornAgainDevice_wrap.cpp
index cf6146414700ecf6c3a811ff41b7129b4732df12..82a5cfbaa2cc9090bc78bb55da4991d95c490f36 100644
--- a/auto/Wrap/libBornAgainDevice_wrap.cpp
+++ b/auto/Wrap/libBornAgainDevice_wrap.cpp
@@ -41434,35 +41434,6 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_IHistogram_getBinContent__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
-  PyObject *resultobj = 0;
-  IHistogram *arg1 = (IHistogram *) 0 ;
-  size_t arg2 ;
-  void *argp1 = 0 ;
-  int res1 = 0 ;
-  size_t val2 ;
-  int ecode2 = 0 ;
-  double result;
-  
-  if ((nobjs < 2) || (nobjs > 2)) SWIG_fail;
-  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_IHistogram, 0 |  0 );
-  if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IHistogram_getBinContent" "', argument " "1"" of type '" "IHistogram const *""'"); 
-  }
-  arg1 = reinterpret_cast< IHistogram * >(argp1);
-  ecode2 = SWIG_AsVal_size_t(swig_obj[1], &val2);
-  if (!SWIG_IsOK(ecode2)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "IHistogram_getBinContent" "', argument " "2"" of type '" "size_t""'");
-  } 
-  arg2 = static_cast< size_t >(val2);
-  result = (double)((IHistogram const *)arg1)->getBinContent(arg2);
-  resultobj = SWIG_From_double(static_cast< double >(result));
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
 SWIGINTERN PyObject *_wrap_IHistogram_getData__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
   IHistogram *arg1 = (IHistogram *) 0 ;
@@ -41541,6 +41512,35 @@ fail:
 }
 
 
+SWIGINTERN PyObject *_wrap_IHistogram_getBinContent__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
+  PyObject *resultobj = 0;
+  IHistogram *arg1 = (IHistogram *) 0 ;
+  size_t arg2 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  size_t val2 ;
+  int ecode2 = 0 ;
+  double result;
+  
+  if ((nobjs < 2) || (nobjs > 2)) SWIG_fail;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_IHistogram, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IHistogram_getBinContent" "', argument " "1"" of type '" "IHistogram const *""'"); 
+  }
+  arg1 = reinterpret_cast< IHistogram * >(argp1);
+  ecode2 = SWIG_AsVal_size_t(swig_obj[1], &val2);
+  if (!SWIG_IsOK(ecode2)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "IHistogram_getBinContent" "', argument " "2"" of type '" "size_t""'");
+  } 
+  arg2 = static_cast< size_t >(val2);
+  result = (double)((IHistogram const *)arg1)->getBinContent(arg2);
+  resultobj = SWIG_From_double(static_cast< double >(result));
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
 SWIGINTERN PyObject *_wrap_IHistogram_getBinContent__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
   IHistogram *arg1 = (IHistogram *) 0 ;
@@ -48185,7 +48185,7 @@ static PyMethodDef SwigMethods[] = {
 	 { "IsGISAXSDetector_swiginit", IsGISAXSDetector_swiginit, METH_VARARGS, NULL},
 	 { "delete_IHistogram", _wrap_delete_IHistogram, METH_O, "\n"
 		"delete_IHistogram(IHistogram self)\n"
-		"virtual IHistogram::~IHistogram()\n"
+		"virtual IHistogram::~IHistogram()=default\n"
 		"\n"
 		""},
 	 { "IHistogram_clone", _wrap_IHistogram_clone, METH_O, "\n"
@@ -48308,7 +48308,7 @@ static PyMethodDef SwigMethods[] = {
 	 { "IHistogram_getData", _wrap_IHistogram_getData, METH_VARARGS, "\n"
 		"IHistogram_getData(IHistogram self) -> OutputData< CumulativeValue > const\n"
 		"IHistogram_getData(IHistogram self) -> OutputData< CumulativeValue > &\n"
-		"OutputData<CumulativeValue>& IHistogram::getData()\n"
+		"OutputData< CumulativeValue > & IHistogram::getData()\n"
 		"\n"
 		""},
 	 { "IHistogram_getBinContent", _wrap_IHistogram_getBinContent, METH_VARARGS, "\n"
@@ -48713,7 +48713,7 @@ static PyMethodDef SwigMethods[] = {
 		""},
 	 { "SimulationResult_max", _wrap_SimulationResult_max, METH_O, "\n"
 		"SimulationResult_max(SimulationResult self) -> double\n"
-		"double & SimulationResult::max() const\n"
+		"double SimulationResult::max() const\n"
 		"\n"
 		""},
 	 { "SimulationResult_empty", _wrap_SimulationResult_empty, METH_O, "\n"
diff --git a/auto/Wrap/libBornAgainParam.py b/auto/Wrap/libBornAgainParam.py
index 8db6b9400fbc59f5fad98eba85ba6a7698d75e69..b7a0664e24659993b570e12491e5a7d78e957473 100644
--- a/auto/Wrap/libBornAgainParam.py
+++ b/auto/Wrap/libBornAgainParam.py
@@ -3455,14 +3455,6 @@ def VisitNodesPreorder(node, visitor):
 
     """
     return _libBornAgainParam.VisitNodesPreorder(node, visitor)
-
-def VisitNodesPostorder(node, visitor):
-    r"""
-    VisitNodesPostorder(INode node, INodeVisitor visitor)
-    void VisitNodesPostorder(const INode &node, INodeVisitor &visitor)
-
-    """
-    return _libBornAgainParam.VisitNodesPostorder(node, visitor)
 class IDistribution1D(libBornAgainBase.ICloneable, INode):
     r"""
 
diff --git a/auto/Wrap/libBornAgainParam_wrap.cpp b/auto/Wrap/libBornAgainParam_wrap.cpp
index d6478b5d6b6a28c8c5192047ef26d11aa717bc88..d701261cbe9aab7319ca57e3e54851d73f19597a 100644
--- a/auto/Wrap/libBornAgainParam_wrap.cpp
+++ b/auto/Wrap/libBornAgainParam_wrap.cpp
@@ -42336,41 +42336,6 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_VisitNodesPostorder(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
-  PyObject *resultobj = 0;
-  INode *arg1 = 0 ;
-  INodeVisitor *arg2 = 0 ;
-  void *argp1 = 0 ;
-  int res1 = 0 ;
-  void *argp2 = 0 ;
-  int res2 = 0 ;
-  PyObject *swig_obj[2] ;
-  
-  if (!SWIG_Python_UnpackTuple(args, "VisitNodesPostorder", 2, 2, swig_obj)) SWIG_fail;
-  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1, SWIGTYPE_p_INode,  0  | 0);
-  if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "VisitNodesPostorder" "', argument " "1"" of type '" "INode const &""'"); 
-  }
-  if (!argp1) {
-    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "VisitNodesPostorder" "', argument " "1"" of type '" "INode const &""'"); 
-  }
-  arg1 = reinterpret_cast< INode * >(argp1);
-  res2 = SWIG_ConvertPtr(swig_obj[1], &argp2, SWIGTYPE_p_INodeVisitor,  0 );
-  if (!SWIG_IsOK(res2)) {
-    SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "VisitNodesPostorder" "', argument " "2"" of type '" "INodeVisitor &""'"); 
-  }
-  if (!argp2) {
-    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "VisitNodesPostorder" "', argument " "2"" of type '" "INodeVisitor &""'"); 
-  }
-  arg2 = reinterpret_cast< INodeVisitor * >(argp2);
-  VisitNodesPostorder((INode const &)*arg1,*arg2);
-  resultobj = SWIG_Py_Void();
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
 SWIGINTERN PyObject *_wrap_IDistribution1D_clone(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
   IDistribution1D *arg1 = (IDistribution1D *) 0 ;
@@ -51374,11 +51339,6 @@ static PyMethodDef SwigMethods[] = {
 		"void VisitNodesPreorder(const INode &node, INodeVisitor &visitor)\n"
 		"\n"
 		""},
-	 { "VisitNodesPostorder", _wrap_VisitNodesPostorder, METH_VARARGS, "\n"
-		"VisitNodesPostorder(INode node, INodeVisitor visitor)\n"
-		"void VisitNodesPostorder(const INode &node, INodeVisitor &visitor)\n"
-		"\n"
-		""},
 	 { "IDistribution1D_clone", _wrap_IDistribution1D_clone, METH_O, "\n"
 		"IDistribution1D_clone(IDistribution1D self) -> IDistribution1D\n"
 		"virtual IDistribution1D* IDistribution1D::clone() const =0\n"
diff --git a/auto/Wrap/libBornAgainSample.py b/auto/Wrap/libBornAgainSample.py
index 5634d795bb7def527b8bf158ac1adde9fcd8b1b8..3c6988c5fda786e0333741bb43737b8e301be141 100644
--- a/auto/Wrap/libBornAgainSample.py
+++ b/auto/Wrap/libBornAgainSample.py
@@ -3602,11 +3602,6 @@ class IRotation(libBornAgainBase.ICloneable, libBornAgainParam.INode):
         r"""createRotation(Transform3D const & transform) -> IRotation"""
         return _libBornAgainSample.IRotation_createRotation(transform)
 
-    @staticmethod
-    def createIdentity():
-        r"""createIdentity() -> IRotation"""
-        return _libBornAgainSample.IRotation_createIdentity()
-
     def clone(self):
         r"""
         clone(IRotation self) -> IRotation
@@ -3669,10 +3664,6 @@ def IRotation_createRotation(transform):
     r"""IRotation_createRotation(Transform3D const & transform) -> IRotation"""
     return _libBornAgainSample.IRotation_createRotation(transform)
 
-def IRotation_createIdentity():
-    r"""IRotation_createIdentity() -> IRotation"""
-    return _libBornAgainSample.IRotation_createIdentity()
-
 
 def createProduct(left, right):
     r"""
diff --git a/auto/Wrap/libBornAgainSample_wrap.cpp b/auto/Wrap/libBornAgainSample_wrap.cpp
index 79ad88f16c36bc7469273ae4da6bf47f95721dd0..b4b7786a4ee479f798afaf55032573558cb1154d 100644
--- a/auto/Wrap/libBornAgainSample_wrap.cpp
+++ b/auto/Wrap/libBornAgainSample_wrap.cpp
@@ -40412,19 +40412,6 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_IRotation_createIdentity(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
-  PyObject *resultobj = 0;
-  IRotation *result = 0 ;
-  
-  if (!SWIG_Python_UnpackTuple(args, "IRotation_createIdentity", 0, 0, 0)) SWIG_fail;
-  result = (IRotation *)IRotation::createIdentity();
-  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_IRotation, 0 |  0 );
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
 SWIGINTERN PyObject *_wrap_IRotation_clone(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
   IRotation *arg1 = (IRotation *) 0 ;
@@ -72466,7 +72453,6 @@ static PyMethodDef SwigMethods[] = {
 		""},
 	 { "IFormFactorDecorator_swigregister", IFormFactorDecorator_swigregister, METH_O, NULL},
 	 { "IRotation_createRotation", _wrap_IRotation_createRotation, METH_O, "IRotation_createRotation(Transform3D const & transform) -> IRotation"},
-	 { "IRotation_createIdentity", _wrap_IRotation_createIdentity, METH_NOARGS, "IRotation_createIdentity() -> IRotation"},
 	 { "IRotation_clone", _wrap_IRotation_clone, METH_O, "\n"
 		"IRotation_clone(IRotation self) -> IRotation\n"
 		"virtual IRotation* IRotation::clone() const =0\n"
diff --git a/cmake/commons/CoverageFunction.cmake b/cmake/commons/CoverageFunction.cmake
index 17f3aaebdec671fe0bd9dfae5741a48e9954985d..c5579044cf6581b74a308efaec4a65be96bfc5dd 100644
--- a/cmake/commons/CoverageFunction.cmake
+++ b/cmake/commons/CoverageFunction.cmake
@@ -22,31 +22,31 @@ find_program(GENHTML_COMMAND genhtml)
 if(GCOV_COMMAND)
     set(CMAKE_GCOV_FOUND TRUE)
 else()
-    message(ERROR "could not locate gcov executable, cannot add code coverage custom target!")
+    message(FATAL_ERROR "could not locate gcov executable, cannot add code coverage custom target!")
 endif()
 
 if(LCOV_COMMAND)
     set(CMAKE_LCOV_FOUND TRUE)
 else()
-    message(ERROR "could not locate lcov executable, cannot add code coverage custom target!")
+    message(FATAL_ERROR "could not locate lcov executable, cannot add code coverage custom target!")
 endif()
 
 if(GENHTML_COMMAND)
     set(CMAKE_GENHTML_FOUND TRUE)
 else()
-    message(ERROR "could not locate genhtml executable, cannot add code coverage custom target!")
+    message(FATAL_ERROR "could not locate genhtml executable, cannot add code coverage custom target!")
 endif()
 
 if(CMAKE_C_COMPILER_ID STREQUAL "GNU")
     string(APPEND CMAKE_C_FLAGS " -coverage")
 else()
-    message(ERROR "C compiler is not GNU; don't know how to set code coverage flags!")
+    message(FATAL_ERROR "C compiler is not GNU; don't know how to set code coverage flags!")
 endif()
 
 if(CMAKE_CXX_COMPILER_ID STREQUAL "GNU")
     string(APPEND CMAKE_CXX_FLAGS " -coverage")
 else()
-    message(ERROR "CXX compiler is not GNU; don't know how to set code coverage flags!")
+    message(FATAL_ERROR "CXX compiler is not GNU; don't know how to set code coverage flags!")
 endif()
 
 # function to add a coverage target