//#####################################################################
// 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_HELPERS__
#define __INTERPOLATION_HELPERS__

#include <Core/Arrays/ARRAY.h>
#include <Core/Grids/GRID.h>
#include <Core/Log/DEBUG_UTILITIES.h>
namespace PhysBAM
{
namespace INTERPOLATION_HELPERS
{
// stencil_size: 2=linear, 3=quadratic, 4=cubic, ...
// returned index is lowest indexed value to be used
template<int stencil_size,bool const_extrap,bool is_mac,class IV,class TV>
static std::pair<IV,TV> Clamped_Index(const GRID<TV>& grid,GRID_TYPE<is_mac>,
    const RANGE<IV>& domain,const TV& X)
{
    typedef typename TV::SCALAR T;
    TV Y=(X-grid.domain.min_corner)*grid.one_over_dX-(T).5*(stencil_size-2+is_mac);
    IV Z=clamp(IV(floor(Y)),domain.min_corner,domain.max_corner-stencil_size);
    TV U=Y-TV(Z);
    if constexpr(const_extrap) U=clamp(U,TV(),TV()+1);
    return {Z,U};
}

// stencil_size: 2=linear, 3=quadratic, 4=cubic, ...
// returned index is lowest indexed value to be used
template<int stencil_size,bool const_extrap,class IV,class TV>
static std::pair<IV,TV> Clamped_Index_Face(const GRID<TV>& grid,
    const RANGE<IV>& domain,const TV& X,int axis)
{
    typedef typename TV::SCALAR T;
    TV s=TV()+(T).5*(stencil_size-1);
    TV Y=(X-grid.domain.min_corner)*grid.one_over_dX-s.Add_Axis(axis,-(T).5);
    IV Z=clamp(IV(floor(Y)),domain.min_corner,domain.max_corner-stencil_size);
    TV U=Y-TV(Z);
    if constexpr(const_extrap) U=clamp(U,TV(),TV()+1);
    return {Z,U};
}

// stencil_size: 2=linear, 3=quadratic, 4=cubic, ...
// returned index is lowest indexed value to be used
template<int stencil_size,bool const_extrap,class IV,class TV>
static VECTOR<std::pair<IV,TV>,TV::m> Clamped_Index_Face(const GRID<TV>& grid,
    const RANGE<IV>& domain,const TV& X)
{
    typedef typename TV::SCALAR T;
    VECTOR<std::pair<IV,TV>,TV::m> p;
    TV s=TV()+(T).5*(stencil_size-1);
    TV W=(X-grid.domain.min_corner)*grid.one_over_dX-s;
    for(int i=0;i<TV::m;i++)
    {
        TV Y=W.Add_Axis(i,(T).5);
        IV Z=clamp(IV(floor(Y)),domain.min_corner,domain.max_corner.Add_Axis(i,1)-stencil_size);
        TV U=Y-TV(Z);
        if constexpr(const_extrap) U=clamp(U,TV(),TV()+1);
        p(i).first=Z;
        p(i).second=U;
    }
    return p;
}
}
}
#endif
