guilhem_lavaux/cosmotool
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Added support for log-log interpolation

Browse Source
This commit is contained in:
Guilhem Lavaux 2009-02-24 11:46:18 -06:00
parent 298acdfd11
commit 7b67b408ec
2 changed files with 24 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

27
src/interpolate.cpp
View File

@ -1,3 +1,4 @@
#include <cmath>
#include <algorithm> #include <algorithm>
#include "interpolate.hpp" #include "interpolate.hpp"
#include <vector> #include <vector>
@ -8,12 +9,15 @@
using namespace std; using namespace std;
using namespace CosmoTool; using namespace CosmoTool;
Interpolate::Interpolate(double *xs, double *values, uint32_t N, bool autofree) Interpolate::Interpolate(double *xs, double *values, uint32_t N, bool autofree,
bool logx, bool logy)
{ {
spline = gsl_spline_alloc (gsl_interp_linear, N); spline = gsl_spline_alloc (gsl_interp_linear, N);
gsl_spline_init (spline, xs, values, N); gsl_spline_init (spline, xs, values, N);
accel_interp = gsl_interp_accel_alloc(); accel_interp = gsl_interp_accel_alloc();
this->logx = logx;
this->logy = logy;
this->autoFree = autofree; this->autoFree = autofree;
} }
@ -29,10 +33,19 @@ double Interpolate::compute(double x)
throw (InvalidRangeException) throw (InvalidRangeException)
{ {
double y; double y;
if (logx)
x = log(x);
int err = gsl_spline_eval_e(spline, x, accel_interp, &y); int err = gsl_spline_eval_e(spline, x, accel_interp, &y);
if (err) if (err)
throw InvalidRangeException("Interpolate argument outside range"); throw InvalidRangeException("Interpolate argument outside range");
if (logy)
return exp(y);
else
return y;
} }
uint32_t Interpolate::getNumPoints() const uint32_t Interpolate::getNumPoints() const
@ -66,7 +79,7 @@ const Interpolate& Interpolate::operator=(const Interpolate& a)
double Interpolate::getMaxX() const double Interpolate::getMaxX() const
{ {
return spline->x[spline->size-1]; return exp(spline->x[spline->size-1]);
} }
typedef struct { typedef struct {
@ -136,7 +149,8 @@ Interpolate CosmoTool::buildInterpolateFromFile(const char *fname)
return inter; return inter;
} }
Interpolate CosmoTool::buildInterpolateFromColumns(const char *fname, uint32_t col1, uint32_t col2) Interpolate CosmoTool::buildInterpolateFromColumns(const char *fname, uint32_t col1, uint32_t col2, bool logx,
bool logy)
throw (NoSuchFileException,InvalidRangeException) throw (NoSuchFileException,InvalidRangeException)
{ {
vector<MyPair> allData; vector<MyPair> allData;
@ -198,15 +212,14 @@ Interpolate CosmoTool::buildInterpolateFromColumns(const char *fname, uint32_t c
for (uint32_t i = 0; i < allData.size(); i++) for (uint32_t i = 0; i < allData.size(); i++)
{ {
x[i] = allData[i].a; x[i] = logx ? log(allData[i].a) : allData[i].a;
y[i] = allData[i].b; y[i] = logy ? log(allData[i].b) : allData[i].b;
} }
Interpolate inter = Interpolate(x, y, allData.size(), true); Interpolate inter = Interpolate(x, y, allData.size(), true, logx, logy);
delete[] x; delete[] x;
delete[] y; delete[] y;
return inter; return inter;
} }

6
src/interpolate.hpp
View File

@ -15,7 +15,8 @@ namespace CosmoTool
{ {
public: public:
Interpolate() : spline(0), accel_interp(0) { } Interpolate() : spline(0), accel_interp(0) { }
Interpolate(double *xs, double *values, uint32_t N, bool autofree = false); Interpolate(double *xs, double *values, uint32_t N, bool autofree = false,
bool logx = false, bool logy = false);
~Interpolate(); ~Interpolate();
double compute(double x) double compute(double x)
@ -29,13 +30,14 @@ namespace CosmoTool
gsl_interp_accel *accel_interp; gsl_interp_accel *accel_interp;
gsl_spline *spline; gsl_spline *spline;
bool autoFree; bool autoFree;
bool logx, logy;
}; };
typedef std::vector< std::pair<double,double> > InterpolatePairs; typedef std::vector< std::pair<double,double> > InterpolatePairs;
Interpolate buildInterpolateFromFile(const char *fname) Interpolate buildInterpolateFromFile(const char *fname)
throw (NoSuchFileException); throw (NoSuchFileException);
Interpolate buildInterpolateFromColumns(const char *fname, uint32_t col1, uint32_t col2) Interpolate buildInterpolateFromColumns(const char *fname, uint32_t col1, uint32_t col2, bool logx = false, bool logy = false)
throw (NoSuchFileException,InvalidRangeException); throw (NoSuchFileException,InvalidRangeException);
Interpolate buildFromVector(const InterpolatePairs& v); Interpolate buildFromVector(const InterpolatePairs& v);