diff --git a/NumericalConstants.cc b/NumericalConstants.cc
index eacb3434..de53e4a0 100644
--- a/NumericalConstants.cc
+++ b/NumericalConstants.cc
@@ -19,6 +19,7 @@
 
 #include <cstdlib>
 #include <cassert>
+#include <cmath>
 #include <iostream>
 
 #include "NumericalConstants.hh"
@@ -26,16 +27,21 @@
 int
 NumericalConstants::AddConstant(const string &iConst)
 {
-  map<string, int>::iterator iter = numConstantsIndex.find(iConst);
+  map<string, int>::const_iterator iter = numConstantsIndex.find(iConst);
 
   if (iter != numConstantsIndex.end())
     return iter->second;
 
-  assert(atof(iConst.c_str()) >= 0);
-
   int id = (int) mNumericalConstants.size();
   mNumericalConstants.push_back(iConst);
   numConstantsIndex[iConst] = id;
+
+  char *endptr;
+  double val = strtod(iConst.c_str(), &endptr);
+  assert(endptr != iConst.c_str()); // Check that the conversion succeeded
+  assert(val >= 0 || isnan(val)); // Check we have a positive constant or a NaN
+  double_vals.push_back(val);
+
   return id;
 }
 
@@ -47,7 +53,8 @@ NumericalConstants::get(int ID) const
 }
 
 double
-NumericalConstants::getDouble(int iID) const
+NumericalConstants::getDouble(int ID) const
 {
-  return(atof(get(iID).c_str()));
+  assert(ID >= 0 && ID < (int) double_vals.size());
+  return(double_vals[ID]);
 }
diff --git a/NumericalConstants.hh b/NumericalConstants.hh
index de771270..49adc774 100644
--- a/NumericalConstants.hh
+++ b/NumericalConstants.hh
@@ -32,15 +32,17 @@ class NumericalConstants
 private:
   //! Vector of numerical constants
   vector<string> mNumericalConstants;
+  //! Double values of these constants
+  vector<double> double_vals;
   //! Map matching constants to their id
   map<string, int> numConstantsIndex;
 public:
   //! Adds a constant and returns its ID
   int AddConstant(const string &iConst);
   //! Get a constant in string form
-  string get(int iID) const;
+  string get(int ID) const;
   //! Get a constant in double form
-  double getDouble(int iID) const;
+  double getDouble(int ID) const;
 };
 
 #endif