NumericalApproximation.h

#ifndef NumericalApproximation_h
#define NumericalApproximation_h

#include <string>
#include <vector>
#include <map>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include "IException.h"

using namespace std;

namespace Isis {
  class NumericalApproximation {
    public:
      enum InterpType { Linear,               
                        Polynomial,           
                        PolynomialNeville,    
                        CubicNatural,         
                        CubicClamped,         
                        CubicNatPeriodic,     
                        CubicNeighborhood,    
                        CubicHermite,         
                        Akima,                
                        AkimaPeriodic         
                      };

      // CONSTRUCTORS
      NumericalApproximation(const NumericalApproximation::InterpType &itype = CubicNatural);
      NumericalApproximation(unsigned int n, double *x, double *y,
                             const NumericalApproximation::InterpType &itype = CubicNatural);
      NumericalApproximation(const vector <double> &x, const vector <double> &y,
                             const NumericalApproximation::InterpType &itype = CubicNatural);
      NumericalApproximation(const NumericalApproximation &dint);
      // ASSIGNMENT OPERATOR
      NumericalApproximation &operator=(const NumericalApproximation &numApMeth);
      // DESTRUCTOR
      virtual ~NumericalApproximation();


      // ACCESSOR METHODS FOR OBJECT PROPERTIES
      string Name() const;
      inline InterpType   InterpolationType() {
        return (p_itype);
      }
      int MinPoints();
      int MinPoints(NumericalApproximation::InterpType itype);

      // ADD DATA TO OBJECT
      void AddData(const double x, const double y);
      void AddData(unsigned int n, double *x, double *y);
      void AddData(const vector <double> &x, const vector <double> &y);
      void SetCubicClampedEndptDeriv(const double yp1, const double ypn);
      void AddCubicHermiteDeriv(unsigned int n, double *fprimeOfx);
      void AddCubicHermiteDeriv(const vector <double> &fprimeOfx);
      void AddCubicHermiteDeriv(const double fprimeOfx);

      //ACCESSOR METHODS AFTER DATA IS ENTERED
      double DomainMinimum();
      double DomainMaximum();
      bool Contains(double x);
      inline unsigned int Size() {
        return(p_x.size());
      }

      enum ExtrapType { ThrowError,     
                        Extrapolate,    
                        NearestEndpoint 
                      };
      // INTERPOLATION AND EXTRAPOLATION RESULTS
      double          Evaluate(const double          a, const ExtrapType &etype = ThrowError);
      vector <double> Evaluate(const vector <double> &a, const ExtrapType &etype = ThrowError);
      vector <double> PolynomialNevilleErrorEstimate();
      vector <double> CubicClampedSecondDerivatives();
      double EvaluateCubicHermiteFirstDeriv(const double a);
      double EvaluateCubicHermiteSecDeriv(const double a);

      // DIFFERENTIATION METHODS
      double GslFirstDerivative(const double a);
      double BackwardFirstDifference(const double a, const unsigned int n = 3,
                                     const double h = 0.1);
      double ForwardFirstDifference(const double a, const unsigned int n = 3,
                                    const double h = 0.1);
      double CenterFirstDifference(const double a, const unsigned int n = 5,
                                   const double h = 0.1);

      double GslSecondDerivative(const double a);
      double BackwardSecondDifference(const double a, const unsigned int n = 3,
                                      const double h = 0.1);
      double ForwardSecondDifference(const double a, const unsigned int n = 3,
                                     const double h = 0.1);
      double CenterSecondDifference(const double a, const unsigned int n = 5,
                                    const double h = 0.1);

      // INTERGRATION METHODS
      double GslIntegral(const double a, const double b);
      double TrapezoidalRule(const double a, const double b);
      double Simpsons3PointRule(const double a, const double b);
      double Simpsons4PointRule(const double a, const double b);
      double BoolesRule(const double a, const double b);
      double RefineExtendedTrap(double a, double b, double s, unsigned int n);
      double RombergsMethod(double a, double b);

      // ASSIGNMENT OPERATORS
      void Reset();
      void Reset(NumericalApproximation::InterpType itype);
      void SetInterpType(NumericalApproximation::InterpType itype);

    protected:
      // == CLASS VARIABLES
      // VARIABLES FOR ALL INTERP TYPES
      InterpType        p_itype;                            
      vector<double>    p_x;                                
      vector<double>    p_y;                                
      bool              p_dataValidated;                    
      // GSL INTERP VARIABLES
      typedef const gsl_interp_type *InterpFunctor;         
      typedef map<InterpType, InterpFunctor> FunctorList;   
      typedef FunctorList::const_iterator FunctorConstIter; 
      gsl_interp_accel *p_acc;                              
      gsl_spline       *p_interp;                           
      static FunctorList p_interpFunctors;                   
      // CUBIC CLAMPED VARIABLES
      bool              p_clampedEndptsSet;                 
      bool              p_clampedComputed;                  
      double            p_clampedDerivFirstPt;              
      double            p_clampedDerivLastPt;               
      vector<double>    p_clampedSecondDerivs;              
      // POLYNOMIAL NEVILLE VARIABLES
      vector <double>   p_polyNevError;                     
      // CUBIC HERMITE VARIABLES
      vector<double>    p_fprimeOfx;                        


      // == PROTECTED METHODS
      // CREATING, DESTROYING OBJECT
      void Init(NumericalApproximation::InterpType itype);
      bool GslInterpType(NumericalApproximation::InterpType itype) const;
      void GslAllocation(unsigned int npoints);
      void GslDeallocation();
      InterpFunctor GslFunctor(NumericalApproximation::InterpType itype) const;
      // VERIFICATION METHODS
      void GslIntegrityCheck(int gsl_status, const char *src,
                             int line);
      void ValidateDataSet();
      bool InsideDomain(const double a);
      // COMPUTATION AND EVALUATION METHODS
      bool   GslComputed() const;
      void   ComputeGsl();
      void   ComputeCubicClamped();
      double ValueToExtrapolate(const double a, const ExtrapType &etype);
      double          EvaluateCubicNeighborhood(const double a);
      vector <double> EvaluateCubicNeighborhood(const vector <double> &a, const ExtrapType &etype);
      double          EvaluateCubicClamped(const double a);
      double          EvaluateCubicHermite(const double a);
      double          EvaluatePolynomialNeville(const double a);
      vector <double> EvaluateForIntegration(const double a, const double b,
                                             const unsigned int n);
      int FindIntervalLowerIndex(const double a);
      // STANDARDIZE ERRORS
      void ReportException(IException::ErrorType type, const string &method,
                           const string &message, const char *filesrc,
                           int lineno) const;
  };
};

#endif