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
{
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.Size()[0] /*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.Size()[0] /*Extension()[0]*/, myDelimiter );
#endif
        }
 
        index_t size() const
        {
                return M.Size().maxidx();
        }
 
 
        Type eval() const
        {
        //      printf("\nInterpolate<DerivativeDimension=%d> Eval\n", DerivativeDimension);
 
                return IpolDerivative<Type, Interpol, 0==DerivativeDimension>::
			interpolate(M, Point[0],  M.Size()[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.Size()[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 */