function [EdgeVelocities] = AdvectionVelocity(AdvectionField, NatNeighborPairs, PosX, PosY, I, J, NumInsideNodes, AdvectionDiscrModel, varargin)

% Evaluates the advection field function at the edges of the control
% volumes using the input variables. AdvectionDiscrModel could be
% 'CentralDifference' method only.
%
% Two important notes: 
%   1- The independent variables should be defined inside the domain only.
%   2- The AdvectionField should be passed as a 'cell' with two elements
%   (2D), each a function or a constant.
%   2- The functions should be defined using Element-By-Element notations.
%
% For example, if advection velocity is constant [1 0.5], then
% AdvectionField should be specified as:  AdvectionField = {1 0.5}, and
% called as:
% [EdgeVelocities] = AdvectionVelocity(AdvectionField, NatNeighborPairs, ...
%      PosX, PosY, I, J, NumInsideNodes, 'CentralDifference')
%
%
% Or, if we have a non-constant advection field that depends on coordinates
% [X -Y], then AdvectionField should be specified as a cell function:  
% AdvectionField = {@(x,y) x; @(x,y) -y}, and called as:
%
% X_local = PosX(InsideOldMask);  Y_local = PosY(InsideOldMask);
% [EdgeVelocities] = AdvectionVelocity(AdvectionField, NatNeighborPairs, ...
%      PosX, PosY, I, J, NumInsideNodes, 'CentralDifference',X_local,Y_local)


% Create the edge unit normal vectors
II = NatNeighborPairs(:,1);
JJ = NatNeighborPairs(:,2);
NormalsX = PosX(JJ)-PosX(II);
NormalsY = PosY(JJ)-PosY(II);
Lengths = sqrt(NormalsX.*NormalsX + NormalsY.*NormalsY);
Normals = [NormalsX./Lengths NormalsY./Lengths];

% evaluate velocity functions at the cell centers
AdvectionFieldCell = [AdvectionField{1}(varargin{:}).*ones(NumInsideNodes,1) AdvectionField{2}(varargin{:}).*ones(NumInsideNodes,1)];

% use normal vectors and cell velocities to approximate velocities at the
% edges of the CVs

switch AdvectionDiscrModel
  case 'CentralDifference'
    EdgeVelocities = sum((AdvectionFieldCell(I,:)+AdvectionFieldCell(J,:))/2.*Normals,2);
  otherwise
    
end

end % end function 'AdvectionVelocity'