From 5b3bef64b19bc6df105943b4d8b5023ae515a128 Mon Sep 17 00:00:00 2001
From: Your Name <guilhem.lavaux@gmail.com>
Date: Wed, 9 Feb 2011 14:32:01 -0500
Subject: [PATCH] Added a functor to interface with newton

---
 src/interpolate.cpp | 26 ++++++++++++++++++++++++++
 src/interpolate.hpp | 15 +++++++++++++++
 2 files changed, 41 insertions(+)

diff --git a/src/interpolate.cpp b/src/interpolate.cpp
index bf418d7..3693b4e 100644
--- a/src/interpolate.cpp
+++ b/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;
diff --git a/src/interpolate.hpp b/src/interpolate.hpp
index 08a9a88..ba839c4 100644
--- a/src/interpolate.hpp
+++ b/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;
@@ -41,6 +44,18 @@ namespace CosmoTool
   Interpolate buildInterpolateFromColumns(const char *fname, uint32_t col1, uint32_t col2, bool logx = false, bool logy = false)
     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;
+  };
   
 };