replace double_epsilon by direct use of <limits>

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 {

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

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

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

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)

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

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

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

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

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

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.

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

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

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

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