//#####################################################################
// Class ADVECTION_SECOND_ORDER_SEMI_LAGRANGIAN
//#####################################################################
#ifndef __ADVECTION_SECOND_ORDER_SEMI_LAGRANGIAN__
#define __ADVECTION_SECOND_ORDER_SEMI_LAGRANGIAN__

#include <Core/Advection/ADVECTION.h>
#include <Core/Grids/NODE_ITERATOR.h>

namespace PhysBAM {

template<class TV, class T2, class T_Z_INTERPOLATION, class T_V_INTERPOLATION, class T_FACE_ARRAY=ARRAY<T2,FACE_INDEX<TV::m> >>
class ADVECTION_SECOND_ORDER_SEMI_LAGRANGIAN : public ADVECTION<TV,T2,T_FACE_ARRAY> 
{
    typedef typename TV::SCALAR T;typedef VECTOR<int,TV::m> IV;

public:
    void Update_Advection_Equation_Node(const GRID<TV>& grid, ARRAY<T2,IV>& Z, const ARRAY<T2,IV>& Z_ghost, 
            const ARRAY<TV,IV>& V, const T dt, const T time,
            const ARRAY<T2,IV>* Z_min_ghost=0, const ARRAY<T2,IV>* Z_max_ghost=0, ARRAY<T2,IV>* Z_min=0, ARRAY<T2,IV>* Z_max=0,
        ARRAY<T2,FACE_INDEX<TV::m> >* flux=0)
    {
        T_Z_INTERPOLATION z_interpolation;
        T_V_INTERPOLATION v_interpolation;

        for (NODE_ITERATOR<TV> it(grid); it.Valid(); it.Next()) {
            IV node = it.Node_Index();
            TV X = it.Location();
            TV X1 = X - (T)0.5*dt*V(node);

            TV V1 = v_interpolation.Interpolate_Node(grid, V, X1); 

            TV X2 = X - dt * V1;

            Z(node) = z_interpolation.Interpolate_Node(grid, Z_ghost, X2);

            if (Z_min && Z_max) {
                VECTOR<T2,2> extrema=z_interpolation.Extrema_Node(grid,*Z_min_ghost,*Z_max_ghost,X2);
                (*Z_min)(node)=extrema.x;
                (*Z_max)(node)=extrema.y;
            }
        }

    }


    void Update_Advection_Equation_Cell_Lookup(const GRID<TV>& grid,ARRAY<T2,IV>& Z,const ARRAY<T2,IV>& Z_ghost,
        const T_FACE_ARRAY& face_velocities,const T dt,const T time,
        const ARRAY<T2,IV>* Z_min_ghost,const ARRAY<T2,IV>* Z_max_ghost,ARRAY<T2,IV>* Z_min,ARRAY<T2,IV>* Z_max,
        ARRAY<T2,FACE_INDEX<TV::m> >* flux)
    {
        T_Z_INTERPOLATION z_interpolation;
        T_V_INTERPOLATION v_interpolation;

        for(CELL_ITERATOR<TV> it(grid);it.Valid();it.Next()){
            IV cell=it.Cell_Index();

            TV X = it.Location();
            TV X1 = X - (T)0.5*dt*v_interpolation.Interpolate_Face(grid, face_velocities, X); 
            TV V1 = v_interpolation.Interpolate_Face(grid, face_velocities, X1); 

            TV X2 = X - dt * V1;

            Z(cell)=z_interpolation.Interpolate_Cell(grid,Z_ghost,X2);

            if(Z_min && Z_max) {
                VECTOR<T2,2> extrema=z_interpolation.Extrema_Cell(grid,*Z_min_ghost,*Z_max_ghost,X2);
                (*Z_min)(cell)=extrema.x;
                (*Z_max)(cell)=extrema.y;
            }
        }
    }

    void Update_Advection_Equation_Face_Lookup(const GRID<TV>& grid,ARRAY<T,FACE_INDEX<TV::m> >& Z,const T_FACE_ARRAY& Z_ghost,
            const T_FACE_ARRAY& face_velocities,const T dt,const T time,
            const T_FACE_ARRAY* Z_min_ghost,const T_FACE_ARRAY* Z_max_ghost,ARRAY<T,FACE_INDEX<TV::m> >* Z_min,ARRAY<T,FACE_INDEX<TV::m> >* Z_max,
        VECTOR<ARRAY<T2,FACE_INDEX<TV::m> >,TV::m>* flux)
    {
        T_V_INTERPOLATION v_interpolation;
        T_Z_INTERPOLATION z_interpolation;

        for(FACE_ITERATOR<TV> it(grid);it.Valid();it.Next()){
            IV face=it.Face_Index();
            int axis=it.Axis();

            TV X = it.Location();
            TV X1 = X - (T)0.5*dt*v_interpolation.Interpolate_Face(grid, face_velocities, X);
            TV V1 = v_interpolation.Interpolate_Face(grid, face_velocities, X1); // Should this be the compoent one?

            TV X2 = X - dt * V1;

            Z.Component(axis)(face)=z_interpolation.Interpolate_Face(grid,Starting_Point_Face(Z_ghost,axis,face),X2,axis);

            if(Z_min && Z_max) {
                VECTOR<T,2> extrema=z_interpolation.Extrema_Face(grid,*Starting_Point_Face(Z_max_ghost,axis,face),*Starting_Point_Face(Z_max_ghost,axis,face),X2,axis);
                (*Z_min).Component(axis)(face)=extrema.x;(*Z_max).Component(axis)(face)=extrema.y;
            }
        }
    }
};

}

#endif 
