00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032 #include "GSLInterpolator.h"
00033
00034
00035 #include <cassert>
00036
00037 namespace ROOT {
00038 namespace Math {
00039
00040
00041 GSLInterpolator::GSLInterpolator (unsigned int size, Interpolation::Type type) :
00042 fAccel(0),
00043 fSpline(0)
00044 {
00045
00046
00047 switch ( type )
00048 {
00049 case ROOT::Math::Interpolation::kLINEAR :
00050 fInterpType = gsl_interp_linear;
00051 break ;
00052 case ROOT::Math::Interpolation::kPOLYNOMIAL :
00053 fInterpType = gsl_interp_polynomial;
00054 break ;
00055
00056 case ROOT::Math::Interpolation::kCSPLINE :
00057 fInterpType = gsl_interp_cspline ;
00058 break ;
00059 case ROOT::Math::Interpolation::kCSPLINE_PERIODIC :
00060 fInterpType = gsl_interp_cspline_periodic ;
00061 break ;
00062 case ROOT::Math::Interpolation::kAKIMA :
00063 fInterpType = gsl_interp_akima;
00064 break ;
00065 case ROOT::Math::Interpolation::kAKIMA_PERIODIC :
00066 fInterpType = gsl_interp_akima_periodic;
00067 break ;
00068 default :
00069
00070 fInterpType = gsl_interp_cspline;
00071 break ;
00072 }
00073
00074
00075 if (size >= fInterpType->min_size)
00076 fSpline = gsl_spline_alloc( fInterpType, size);
00077
00078 }
00079
00080 bool GSLInterpolator::Init(unsigned int size, const double *x, const double * y) {
00081
00082
00083 if (fSpline == 0)
00084 fSpline = gsl_spline_alloc( fInterpType, size);
00085
00086 else {
00087 gsl_interp * interp = fSpline->interp;
00088 if (size != interp->size) {
00089
00090 gsl_spline_free(fSpline);
00091 fSpline = gsl_spline_alloc( fInterpType, size);
00092
00093 }
00094 }
00095 if (!fSpline) return false;
00096
00097 int iret = gsl_spline_init( fSpline , x , y , size );
00098 if (iret != 0) return false;
00099
00100 fAccel = gsl_interp_accel_alloc() ;
00101
00102
00103
00104 assert (fSpline != 0);
00105 assert (fAccel != 0);
00106 return true;
00107 }
00108
00109 GSLInterpolator::~GSLInterpolator()
00110 {
00111
00112 if (fSpline != 0) gsl_spline_free(fSpline);
00113 if (fAccel != 0) gsl_interp_accel_free( fAccel);
00114 }
00115
00116 GSLInterpolator::GSLInterpolator(const GSLInterpolator &)
00117 {
00118
00119 }
00120
00121 GSLInterpolator & GSLInterpolator::operator = (const GSLInterpolator &rhs)
00122 {
00123
00124 if (this == &rhs) return *this;
00125
00126 return *this;
00127 }
00128
00129
00130 }
00131 }