Interpolate.hpp

//
// $Id: Interpolate.hpp,v 1.5 2007/12/28 13:50:49 werner Exp $
//
#ifndef __FIBER_ARRAY_INTERPOLATE_HPP
#define __FIBER_ARRAY_INTERPOLATE_HPP "Created 3.07.2004 21:42:21 by bzfbenge"

#include <iostream>

#include "VVector.hpp"
#include "MultiArray.hpp"
#include "IpolDelimiter.hpp"

namespace Fiber
{

        using namespace std;

template <class T, class Interpol1D, bool Derive=false>
class   IpolDerivative
{
public:
        template <class Storage1D, class Limiter>
static  T interpolate(const Storage1D&Data, double t, index_t size, const Limiter&L)
        {
//              printf("IpolDerivative, non ipol, at %.20lg\n", t);
                return Interpol1D::interpolate(Data,t,size,L);
        }
};

template <class T, class Interpol1D>
class   IpolDerivative<T, Interpol1D, true>
{
public:
        template <class Storage1D, class Limiter>
static  T interpolate(const Storage1D&Data, double t, index_t size, const Limiter&L)
        {
                return Interpol1D::derivative(Data,t,size,L);
        }
};


template <int N, class Type, class Interpol, class CoordinateType = double,
                 class Delimiter = NoDelimiter<Type>, int DerivativeDimension=-1 >
struct Interpolate;

template <class Type, class Interpol, class CoordinateType, class Delimiter, int DerivativeDimension>
struct Interpolate<1, Type, Interpol, CoordinateType, Delimiter, DerivativeDimension>
{
        enum { Dims = 1 };

        typedef Type value_type;
        typedef MultiArray<Dims, Type> MultiArray_t;

        typedef FixedArray<CoordinateType,Dims> Point_t;

private:
        const MultiArray_t&M;
        const Point_t     &Point;

public: Delimiter myDelimiter;

        Interpolate(const MultiArray_t&m,
                    const Point_t&point)
        : M(m), Point(point)
        {}

        Interpolate(const MultiArray_t&m,
                    const Point_t&point,
                    const Delimiter&myDeli)
        : M(m), Point(point), myDelimiter(myDeli)
        {}

        Type operator[](index_t i) const
        {
        //      printf("\nInterpolate<1D, DerivativeDimension=%d>: Slice %d\n", DerivativeDimension, i);
#if 0
                Type T = IpolDerivative<Type, Interpol, 0==DerivativeDimension>::
			 interpolate(M, Point[0],  M.Extension()[0], myDelimiter ); 
                cout << "Interpolate<1D> at " << Point << " yields " << T << endl;
                return T;
#else
                return IpolDerivative<Type, Interpol, 0==DerivativeDimension>::
			interpolate(M, Point[0],  M.Extension()[0], myDelimiter );
#endif
        }

        index_t size() const
        {
                return M.Extension().maxidx();
        }


        Type eval() const
        {
        //      printf("\nInterpolate<DerivativeDimension=%d> Eval\n", DerivativeDimension);

                return IpolDerivative<Type, Interpol, 0==DerivativeDimension>::
			interpolate(M, Point[0],  M.Extension()[0], myDelimiter );
        }
};


template <class Type, class Interpol, class CoordinateType, class Delimiter, int DerivativeDimension>
struct Interpolate<2, Type, Interpol, CoordinateType, Delimiter, DerivativeDimension>
{
        enum { Dims = 2 };

        typedef Type value_type;
        typedef MultiArray<Dims, Type> MultiArray_t;
        typedef FixedArray<CoordinateType, Dims> Point_t;
        typedef typename MultiArray_t::SliceStorage_t SliceStorage_t;

private:
        const MultiArray_t &M;
        const Point_t      &Point;

public: Delimiter myDelimiter;

        Interpolate(const MultiArray_t&m,
                    const Point_t     &point)
        : M(m), Point(point)
        {}

        Interpolate(const MultiArray_t&m,
                    const Point_t     &point,
                    const Delimiter&myDeli)
        : M(m), Point(point), myDelimiter(myDeli)
        {}

        Type operator[](index_t i) const
        {
        MultiArray<1, Type> Sub = M[i]; 

        //      printf("\nInterpolate<2D, DerivativeDimension=%d, Type=%s>: slice %d, at %lg\n", 
        //             DerivativeDimension, typeid(Type).name(), i, double(Point[0]) );
#if 0
                Type T = IpolDerivative<Type, Interpol, 0==DerivativeDimension>::
			 interpolate(Sub, Point[0],  M.Size()[0], myDelimiter ); 

                cout << "Interpolate2D<" << typeid(Interpol).name() << "> at " << Point << "  yields " << T << endl; 
                cout << " IArray " << M << endl;

                return T;
#else
                return IpolDerivative<Type, Interpol, 0==DerivativeDimension>::
			interpolate(Sub, Point[0],  M.Size()[0], myDelimiter );
#endif
        }

        index_t size() const
        {
                return M.Size()[1];
        }

        Type eval() const
        {
//              printf("\nInterpolate<DerivativeDimension=%d> Eval\n", DerivativeDimension);

                return IpolDerivative<Type, Interpol, 1==DerivativeDimension>::
			interpolate(*this, Point[1],  M.Size()[1], myDelimiter );
        }
};

template <int N, class Type, class Interpol, class CoordinateType, 
          class Delimiter, int DerivativeDimension>
struct Interpolate
{
        enum { Dims = N };

        typedef Type value_type;
        typedef MultiArray<N, Type> MultiArray_t;
        typedef typename MultiArray_t::Hyperslab_t SliceStorage_t;
        typedef FixedArray<CoordinateType, N> Point_t;
        typedef FixedArray<CoordinateType, N-1> SubPoint_t;

private:
        const MultiArray<N, Type>&M;
        const Point_t&Point;

public: Delimiter myDelimiter;

        const MultiArray<N, Type>&MArray()
        {
                return M;
        }


        Interpolate(const MultiArray_t&m,
                    const Point_t&point)
        : M(m), Point(point)
        {}

        Interpolate(const MultiArray_t&m,
                    const Eagle::VVector<N , CoordinateType>&point)
        : M(m), Point(point.vec() )
        {}

        Interpolate(const MultiArray_t&m,
                    const FixedArray<CoordinateType,N>&point,
                    const Delimiter&myDeli)
        : M(m), Point(point), myDelimiter(myDeli) 
        {}

        Interpolate(const MultiArray_t&m,
                    const Eagle::VVector<N , CoordinateType>&point,
                    const Delimiter&myDeli)
        : M(m), Point(point.vec() ), myDelimiter(myDeli) 
        {}


        Type operator[](index_t i) const
        {
        MultiArray <N-1, Type> Slice = M[i]; 

        const SubPoint_t&subP = Point.get_subdim();

        Interpolate<N-1, Type, Interpol, CoordinateType, Delimiter, DerivativeDimension> 

                LI( Slice, subP, myDelimiter ); 

//              printf("\nInterpolate<%dD (rec), DerivativeDimension=%d>: Slice %d\n", N, DerivativeDimension, i);

                return IpolDerivative<Type, Interpol, N-2==DerivativeDimension>::
				interpolate(LI, Point[N-2],  M.Size()[N-2], this->myDelimiter );
        }

        index_t size() const
        {
                return M.Extension()[N-1];
        }

        Type eval() const
        {
//              printf("\nInterpolateN<DerivativeDimension=%d> Eval at %lg\n", DerivativeDimension, double(Point[N-1]) );
                return IpolDerivative<Type, Interpol, N-1==DerivativeDimension>::
			interpolate(*this, Point[N-1],  M.Size()[N-1], myDelimiter );
        }

        Type operator()() const
        {
                return eval();
        }
};


} /* namespace Fiber */ 

#endif /* __FIBER_ARRAY_INTERPOLATE_HPP */