Interp1D.H

Go to the documentation of this file.

 
#ifndef _INTERP1D_H_
#define _INTERP1D_H_
 
#include <cmath>
#include <gsl/gsl_interp.h>
#include "ErrCodes.H"
#include "ThreadVec.H"
#include "Types.H"
 
namespace criptic {
 
  typedef enum {
    extrapError,   
    extrapLinear,  
    extrapConst,   
    extrapFixedVal 
  } extrapType;
  
  class Interp1D {
 
  public:
 
    Interp1D(const double* x_, const double* y_, const size_t n_,
         const extrapType extrapLo_ = extrapError,
         const extrapType extrapHi_ = extrapError,
         const bool xLog_ = false, const bool yLog_ = false,
         const double yLo_ = 0.0, const double yHi_ = 0.0) :
      x(x_),
      y(y_),
      n(n_),
      extrapLo(extrapLo_),
      extrapHi(extrapHi_),
      xLog(xLog_),
      yLog(yLog_),
      yLo(yLo_),
      yHi(yHi_)
    {
 
      // Allocate and initialize GSL objects
      interp = gsl_interp_alloc(gsl_interp_steffen, n);
      gsl_interp_init(interp, x, y, n);
      for (IdxType i = 0; i < acc.size(); i++)
    acc[i] = gsl_interp_accel_alloc();
 
      // Store slopes
      if (extrapLo == extrapLinear)
    mLo = (y[1] - y[0]) / (x[1] - x[0]);
      if (extrapHi == extrapLinear)
    mHi = (y[n-1] - y[n-2]) / (x[n-1] - x[n-2]);
    }      
 
    ~Interp1D() {
      // Free memory
      gsl_interp_free(interp);
      for (IdxType i = 0; i < acc.size(); i++)
    gsl_interp_accel_free(acc[i]);
    }
 
    double xMin() const { return xLog ? std::pow(10., x[0]) : x[0]; }
    
    double xMax() const { return xLog ? std::pow(10., x[n-1]) : x[n-1]; }
 
    double operator()(const Real& x_) const {
 
      // Take log if requested, and catch case where taking the log
      // and then its inverse causes a value of x that should be at
      // the edge of the table to move slighlty off it; to handle
      // this, if xLog is set, we check for near-equality between x1
      // and the table limits, and coerce x1 back onto the table if it
      // is off by less than some tolerance.
      double x1 = xLog ? std::log10(x_) : x_;
      if (xLog) {
    if (x1 < x[0] && x[0] - x1 < 1.0e-10) x1 = x[0];
    else if (!(x1 < x[n-1]) && x1 - x[n-1] < tol) x1 = x[n-1];
      }
 
      // Compute interpolation
      double y1 = 0.0;
      if (x1 < x[0]) {
    // Off table low
    switch (extrapLo) {
    case extrapLinear: {
      y1 = y[0] + mLo * (x1 - x[0]);
      break;
    }
    case extrapFixedVal: {
      return yLo;
    }
    case extrapConst: {
      y1 = y[0];
      break;
    }
    case extrapError: {
      throw(errBadXLoExtrapolation);
    }
    }
      } else if (x1 > x[n-1]) {
    // Off table high
    switch (extrapHi) {
    case extrapLinear: {
      y1 = y[n-1] + mHi * (x1 - x[n-1]);
      break;
    }
    case extrapFixedVal: {
      return yHi;
    }
    case extrapConst: {
      y1 = y[n-1];
      break;
    }
    case extrapError: {
      throw(errBadXHiExtrapolation);
    }
    }   
      } else {
    // Within table
    y1 = gsl_interp_eval(interp, x, y, x1, acc());
      }
      if (yLog) {
    return std::pow(10., y1);
      } else {
    return y1;
      }
    }
 
  private:
 
    const double *x;     
    const double *y;     
    const size_t n;      
    const extrapType extrapLo; 
    const extrapType extrapHi; 
    const bool xLog;     
    const bool yLog;     
    const double yLo;    
    const double yHi;    
    double mLo;          
    double mHi;          
    static constexpr double tol = 1.0e-8;  
    // GSL machinery
    gsl_interp *interp;                
    ThreadVec<gsl_interp_accel *> acc; 
  };
 
}
 
#endif
// _INTERP1D_H_