Added some more methods to Interpolate (constness of compute and getYi)

This commit is contained in:
Guilhem Lavaux 2012-06-04 22:21:06 -04:00
parent 6d410dc129
commit c5464dce31
3 changed files with 36 additions and 9 deletions

View File

@ -1,8 +1,9 @@
#include <fstream>
#include <iostream> #include <iostream>
#include "dinterpolate.hpp" #include "dinterpolate.hpp"
#define NX 10 #define NX 30
#define NY 10 #define NY 30
using namespace std; using namespace std;
using namespace CosmoTool; using namespace CosmoTool;
@ -11,18 +12,19 @@ typedef DelaunayInterpolate<double,double,2> myTriangle;
int main() int main()
{ {
myTriangle::CoordType pos[] = { {0,0}, {1,0}, {0,1}, {1, 1} } ; myTriangle::CoordType pos[] = { {0,0}, {1,0}, {0,1}, {1, 1}, {2, 0}, {2, 1} } ;
double vals[] = { 0, 1, 1, 0 }; double vals[] = { 0, 1, 1, 0, -0.5, 2.0 };
uint32_t simplex[] = { 0, 1, 2, 3, 1, 2 }; 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 iy = 0; iy <= NY; iy++) {
for (uint32_t ix = 0; ix <= NX; ix++) { for (uint32_t ix = 0; ix <= NX; ix++) {
myTriangle::CoordType inter = { ix *1.0/ NX, iy *1.0/NY }; myTriangle::CoordType inter = { ix *2.0/ NX, iy *1.0/NY };
cout << inter[1] << " " << inter[0] << " " << t.computeValue(inter) << endl; f << inter[1] << " " << inter[0] << " " << t.computeValue(inter) << endl;
} }
cout << endl; f << endl;
} }
return 0; return 0;
} }

View File

@ -49,6 +49,28 @@ double Interpolate::compute(double x)
return y; 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) double Interpolate::derivative(double x)
throw (InvalidRangeException) throw (InvalidRangeException)
{ {

View File

@ -21,6 +21,8 @@ namespace CosmoTool
double compute(double x) double compute(double x)
throw (InvalidRangeException); throw (InvalidRangeException);
double compute(double x) const
throw (InvalidRangeException);
double derivative(double x) double derivative(double x)
throw (InvalidRangeException); throw (InvalidRangeException);
@ -30,6 +32,7 @@ namespace CosmoTool
void fillWithXY(double *x, double *y) const; void fillWithXY(double *x, double *y) const;
double getMaxX() const; double getMaxX() const;
double getXi(int i) const { return spline->x[i]; } double getXi(int i) const { return spline->x[i]; }
double getYi(int i) const { return spline->y[i]; }
protected: protected:
gsl_interp_accel *accel_interp; gsl_interp_accel *accel_interp;
gsl_spline *spline; gsl_spline *spline;