From c5464dce3133fa4e7a3370422d29de1a992ce6ad Mon Sep 17 00:00:00 2001
From: Guilhem Lavaux <lavaux@uwaterloo.ca>
Date: Mon, 4 Jun 2012 22:21:06 -0400
Subject: [PATCH] Added some more methods to Interpolate (constness of compute
 and getYi)

---
 sample/testDelaunay.cpp | 20 +++++++++++---------
 src/interpolate.cpp     | 22 ++++++++++++++++++++++
 src/interpolate.hpp     |  3 +++
 3 files changed, 36 insertions(+), 9 deletions(-)

diff --git a/sample/testDelaunay.cpp b/sample/testDelaunay.cpp
index bc2f6a7..fe14553 100644
--- a/sample/testDelaunay.cpp
+++ b/sample/testDelaunay.cpp
@@ -1,8 +1,9 @@
+#include <fstream>
 #include <iostream>
 #include "dinterpolate.hpp"
 
-#define NX 10
-#define NY 10
+#define NX 30
+#define NY 30
 
 using namespace std;
 using namespace CosmoTool;
@@ -11,18 +12,19 @@ typedef DelaunayInterpolate<double,double,2> myTriangle;
 
 int main()
 {
-  myTriangle::CoordType pos[] = { {0,0}, {1,0}, {0,1}, {1, 1} } ;
-  double vals[] = { 0, 1, 1, 0 };
-  uint32_t simplex[] = { 0, 1, 2, 3, 1, 2 };
+  myTriangle::CoordType pos[] = { {0,0}, {1,0}, {0,1}, {1, 1}, {2, 0}, {2, 1} } ;
+  double vals[] = { 0, 1, 1, 0, -0.5, 2.0 };
+  uint32_t simplex[] = { 0, 1, 2, 3, 1, 2, 1, 3, 4, 4, 5, 3 };
 
-  myTriangle t(&pos[0], &vals[0], &simplex[0], 4, 2);
+  myTriangle t(&pos[0], &vals[0], &simplex[0], 6, 4);
+  ofstream f("output.txt");
 
   for (uint32_t iy = 0; iy <= NY; iy++) {
     for (uint32_t ix = 0; ix <= NX; ix++) {
-      myTriangle::CoordType inter = { ix *1.0/ NX, iy *1.0/NY };
-      cout << inter[1] << " " << inter[0] << " " << t.computeValue(inter) << endl;
+      myTriangle::CoordType inter = { ix *2.0/ NX, iy *1.0/NY };
+      f << inter[1] << " " << inter[0] << " " << t.computeValue(inter) << endl;
     }
-    cout << endl;
+    f << endl;
   }
   return 0;
 }
diff --git a/src/interpolate.cpp b/src/interpolate.cpp
index d9e3e4f..9a9fc70 100644
--- a/src/interpolate.cpp
+++ b/src/interpolate.cpp
@@ -49,6 +49,28 @@ double Interpolate::compute(double x)
     return y;
 }
 
+double Interpolate::compute(double x) const
+  throw (InvalidRangeException)
+{
+  double y;
+
+  if (logx)
+    x = log(x);
+
+  int err = gsl_spline_eval_e(spline, x, 0, &y);
+
+  if (err)
+    throw InvalidRangeException("Interpolate argument outside range");
+
+  if (logy)
+    return exp(y);
+  else
+    return y;
+}
+
+
+
+
 double Interpolate::derivative(double x)
   throw (InvalidRangeException)
 {
diff --git a/src/interpolate.hpp b/src/interpolate.hpp
index ba839c4..a0ffacc 100644
--- a/src/interpolate.hpp
+++ b/src/interpolate.hpp
@@ -21,6 +21,8 @@ namespace CosmoTool
 
     double compute(double x)
       throw (InvalidRangeException);
+    double compute(double x) const
+      throw (InvalidRangeException);
     double derivative(double x)
       throw (InvalidRangeException);
 
@@ -30,6 +32,7 @@ namespace CosmoTool
     void fillWithXY(double *x, double *y) const;
     double getMaxX() const;
     double getXi(int i) const { return spline->x[i]; }
+    double getYi(int i) const { return spline->y[i]; }
   protected:
     gsl_interp_accel *accel_interp;
     gsl_spline *spline;