Added a functor to interface with newton

src/interpolate.cpp

@@ -1,3 +1,4 @@
+#include <cassert>
 #include <cmath>
 #include <algorithm>
 #include "interpolate.hpp"
@@ -48,6 +49,31 @@ double Interpolate::compute(double x)
     return y;
 }
 
+double Interpolate::derivative(double x)
+  throw (InvalidRangeException)
+{
+  double y, dy, x0 = x;
+  
+  if (logx)
+    x0 = log(x0);
+
+  int err = gsl_spline_eval_deriv_e(spline, x0, accel_interp, &dy);
+
+  if (err)
+    throw InvalidRangeException("Interpolate argument outside range");
+
+  if (logy)
+    {
+      int err = gsl_spline_eval_e(spline, x0, accel_interp, &y);
+      
+      assert(err == 0);
+      
+      return dy*exp(y)/x0;
+    }
+  else
+    return dy;
+}
+
 uint32_t Interpolate::getNumPoints() const
 {
   return spline->size;

src/interpolate.hpp

@@ -21,6 +21,9 @@ namespace CosmoTool
 
     double compute(double x)
       throw (InvalidRangeException);
+    double derivative(double x)
+      throw (InvalidRangeException);
+
     const Interpolate& operator=(const Interpolate& a);
 
     uint32_t getNumPoints() const;
@@ -42,6 +45,18 @@ namespace CosmoTool
     throw (NoSuchFileException,InvalidRangeException);
   Interpolate buildFromVector(const InterpolatePairs& v);
 
+
+  class FunctorInterpolate
+  {
+  public:
+    FunctorInterpolate(Interpolate& i) : m_i(i) {}    
+
+    double eval(double x) { return m_i.compute(x); }
+    double derivative(double x) { return m_i.derivative(x); }
+  private:
+    Interpolate& m_i;
+  };
+  
 };
 
 #endif