diff --git a/Core/Aggregate/FTDecayFunctions.cpp b/Core/Aggregate/FTDecayFunctions.cpp
index b1b7b4b34aeaf93b63791639ae7cc06fddd96723..ef1562d199daf84dae6a25f34da60065123a3324 100644
--- a/Core/Aggregate/FTDecayFunctions.cpp
+++ b/Core/Aggregate/FTDecayFunctions.cpp
@@ -116,7 +116,7 @@ FTDecayFunction1DCosine::FTDecayFunction1DCosine(double omega)
 double FTDecayFunction1DCosine::evaluate(double q) const
 {
     double qw = std::abs(q*m_omega);
-    if (std::abs(qw/Pi::PI-1.0) < Numeric::double_epsilon) {
+    if (std::abs(qw/Pi::PI-1.0) < std::numeric_limits<double>::epsilon()) {
         return m_omega/2.0;
     }
     else {
diff --git a/Core/Aggregate/FTDistributions1D.cpp b/Core/Aggregate/FTDistributions1D.cpp
index d7e32203f7001ac0170e7f47e5efcf82fb818e22..186d69662df7e69a334f5b59a5c482264e6d7277 100644
--- a/Core/Aggregate/FTDistributions1D.cpp
+++ b/Core/Aggregate/FTDistributions1D.cpp
@@ -17,10 +17,10 @@
 #include "BornAgainNamespace.h"
 #include "IntegratorReal.h"
 #include "MathFunctions.h"
-#include "Numeric.h"
 #include "ParameterPool.h"
 #include "Pi.h"
 #include "RealParameter.h"
+#include <limits>
 
 void IFTDistribution1D::print(std::ostream& ostr) const
 {
@@ -94,7 +94,7 @@ FTDistribution1DCosine::FTDistribution1DCosine(double omega)
 double FTDistribution1DCosine::evaluate(double q) const
 {
     double qw = std::abs(q*m_omega);
-    if (std::abs(qw/Pi::PI-1.0) < Numeric::double_epsilon)
+    if (std::abs(qw/Pi::PI-1.0) < std::numeric_limits<double>::epsilon())
         return 0.5;
     return MathFunctions::sinc(qw)/(1.0-qw*qw/Pi::PI/Pi::PI);
 }
diff --git a/Core/Aggregate/FTDistributions2D.cpp b/Core/Aggregate/FTDistributions2D.cpp
index d9666deb903d6318cdf09893b8cd8a8d058a3721..ae72d771eb21afaef6abceb87ed0e7310aeae1a8 100644
--- a/Core/Aggregate/FTDistributions2D.cpp
+++ b/Core/Aggregate/FTDistributions2D.cpp
@@ -17,10 +17,10 @@
 #include "BornAgainNamespace.h"
 #include "IntegratorReal.h"
 #include "MathFunctions.h"
-#include "Numeric.h"
 #include "ParameterPool.h"
 #include "Pi.h"
 #include "RealParameter.h"
+#include <limits>
 
 IFTDistribution2D::IFTDistribution2D(
     double coherence_length_x, double coherence_length_y, double gamma, double delta)
@@ -100,7 +100,7 @@ FTDistribution2DCone::FTDistribution2DCone(
 double FTDistribution2DCone::evaluate(double qx, double qy) const
 {
     double scaled_q = std::sqrt(sumsq(qx,qy));
-    if (scaled_q<Numeric::double_epsilon)
+    if (scaled_q<std::numeric_limits<double>::epsilon())
         return 1.0 - 3.0*scaled_q*scaled_q/40.0;
     auto integrator = make_integrator_real(this, &FTDistribution2DCone::coneIntegrand2);
     double integral = integrator->integrate(0.0, scaled_q);
diff --git a/Core/Aggregate/InterferenceFunction2DParaCrystal.cpp b/Core/Aggregate/InterferenceFunction2DParaCrystal.cpp
index 3d64448d67e84538c77b876544cfff26b90f34d2..31c87aff2d51dd7e55289b1ff667475e27d27664 100644
--- a/Core/Aggregate/InterferenceFunction2DParaCrystal.cpp
+++ b/Core/Aggregate/InterferenceFunction2DParaCrystal.cpp
@@ -18,7 +18,6 @@
 #include "Exceptions.h"
 #include "ISampleVisitor.h"
 #include "IntegratorReal.h"
-#include "Numeric.h"
 #include "ParameterPool.h"
 #include "Pi.h"
 #include "RealParameter.h"
@@ -202,7 +201,7 @@ double InterferenceFunction2DParaCrystal::interference1D(
     if (n<1) {
         result = ((1.0 + fp)/(1.0 - fp)).real();
     } else {
-        if (std::norm(1.0-fp) < Numeric::double_epsilon )
+        if (std::norm(1.0-fp) < std::numeric_limits<double>::epsilon() )
             result = nd;
         // for (1-fp)*nd small, take the series expansion to second order in nd*(1-fp)
         else if (std::abs(1.0-fp)*nd < 2e-4) {
diff --git a/Core/Aggregate/InterferenceFunctionRadialParaCrystal.cpp b/Core/Aggregate/InterferenceFunctionRadialParaCrystal.cpp
index 20be390db46258ebcc6f898fc76d771f0e2bd163..aeaae8e8d91974386632275e429e290b2e77ac30 100644
--- a/Core/Aggregate/InterferenceFunctionRadialParaCrystal.cpp
+++ b/Core/Aggregate/InterferenceFunctionRadialParaCrystal.cpp
@@ -17,7 +17,6 @@
 #include "BornAgainNamespace.h"
 #include "Exceptions.h"
 #include "ISampleVisitor.h"
-#include "Numeric.h"
 #include "ParameterPool.h"
 #include "RealParameter.h"
 #include <limits>
@@ -83,7 +82,7 @@ double InterferenceFunctionRadialParaCrystal::evaluate(const kvector_t q) const
     if (n<1) {
         result = ((1.0 + fp)/(1.0 - fp)).real();
     } else {
-        if (std::norm(1.0-fp) < Numeric::double_epsilon ) {
+        if (std::norm(1.0-fp) < std::numeric_limits<double>::epsilon() ) {
             result = nd;
         }
         // for (1-fp)*nd small, take the series expansion to second order in nd*(1-fp)
diff --git a/Core/HardParticle/FormFactorCone.cpp b/Core/HardParticle/FormFactorCone.cpp
index 22c0fcdafa4df1fc2e2f22c7f60178479fe1033b..27c012667b20b4744b8c05de0bcdb1c3e774815b 100644
--- a/Core/HardParticle/FormFactorCone.cpp
+++ b/Core/HardParticle/FormFactorCone.cpp
@@ -17,9 +17,9 @@
 #include "BornAgainNamespace.h"
 #include "Exceptions.h"
 #include "MathFunctions.h"
-#include "Numeric.h"
 #include "Pi.h"
 #include "RealParameter.h"
+#include <limits>
 
 //! @param radius of circular base
 //! @param height of frustum
@@ -58,7 +58,7 @@ complex_t FormFactorCone::Integrand(double Z) const
 complex_t FormFactorCone::evaluate_for_q(const cvector_t q) const
 {
     m_q = q;
-    if ( std::abs(m_q.mag()) < Numeric::double_epsilon) {
+    if ( std::abs(m_q.mag()) < std::numeric_limits<double>::epsilon()) {
         double R = m_radius;
         double H = m_height;
         double tga = std::tan(m_alpha);
diff --git a/Core/HardParticle/FormFactorFullSpheroid.cpp b/Core/HardParticle/FormFactorFullSpheroid.cpp
index 8e253543fd276c48265e1e3994be9a94cbb7586a..9cdfd98d94cb23532dad69c4fa0c8688b9bfa530 100644
--- a/Core/HardParticle/FormFactorFullSpheroid.cpp
+++ b/Core/HardParticle/FormFactorFullSpheroid.cpp
@@ -16,9 +16,9 @@
 #include "FormFactorFullSpheroid.h"
 #include "BornAgainNamespace.h"
 #include "MathFunctions.h"
-#include "Numeric.h"
 #include "Pi.h"
 #include "RealParameter.h"
+#include <limits>
 
 //! @param radius of the two equal axes
 //! @param height total height of the spheroid, i.e. twice the radius of the third axis
@@ -51,10 +51,8 @@ complex_t FormFactorFullSpheroid::evaluate_for_q(const cvector_t q) const
     double R = m_radius;
     m_q = q;
 
-    if (std::abs(m_q.mag()) <= Numeric::double_epsilon) {
+    if (std::abs(m_q.mag()) <= std::numeric_limits<double>::epsilon())
         return Pi::PI2*R*R*H/3.;
-    } else {
-        complex_t qzH_half = H/2*q.z();
-        return 4 * Pi::PI * mP_integrator->integrate(0.0, H/2.0) * exp_I(qzH_half);
-    }
+    complex_t qzH_half = H/2*q.z();
+    return 4 * Pi::PI * mP_integrator->integrate(0.0, H/2.0) * exp_I(qzH_half);
 }
diff --git a/Core/HardParticle/FormFactorHemiEllipsoid.cpp b/Core/HardParticle/FormFactorHemiEllipsoid.cpp
index 7baa217a8034b1c8e246df8fb08cc13dbf960f36..51da0225f0caa0d7609eaba4481af18ed2500d1a 100644
--- a/Core/HardParticle/FormFactorHemiEllipsoid.cpp
+++ b/Core/HardParticle/FormFactorHemiEllipsoid.cpp
@@ -16,9 +16,9 @@
 #include "FormFactorHemiEllipsoid.h"
 #include "BornAgainNamespace.h"
 #include "MathFunctions.h"
-#include "Numeric.h"
 #include "Pi.h"
 #include "RealParameter.h"
+#include <limits>
 
 //! @param radius_x half length of one horizontal main axes
 //! @param radius_y half length of the other horizontal main axes
@@ -64,9 +64,7 @@ complex_t FormFactorHemiEllipsoid::evaluate_for_q(const cvector_t q) const
      double W = m_radius_y;
      double H = m_height;
 
-     if (std::abs(m_q.mag()) <= Numeric::double_epsilon) {
+     if (std::abs(m_q.mag()) <= std::numeric_limits<double>::epsilon())
          return Pi::PI2*R*W*H/3.;
-     } else {
-         return Pi::PI2*mP_integrator->integrate(0.,H );
-     }
+     return Pi::PI2*mP_integrator->integrate(0.,H );
 }
diff --git a/Core/HardParticle/FormFactorTruncatedSphere.cpp b/Core/HardParticle/FormFactorTruncatedSphere.cpp
index 41ee234b706e4a1eb2beaa60139d854b6d966c08..4e831a63d668c5c3f24c3fd5b47815c7c053eb78 100644
--- a/Core/HardParticle/FormFactorTruncatedSphere.cpp
+++ b/Core/HardParticle/FormFactorTruncatedSphere.cpp
@@ -18,11 +18,11 @@
 #include "Exceptions.h"
 #include "Limits.h"
 #include "MathFunctions.h"
-#include "Numeric.h"
 #include "Pi.h"
 #include "RealParameter.h"
+#include <limits>
 
-using namespace  BornAgain;
+using namespace BornAgain;
 
 FormFactorTruncatedSphere::FormFactorTruncatedSphere(double radius, double height)
     : m_radius(radius), m_height(height)
@@ -63,7 +63,7 @@ complex_t FormFactorTruncatedSphere::Integrand(double Z) const
 complex_t FormFactorTruncatedSphere::evaluate_for_q(const cvector_t q) const
 {
     m_q = q;
-    if ( std::abs(q.mag()) < Numeric::double_epsilon) {
+    if ( std::abs(q.mag()) < std::numeric_limits<double>::epsilon()) {
         double HdivR = m_height/m_radius;
         return Pi::PI/3.*m_radius*m_radius*m_radius
                 *(3.*HdivR -1. - (HdivR - 1.)*(HdivR - 1.)*(HdivR - 1.));
diff --git a/Core/HardParticle/FormFactorTruncatedSpheroid.cpp b/Core/HardParticle/FormFactorTruncatedSpheroid.cpp
index af08c277cd9bab3a791322dff8b61c986403cf32..dbda19c0d21c9e8adf0882c39ba5f60281a26894 100644
--- a/Core/HardParticle/FormFactorTruncatedSpheroid.cpp
+++ b/Core/HardParticle/FormFactorTruncatedSpheroid.cpp
@@ -17,9 +17,9 @@
 #include "BornAgainNamespace.h"
 #include "Exceptions.h"
 #include "MathFunctions.h"
-#include "Numeric.h"
 #include "Pi.h"
 #include "RealParameter.h"
+#include <limits>
 
 FormFactorTruncatedSpheroid::FormFactorTruncatedSpheroid(
     double radius, double height, double height_flattening)
@@ -68,10 +68,8 @@ complex_t FormFactorTruncatedSpheroid::evaluate_for_q(const cvector_t q) const
     double fp = m_height_flattening;
     m_q = q;
 
-    if (std::abs(m_q.mag()) <= Numeric::double_epsilon) {
+    if (std::abs(m_q.mag()) <= std::numeric_limits<double>::epsilon())
         return Pi::PI*R*H*H/fp*(1.-H/(3.*fp*R));
-    } else {
-        complex_t z_part    =  std::exp(complex_t(0.0, 1.0)*m_q.z()*(H-fp*R));
-        return Pi::PI2 * z_part *mP_integrator->integrate(fp*R-H,fp*R );
-    }
+    complex_t z_part    =  std::exp(complex_t(0.0, 1.0)*m_q.z()*(H-fp*R));
+    return Pi::PI2 * z_part *mP_integrator->integrate(fp*R-H,fp*R );
 }
diff --git a/Core/Mask/Line.cpp b/Core/Mask/Line.cpp
index 3b6191d83c44d4debf4cda4b518d123f7719f0dd..5ca4ae9a9f1f67975b5b5191c68b1462350a5411 100644
--- a/Core/Mask/Line.cpp
+++ b/Core/Mask/Line.cpp
@@ -17,6 +17,7 @@
 #include "Bin.h"
 #include "Macros.h"
 #include "Numeric.h"
+#include <limits>
 GCC_DIAG_OFF(unused-parameter)
 #include <boost/geometry.hpp>
 #include <boost/geometry/geometries/point_xy.hpp>
@@ -42,7 +43,7 @@ bool Line::contains(double x, double y) const
 
     double d = distance(p, line);
 
-    return d<Numeric::double_epsilon;
+    return d<std::numeric_limits<double>::epsilon();
 }
 
 // Calculates if line crosses the box made out of our bins.
diff --git a/Core/SoftParticle/FormFactorSphereUniformRadius.cpp b/Core/SoftParticle/FormFactorSphereUniformRadius.cpp
index 0479eb536d7826d5bda134bac745a55404957683..19d6fd888b1f0d9cf0ef8668bd8e75dffa0ef56f 100644
--- a/Core/SoftParticle/FormFactorSphereUniformRadius.cpp
+++ b/Core/SoftParticle/FormFactorSphereUniformRadius.cpp
@@ -16,9 +16,9 @@
 #include "FormFactorSphereUniformRadius.h"
 #include "BornAgainNamespace.h"
 #include "Exceptions.h"
-#include "Numeric.h"
 #include "Pi.h"
 #include "RealParameter.h"
+#include <limits>
 
 FormFactorSphereUniformRadius::FormFactorSphereUniformRadius(double mean,
         double full_width)
@@ -40,13 +40,13 @@ complex_t FormFactorSphereUniformRadius::evaluate_for_q(const cvector_t q) const
     double W = m_full_width;
     double q2 = std::norm(q.x()) + std::norm(q.y()) + std::norm(q.z());
     double q_r = std::sqrt(q2);
-    if (q_r*R < Numeric::double_epsilon)
+    if (q_r*R < std::numeric_limits<double>::epsilon())
         return (4.0*Pi::PI*R*R*R + Pi::PI*R*W*W)/3.0;
     double qR = q_r*R;
     double qW = q_r*W;
     double nominator = 4*Pi::PI*( 4*std::sin(qR)*std::sin(qW/2.0)
-                                     - qW*std::cos(qW/2.0)*std::sin(qR)
-                                     - 2.0*qR*std::cos(qR)*std::sin(qW/2.0) );
+                                  - qW*std::cos(qW/2.0)*std::sin(qR)
+                                  - 2.0*qR*std::cos(qR)*std::sin(qW/2.0) );
     return nominator/(q2*q2*W);
 }
 
diff --git a/Core/Tools/MathFunctions.cpp b/Core/Tools/MathFunctions.cpp
index a429139923ff882158dcd942e6bff31b2e7288a4..71a433b3bb5723b3455fce03387e1785dacc83ff 100644
--- a/Core/Tools/MathFunctions.cpp
+++ b/Core/Tools/MathFunctions.cpp
@@ -14,7 +14,6 @@
 // ************************************************************************** //
 
 #include "MathFunctions.h"
-#include "Numeric.h"
 #include "Pi.h"
 #include <gsl/gsl_sf_bessel.h>
 #include <gsl/gsl_sf_erf.h>
@@ -23,6 +22,7 @@
 #include <fftw3.h>
 #include <chrono>
 #include <cstring>
+#include <limits>
 #include <random>
 #include <stdexcept> // need overlooked by g++ 5.4
 
@@ -75,7 +75,7 @@ complex_t MathFunctions::sinc(const complex_t z)  // Sin(x)/x
 
 complex_t MathFunctions::tanhc(const complex_t z)  // tanh(x)/x
 {
-    if(std::abs(z)<Numeric::double_epsilon)
+    if(std::abs(z)<std::numeric_limits<double>::epsilon())
         return 1.0;
     return std::tanh(z)/z;
 }
@@ -84,7 +84,7 @@ complex_t MathFunctions::Laue(const complex_t z, size_t N) // Exp(iNx/2)*Sin((N+
 {
     if (N==0)
         return 1.0;
-    if(std::abs(z)<Numeric::double_epsilon)
+    if(std::abs(z)<std::numeric_limits<double>::epsilon())
         return N+1.0;
     return exp_I(N/2.0*z)*std::sin(z*(N+1.0)/2.0)/std::sin(z/2.0);
 }
diff --git a/Core/Tools/Numeric.cpp b/Core/Tools/Numeric.cpp
index ec87d422646b1afa752f192768f1021ac7fd33b0..ac1db42e1b15035680dd2bbf2ae027f677e66797 100644
--- a/Core/Tools/Numeric.cpp
+++ b/Core/Tools/Numeric.cpp
@@ -22,22 +22,23 @@
 
 namespace Numeric {
 
-//! Returns true if two doubles agree within epsilon*tolerance
+//! Returns true if two doubles agree within epsilon*tolerance.
 bool areAlmostEqual(double a, double b, double tolerance)
 {
     constexpr double eps = std::numeric_limits<double>::epsilon();
     return std::abs(a-b) <= eps * std::max( tolerance*eps, std::max(1., tolerance)*std::abs(b) );
 }
 
-//! calculates safe relative difference |(a-b)/b|
+//! Returns the safe relative difference, which is |(a-b)/b| except in special cases.
 double get_relative_difference(double a, double b)
 {
+    constexpr double eps = std::numeric_limits<double>::epsilon();
     // return 0.0 if relative error smaller than epsilon
-    if (std::abs(a-b) <= double_epsilon*std::abs(b))
+    if (std::abs(a-b) <= eps*std::abs(b))
         return 0.0;
     // for small numbers, divide by epsilon (to avoid catastrophic cancellation)
-    if (std::abs(b) <= double_epsilon)
-        return std::abs((a-b)/double_epsilon);
+    if (std::abs(b) <= eps)
+        return std::abs((a-b)/eps);
     return std::abs((a-b)/b);
 }
 
diff --git a/Core/Tools/Numeric.h b/Core/Tools/Numeric.h
index f1580ca7cbf664b2de469fc0d74f856b1bd51503..b1f073b03d6d2782e3e96e3e9d305efc7c71c3ac 100644
--- a/Core/Tools/Numeric.h
+++ b/Core/Tools/Numeric.h
@@ -23,12 +23,8 @@
 
 namespace Numeric {
 
-static const double double_epsilon = std::numeric_limits<double>::epsilon();
-
-//! compare two doubles
 bool BA_CORE_API_ areAlmostEqual(double a, double b, double tolerance_factor=1.0);
 
-//! calculates safe relative difference |(a-b)/b|
 double BA_CORE_API_ get_relative_difference(double a, double b);
 
 } // Numeric namespace