Reduce machine epsilon and add a distance epsilon

Reduce the value of the egs++ floating-point epsilon from 10e-10 to
a value that is 4 times (2 bits away) from the true machine epsilon, and
appropriately for single and double precision versions of the code.

Add a distanceEpsilon constant that is independent of the machine
epsilon, and much larger, but still entirely sufficient for physical
dimensions (well below atomic dimensions in double precision). By
default, the base geometry sets its boundaryTolerance to this
constant. Its value is close to the previous value of epsilon (1e-10),
which should curb any change in floating point error triggers in the
geometry library at this point.

Earlier in #177, the epsilon constant was introduced to remove
hard-coded floating point comparison guards. The 1e-10 value was chosen
because it matched most of the hard-coded values in the geometry
library. However, this epsilon is now starting to be used as a general
float guard (see #525 for example), where 1e-10 is not appropriate, as
it effectively reduces the floating-point precision of the computation.
Hence the need for independent epsilon values for computations and for
geometrical purposes.

HEN_HOUSE/egs++/egs_base_geometry.cpp

@@ -430,7 +430,7 @@ EGS_BaseGeometry::EGS_BaseGeometry(const string &Name) : nreg(0), name(Name),
  region_media(0), med(-1), has_rho_scaling(false), rhor(0),
  has_B_scaling(false), has_Ref_rho(false), bfactor(0), rhoRef(1.0),
  nref(0), debug(false), is_convex(true), bproperty(0), bp_array(0),
- boundaryTolerance(epsilon) {
+ boundaryTolerance(distanceEpsilon) {
 
  halfBoundaryTolerance = boundaryTolerance/2.;
  if (!egs_geometries.size()) {

HEN_HOUSE/egs++/egs_functions.h

@@ -49,9 +49,41 @@ using namespace std;
  *
  * The epsilon constant is a small number that can be used
  * when comparing two floating point numbers that may be close
- * in value.
+ * in value. It is set to 4 times the machine epsilon (2 bits away),
+ * to allow for some leeway in floating point error accumulation.
+ * It is meant to be used as a float guard for comparisons
+ * involving numbers that are of order 1.0; typically used to
+ * test if a float x vanishes, as in fabs(x) < epsilon.
  */
-const EGS_Float epsilon = 1e-10;
+#ifdef SINGLE
+ const EGS_Float epsilon = 1.0/(1<<21);
+#else
+ const EGS_Float epsilon = 1.0/(1ULL<<50);
+#endif
+
+/*! \brief The distanceEpsilon constant for physical distance
+/* comparisons
+ *
+ * \ingroup egspp_main
+ *
+ * The distanceEpsilon constant is a small number that can be used
+ * when comparing distances in the egs++ geometry library. It is
+ * independent from the machine epsilon defined above. In double
+ * precision it is set to 1.0/(1<<32), i.e., about 2.3283e-10, which
+ * in centimetres (default units) is smaller than any atomic radius.
+ * In single precision it is ~ 4.7684e-7 cm, in which case the
+ * geometrical precision limit is then roughly 5 nm.
+ *
+ * The distanceEpsilon constant is used as default value for the
+ * boundaryTolerance in the egs++ geometries. It is meant to be used
+ * as a float guard for comparisons involving distances that are on
+ * the order of 1.0 cm.
+ */
+#ifdef SINGLE
+ const EGS_Float distanceEpsilon = 1.0/(1<<21);
+#else
+ const EGS_Float distanceEpsilon = 1.0/(1<<32);
+#endif
 
 /*! \brief The maximum number of iterations for near-infinite loops
  *