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); 

//              cout << "LinearIpol @ "<< t << " of " << size << endl;
//              printf("LinearIpol():%lg of %d\n", t, int(size) ); 

                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]); 

//              cout << "   LinearIpol: [" << i << "]=" << Data[i] << ", [" << i+1 << "]="<<Data[i+1]<<" @ " << t << endl; 

//              printf("LinearIpol(%d,%d)@%lg\n",i,i+1, t);
                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 */