//##################################################################### 
// Copyright 2021, Craig Schroeder.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
#ifndef __INTERPOLATION_SPLINES__
#define __INTERPOLATION_SPLINES__

#include <Core/Vectors/VECTOR.h>
#include <array>

namespace PhysBAM {
template<int k,class T>
VECTOR<T,k+1> Make_Powers(T a)
{
    VECTOR<T,k+1> u;
    T x=1;
    for(int i=0;i<k+1;i++)
    {
        u(i)=x;
        x*=a;
    }
    return u;
}

template<class T,size_t r,size_t d,int k>
VECTOR<T,r> Evaluate_Spline(const std::array<std::array<double,d>,r>& A,const VECTOR<T,k>& u)
{
    static_assert(d<=k,"not enough powers available to evaluate spline.");
    VECTOR<T,r> v;
    for(int j=0;j<d;j++)
    {
        T x=u[j];
        for(int i=0;i<r;i++)
            v(i)+=A[i][j]*x;
    }
    return v;
}

template<class T,size_t r,size_t d,int k>
VECTOR<T,r> Evaluate_Spline_Diff(const std::array<std::array<double,d>,r>& A,const VECTOR<T,k>& u)
{
    static_assert(d-1<=k,"not enough powers available to evaluate spline.");
    VECTOR<T,r> v;
    for(int j=1;j<d;j++)
    {
        T x=j*u[j-1];
        for(int i=0;i<r;i++)
            v(i)+=A[i][j]*x;
    }
    return v;
}

template<class T,size_t r,size_t d,int k>
VECTOR<T,r> Evaluate_Spline_Hess(const std::array<std::array<double,d>,r>& A,const VECTOR<T,k>& u)
{
    static_assert(d-2<=k,"not enough powers available to evaluate spline.");
    VECTOR<T,r> v;
    for(int j=2;j<d;j++)
    {
        T x=j*(j-1)*u[j-2];
        for(int i=0;i<r;i++)
            v(i)+=A[i][j]*x;
    }
    return v;
}

namespace IS {
extern const std::array<std::array<double,2>,2> B1;
extern const std::array<std::array<double,3>,3> B2;
extern const std::array<std::array<double,4>,4> B3;

extern const std::array<std::array<double,3>,2> C2;
extern const std::array<std::array<double,4>,3> C3;
extern const std::array<std::array<double,5>,4> C4;

extern const std::array<std::array<double,4>,2> D3;
extern const std::array<std::array<double,5>,3> D4;
extern const std::array<std::array<double,6>,4> D5;

extern const std::array<std::array<double,5>,3> F4;
extern const std::array<std::array<double,6>,4> F5;

extern const std::array<std::array<double,5>,3> G4;
extern const std::array<std::array<double,6>,4> G5;

extern const std::array<std::array<double,4>,3> H3;
extern const std::array<std::array<double,5>,4> H4;

extern const std::array<std::array<double,4>,4> I3;
extern const std::array<std::array<double,4>,6> I3p;
extern const std::array<std::array<double,4>,8> I3pp;
extern const std::array<std::array<double,4>,4> J3;
extern const std::array<std::array<double,4>,6> J3p;
extern const std::array<std::array<double,4>,8> J3pp;
extern const std::array<std::array<double,2>,6> B1pp;
extern const std::array<std::array<double,2>,8> B1ppp;
}

namespace IS4 {
extern const std::array<std::array<double,3>,5> B2;
extern const std::array<std::array<double,3>,7> B2p;
extern const std::array<std::array<double,5>,7> I4;
extern const std::array<std::array<double,5>,7> J4;
}

}
#endif
