DeriveArray.hpp

//
// $Id: DeriveArray.hpp,v 1.2 2008/01/14 19:05:50 werner Exp $
//
#ifndef __FIBER_DERIVEARRAY_HPP
#define __FIBER_DERIVEARRAY_HPP "Created 30.07.2006 22:30:21 by werner"

#include "MultiArray.hpp"
#include "IpolDelimiter.hpp"
#include "Interpolate.hpp"

namespace Fiber
{

template <class Interpol, int DerivativeDimension, class Delimiter = NoDelimiter<typename Interpol::value_type> >
struct DeriveArray
{
        template <int N, class Type> 
static  void compute(MultiArray<N, Type>&result, const MultiArray<N, Type>&D)
        {
        MultiIndex<N>   idx;
        FixedArray<double, N> point; 
        Interpolate<N, Type, Interpol , double, Delimiter, DerivativeDimension> 
                ipol(D, point); 

                do
                {
                        for(int i=0; i<N; i++)
                                point[i] = idx[i];

                        result[ idx ] = ipol.eval();
                }
                while( idx.inc( D.Size() ) );
        }

        template <int N, class Type> 
static  bool computeNthDerivative(MultiArray<N, Type>&result, const MultiArray<N, Type>&D, int NthDerivative)
        {
                if (NthDerivative == DerivativeDimension)
                {
                        compute( result, D); 
                        return true;
                } 
                return DeriveArray<Interpol, DerivativeDimension-1, Delimiter>::computeNthDerivative(result, D, NthDerivative); 
        }
};

template <class Interpol, class Delimiter>
struct DeriveArray<Interpol, -1, Delimiter>
{
        template <int N, class Type> 
static  void compute(MultiArray<N, Type>&result, const MultiArray<N, Type>&D)
        {}

        template <int N, class Type> 
static  bool computeNthDerivative(MultiArray<N, Type>&result, const MultiArray<N, Type>&D, int)
        {
                return false;
        }
};

} // namespace Fiber

#endif /* __Fiber_DERIVEARRAY_HPP */