Interp2D.H

Go to the documentation of this file.

 
#ifndef _INTERP2D_H_
#define _INTERP2D_H_
 
#include <cmath>
#include <gsl/gsl_interp2d.h>
#include "ErrCodes.H"
#include "Interp1D.H"
#include "ThreadVec.H"
#include "Types.H"
 
namespace criptic {
 
  class Interp2D {
 
  public:
 
    Interp2D(const double* x_, const double* y_, const double* z_,
         const size_t nx_, const size_t ny_,
         const extrapType extrapXLo_ = extrapError,
         const extrapType extrapXHi_ = extrapError,
         const extrapType extrapYLo_ = extrapError,
         const extrapType extrapYHi_ = extrapError,
         const bool xLog_ = false, const bool yLog_ = false,
         const bool zLog_ = false,
         const double zXLo_ = 0.0, const double zXHi_ = 0.0,
         const double zYLo_ = 0.0, const double zYHi_ = 0.0) :
      x(x_),
      y(y_),
      z(z_),
      nx(nx_),
      ny(ny_),
      extrapXLo(extrapXLo_),
      extrapXHi(extrapXHi_),
      extrapYLo(extrapYLo_),
      extrapYHi(extrapYHi_),
      xLog(xLog_),
      yLog(yLog_),
      zLog(zLog_),
      zXLo(zXLo_),
      zXHi(zXHi_),
      zYLo(zYLo_),
      zYHi(zYHi_)
    {
 
      // Allocate and initialize GSL objects
      interp = gsl_interp2d_alloc(gsl_interp2d_bilinear, nx, ny);
      gsl_interp2d_init(interp, x, y, z, nx, ny);
      for (IdxType i = 0; i < accX.size(); i++) {
    accX[i] = gsl_interp_accel_alloc();
    accY[i] = gsl_interp_accel_alloc();
      }
    }      
 
    ~Interp2D() {
      // Free memory
      gsl_interp2d_free(interp);
      for (IdxType i = 0; i < accX.size(); i++) {
    gsl_interp_accel_free(accX[i]);
    gsl_interp_accel_free(accY[i]);
      }
    }
 
    double xMin() const { return xLog ? std::pow(10., x[0]) : x[0];  }
    
    double xMax() const { return xLog ? std::pow(10., x[nx-1]) : x[nx-1];  }
 
    double yMin() const { return yLog ? std::pow(10., y[0]) : y[0]; }
    
    double yMax() const { return yLog ? std::pow(10., y[ny-1]) : y[ny-1]; }
    
    double operator()(const Real& x_, const Real& y_) const {
 
      // Take log if requested
      double x1 = xLog ? std::log10(x_) : x_;
      double y1 = yLog ? std::log10(y_) : y_;
 
      // Catch case where we are logging variables and we input a
      // value that is at the top or bottom of the table, but taking
      // the logarithm leads us to be off the table by a machine
      // precision-level amount; this can cause an error if
      // extrapolation is forbidden. To handle this, if xLog or yLog
      // is set, we check for near-equality between x1 and y1 and the
      // extrema of the table, and coerce the points back onto the
      // table if they are off by a very small amount.
      if (xLog) {
    if (x1 < x[0] && x[0] - x1 < 1.0e-10) x1 = x[0];
    else if (!(x1 < x[nx-1]) && x1 - x[nx-1] < tol) x1 = x[nx-1];
      }
      if (yLog) {
    if (y1 < y[0] && y[0] - y1 < 1.0e-10) y1 = y[0];
    else if (!(y1 < y[ny-1]) && y1 - y[ny-1] < tol) y1 = y[ny-1];
      }
 
      // Find out where we are relative to grid; note that there are 9
      // possible regions, defined by being on the grid and in 8
      // possible zones around it
      double z1 = 0.0;
      if (x1 < x[0]) {  
    
    // Off-table low in x direction; bail if this region is
    // disallowed
    if (extrapXLo == extrapError)
      throw(errBadXLoExtrapolation);
 
    // Now check with region we are in in the y direction
    if (y1 < y[0]) {
 
      // Allowed combinations in this region are Const for both x
      // and y, in which case we just return the (0,0) value,
      // FixedVal for both *if the fixed values are the same*, and
      // linear for both, in which case we linearly extrapolate in
      // both directions; any other combination is an error
      if (extrapXLo == extrapConst &&
          extrapYLo == extrapConst) {
        z1 = z[0];
      } else if (extrapXLo == extrapFixedVal &&
             extrapYLo == extrapFixedVal &&
             zXLo == zYLo) {
        z1 = zXLo;
      } else if (extrapXLo == extrapLinear &&
             extrapYLo == extrapLinear) {
        double mx = (z[1] - z[0]) / (x[1] - x[0]);
        double my = (z[nx] - z[0]) / (y[1] - y[0]);
        z1 = z[0] + mx * (x1 - x[0]) + my * (y1 - y[0]);
      } else {
        throw(errBadXYExtrapolation);
      }
      
    } else if (y1 < y[ny-1]) {
 
      // On-table in y direction
      switch (extrapXLo) {
      case extrapLinear: {
        double zx0 = gsl_interp2d_eval(interp, x, y, z,
                       x[0], y1, accX(), accY());
        double zx1 = gsl_interp2d_eval(interp, x, y, z,
                       x[1], y1, accX(), accY());
        double mx = (zx1 - zx0) / (x[1] - x[0]);
        z1 = zx0 + mx * (x1 - x[0]);
        break;
      }
      case extrapFixedVal:
        return zXLo;
      case extrapConst: {
        z1 = gsl_interp2d_eval(interp, x, y, z,
                   x[0], y1, accX(), accY());
        break;
      }
      case extrapError:
        throw(errBadXLoExtrapolation);
      }
 
    } else {
 
      // Off-table high in y direction; bail if disallowed
      if (extrapYHi == extrapError) {
        throw(errBadYHiExtrapolation);
      }
      
      // Allowed combinations in this region are Const for both x
      // and y, in which case we just return the (0,ny-1) value,
      // FixedVal for both *if the fixed values are the same*, and
      // linear for both, in which case we linearly extrapolate in
      // both directions; any other combination is an error
      if (extrapXLo == extrapConst &&
          extrapYHi == extrapConst) {
        z1 = z[(ny-1)*nx];
      } else if (extrapXLo == extrapFixedVal &&
             extrapYHi == extrapFixedVal &&
             zXLo == zYHi) {
        z1 = zXLo;
      } else if (extrapXLo == extrapLinear &&
             extrapYHi == extrapLinear) {
        double mx = (z[ (ny-1) * nx + 1] - z[ (ny-1) * nx ]) /
          (x[1] - x[0]);
        double my = (z[ (ny-1) * nx ] - z[ (ny-2) * nx ])
          / (y[ny-1] - y[ny-2]);
        z1 = z[ (ny-1) * nx ] + mx * (x1 - x[0]) + my * (y1 - y[ny-1]);
      } else {
        throw(errBadXYExtrapolation);
      }
 
    }
 
      } else if (x1 < x[nx-1]) {
 
    // On-table in x direction; see where we are in y
    if (y1 < y[0]) {
 
      // Off-table low in y
      switch (extrapYLo) {
      case extrapLinear: {
        double zy0 = gsl_interp2d_eval(interp, x, y, z,
                       x1, y[0], accX(), accY());
        double zy1 = gsl_interp2d_eval(interp, x, y, z,
                       x1, y[1], accX(), accY());
        double my = (zy1 - zy0) / (y[1] - y[0]);
        z1 = zy0 + my * (y1 - y[0]);
        break;
      }
      case extrapFixedVal:
        return zYLo;
      case extrapConst: {
        z1 = gsl_interp2d_eval(interp, x, y, z,
                   x1, y[0], accX(), accY());
        break;
      }
      case extrapError: {
        throw(errBadYLoExtrapolation);
      }
      }
 
    } else if (y1 < y[ny-1]) {
 
      // On table in both directions
      z1 = gsl_interp2d_eval(interp, x, y, z,
                 x1, y1, accX(), accY());
 
    } else {
 
      // Off-table high in y direction
      switch (extrapYHi) {
      case extrapLinear: {
        double zy0 = gsl_interp2d_eval(interp, x, y, z,
                       x1, y[ny-2], accX(), accY());
        double zy1 = gsl_interp2d_eval(interp, x, y, z,
                       x1, y[ny-1], accX(), accY());
        double my = (zy1 - zy0) / (y[ny-1] - y[ny-2]);
        z1 = zy1 + my * (y1 - y[0]);
        break;
      }
      case extrapFixedVal:
        return zYHi;
      case extrapConst: {
        z1 = gsl_interp2d_eval(interp, x, y, z,
                   x1, y[ny-1], accX(), accY());
        break;
      }
      case extrapError: {
        throw(errBadYHiExtrapolation);
      }
      }
    }
 
      } else {
 
    // Off-table high in x direction; bail if this region is
    // disallowed
    if (extrapXHi == extrapError)
      throw(errBadXHiExtrapolation);
 
    // Now check with region we are in in the y direction
    if (y1 < y[0]) {
 
      // Off-table low in y direction; bail if disallowed
      if (extrapYLo == extrapError)
        throw(errBadYLoExtrapolation);
      
      // Allowed combinations in this region are Const for both x
      // and y, in which case we just return the (nx-1,0) value,
      // FixedVal for both *if the fixed values are the same*, and
      // linear for both, in which case we linearly extrapolate in
      // both directions; any other combination is an error
      if (extrapXHi == extrapConst &&
          extrapYLo == extrapConst) {
        z1 = z[nx-1];
      } else if (extrapXHi == extrapFixedVal &&
             extrapYLo == extrapFixedVal &&
             zXHi == zYLo) {
        z1 = zXHi;
      } else if (extrapXHi == extrapLinear &&
             extrapYLo == extrapLinear) {
        double mx = (z[nx-1] - z[nx-2]) / (x[nx-1] - x[nx-2]);
        double my = (z[nx-1 + nx] - z[nx-1]) / (y[1] - y[0]);
        z1 = z[nx-1] + mx * (x1 - x[nx-1]) + my * (y1 - y[0]);
      } else {
        throw(errBadXYExtrapolation);
      }
      
    } else if (y1 < y[ny-1]) {
 
      // On-table in y direction
      switch (extrapXHi) {
      case extrapLinear: {
        double zx0 = gsl_interp2d_eval(interp, x, y, z,
                       x[nx-2], y1, accX(), accY());
        double zx1 = gsl_interp2d_eval(interp, x, y, z,
                       x[nx-1], y1, accX(), accY());
        double mx = (zx1 - zx0) / (x[nx-1] - x[nx-2]);
        z1 = zx1 + mx * (x1 - x[nx-1]);
        break;
      }
      case extrapFixedVal:
        return zXHi;
      case extrapConst: {
        z1 = gsl_interp2d_eval(interp, x, y, z,
                   x[nx-1], y1, accX(), accY());
        break;
      }
      case extrapError:
        throw(errBadXHiExtrapolation);
      }
 
    } else {
 
      // Off-table high in y direction; bail if disallowed
      if (extrapYHi == extrapError) {
        throw(errBadYHiExtrapolation);
      }
      
      // Allowed combinations in this region are Const for both x
      // and y, in which case we just return the (0,ny-1) value,
      // FixedVal for both *if the fixed values are the same*, and
      // linear for both, in which case we linearly extrapolate in
      // both directions; any other combination is an error
      if (extrapXHi == extrapConst &&
          extrapYHi == extrapConst) {
        z1 = z[(ny-1)*nx + nx-1];
      } else if (extrapXHi == extrapFixedVal &&
             extrapYHi == extrapFixedVal &&
             zXHi == zYHi) {
        z1 = zXHi;
      } else if (extrapXHi == extrapLinear &&
             extrapYHi == extrapLinear) {
        double mx = (z[ (ny-1) * nx + nx-2] - z[ (ny-1) * nx + nx-1 ]) /
          (x[nx-1] - x[nx-2]);
        double my = (z[ (ny-1) * nx + nx-1 ] - z[ (ny-2) * nx + nx-1])
          / (y[ny-1] - y[ny-2]);
        z1 = z[ (ny-1) * nx + nx-1] +
          mx * (x1 - x[nx-1]) + my * (y1 - y[ny-1]);
      } else {
        throw(errBadXYExtrapolation);
      }
 
    }
    
      }
 
      // Return
      if (zLog) {
    return std::pow(10., z1);
      } else {
    return z1;
      }
    }
 
  private:
 
    const double *x;     
    const double *y;     
    const double *z;     
    const size_t nx;     
    const size_t ny;     
    const extrapType extrapXLo; 
    const extrapType extrapXHi; 
    const extrapType extrapYLo; 
    const extrapType extrapYHi; 
    const bool xLog;     
    const bool yLog;     
    const bool zLog;     
    const double zXLo;    
    const double zXHi;    
    const double zYLo;    
    const double zYHi;    
    static constexpr double tol = 1.0e-8;  
    // GSL machinery
    gsl_interp2d *interp;               
    ThreadVec<gsl_interp_accel *> accX; 
    ThreadVec<gsl_interp_accel *> accY; 
  };
 
}
 
#endif
// _INTERP2D_H_