LinearIpol.hpp

#ifndef __IPOL_LINEAR_HPP
#define __IPOL_LINEAR_HPP
 
//#include <iostream>
#include "IpolDelimiter.hpp"
#include <assert.h>
 
namespace Fiber
{
//using namespace std;
 
template <class T>
struct  LinearIpolZeroDerivativeTrait
{
static  T zero()
        {
                return T({0});
        }
};
 
template <class T, bool Vectorial = false>
class   LinearIpol
{
public:
        typedef T value_type;
 
        template <class Storage1D, class Limiter>
static  T interpolate(const Storage1D&Data, double t, index_t size, const Limiter&L)
        {
                assert(size > 0);
#ifdef  DEBUG_LINEAR_IPOL
                cout << "LinearIpol @ "<< t << " of " << size << endl;
                printf("LinearIpol():%lg of %d\n", t, int(size) );
#endif
                if (t<=0      ) return L.limit(Data[ 0 ]   );
                if (t+1>=size ) return L.limit(Data[size-1]);
 
        index_t i = index_t(t);
                t = t - floor(t);
 
                if (t==0.0)
                        return L.limit(Data[i]); 
 
#ifdef  DEBUG_LINEAR_IPOL
                cout << "   LinearIpol: [" << i << "]=" << Data[i] << ", [" << i+1 << "]="<<Data[i+1]<<" @ " << t << endl;
                printf("LinearIpol(%d,%d)@%lg\n",i,i+1, t); 
#endif 
 
                return L.limit(Data[i]) * (1-t) + L.limit(Data[i+1]) * t;
        }
 
        template <class Storage1D, class Limiter>
static  T derivative(const Storage1D&Data, double t, index_t size, const Limiter&L)
        {
//              assert(size > 1);
 
// need to think about return value here, since T might not always have a constructor with 0 value
// maybe need a type trait here
                if (t<0       ) return LinearIpolZeroDerivativeTrait<T>::zero();
                if (t>=size-1 ) return LinearIpolZeroDerivativeTrait<T>::zero();
 
        index_t i = index_t(t);
//              printf("LinearIpol Derivative(%d,%d)@%lg\n",i,i+1, t);
                return L.limit(Data[i+1]) - L.limit(Data[i]);
        }
};
 
 
template <class T>
class   LinearIpol<T, true>
{
public:
        typedef T value_type;
 
        template <class Storage1D, class Limiter>
static  T interpolate(const Storage1D&Data, double t, index_t size, const Limiter&L)
        {
                assert(size > 0);
 
                if (t<=0      ) return L.limit(Data[ 0 ]   );
                if (t+1>=size ) return L.limit(Data[size-1]);
 
        index_t i = index_t(t);
                t = t - floor(t);
 
                if (t==0.0)
                        return L.limit(Data[i]);
 
                return L.limit(Data[i]) + (L.limit(Data[i+1]) - L.limit(Data[i])) * t;
        }
 
        template <class Storage1D, class Limiter>
static  T derivative(const Storage1D&Data, double t, index_t size, const Limiter&L)
        {
//              assert(size > 1);
 
                //if (t<0       ) return T(0);
                //if (t>=size-1 ) return T(0);
 
        index_t i = index_t(t);
        //      printf("LinearIpol Derivative(%d,%d)@%lg\n",i,i+1, t);
                return L.limit(Data[i+1]) - L.limit(Data[i]);
        }
};
 
} /* namespace Fiber */
 
#endif /* __IPOL_LINEAR_HPP */