function [Fu Fv] = UVadvect_DO_CDS_mex(Nx, Ny, dx, dy, ui, vi, uj, vj, Nodeu, Nodev, idsu, idsv,Nbcs,NbcsDu,NbcsDv)
% function [Fu Fv] = UVadvect_DO_CDS_mex(Nx, Ny, dx, dy, u, v, Nodeu, Nodev, idsu, idsv)
% This function calculates the FV fluxes in a staggered grid for velocity
% advection, and allows for different velocities to be used to determine the
% upwind direction (as required by the DO implementation). Central Difference
% scheme used for advection.
%
% INPUTS:
%       Nx:     Number of interior nodes in the x-direction
%       Ny:     Number of interior nodes in the y-direction
%       dx:     Step size in x direction
%       dy:     Step size in y direction
%       ui:     Vector containing the i'th u mode of velocity, including boundary
%                 data
%       vi:     Vector containing the i'th v mode of velocity, including boundary
%                 data
%       uj:     Vector containing the j'th u mode of velocity, including boundary
%                 data
%       vj:     Vector containing the j'th v mode of velocity, including boundary
%                 data
%       Nodeu:  Is the node numbering matrix for the u velocity
%       Nodev:  Is the node numbering matrix for the v velocity
%       idsu:   Location of active u-velocity nodes in interior
%       idsv:   Location of active v-velocity nodes in interior
%       Nbcs:   Number of boundary conditions
%       NbcsDu: Number of Dirichlet u-velocity boundary conditions
%       NbcsDv: Number of Dirichlet v-velocity boundary conditions
%
% OUTPUTS:
%       Fu:      Output vector:Fu = (Fue-Fuw)*dy/dV+(Fus-Fun)*dx/dV
%       Fv:      Output vector:Fv = (Fve-Fvw)*dy/dV+(Fvs-Fvn)*dx/dV
%
% Authors:   Matt Ueckermann, Themis Sapsis, 
%            Pierre Lermusiaux, Pat Haley for MIT course 2.29
%
% NOTE: In this version of the advection equation, there is no upwind direction.
% This is consistent for the modes which should have no sign attached to them.
%
% Mexed

%% Call C++ function to do advection. Note the casting to int32's (important)
[Fu Fv] = UVadvect_DO_CDS_mx(int32(Nx),int32(Ny),dx,dy,ui,vi,uj,vj,...
    int32(Nodeu),int32(Nodev),int32(Nbcs),int32(NbcsDu),int32(NbcsDv));
Fu=Fu(idsu);
Fv=Fv(idsv);

%Mexed stuff below
%% Initialize variables
% Fu = zeros((Nx-1)* Ny   ,1);
% Fv = zeros( Nx   *(Ny-1),1);
%Using upwind fluxes at boundaries and QUICK scheme in interior
%NOTE: There is an annoying thing where n-s I decreases, whereas in x-y
%coordinates, I would increase.

% 
% 
% %% Calculate fluxes (Vectorized for speed)
% % FU (Calculate fluxes for u-velocity)
% %Fe 
% w=1;e=0;s=0;n=0;
% %Select values on Right and Left of CV boundary (some repetition here)
% UiR = ui(Nodeu(2-n:end-1-n,2+e:end-1+e));
% UiL = ui(Nodeu(2+s:end-1+s,2-w:end-1-w));
% %Now do the j'th mode
% UjR = uj(Nodeu(2-n:end-1-n,2+e:end-1+e));
% UjL = uj(Nodeu(2+s:end-1+s,2-w:end-1-w));
% %Select upwind value (based on j-mode velocity) in flux equation 
% Fu(:) = reshape(Fu,Ny,Nx-1) + dy/(dx*dy)*0.25*(...
%     (UiR+UiL).*(UjR + UjL)...
%   );
% 
% %Fw
% %Only change from above is e=1;
% w=0;e=1;s=0;n=0;
% %Select values on Right and Left of CV boundary (some repetition here)
% UiR = ui(Nodeu(2-n:end-1-n,2+e:end-1+e));
% UiL = ui(Nodeu(2+s:end-1+s,2-w:end-1-w));
% %Now do the j'th mode
% UjR = uj(Nodeu(2-n:end-1-n,2+e:end-1+e));
% UjL = uj(Nodeu(2+s:end-1+s,2-w:end-1-w));
% 
% %Select upwind value (based on j-mode velocity) in flux equation 
% Fu(:) = reshape(Fu,Ny,Nx-1) - dy/(dx*dy)*0.25*(...
%     (UiR+UiL).*(UjR + UjL)...
%   );
% 
% %Fs
% %changes here is using V velocity located at CV wall, and s=1
% w=0;e=0;s=1;n=0;
% %Average the two velocities located on each side of the CV
% V=(vj(Nodev(2-n:end-1+s,2:end-2))+vj(Nodev(2-n:end-1+s,3:end-1)))/2;
% %Find quadratic Right and left velocities as before.
% UiR = ui(Nodeu(2-n:end-1-n,2+e:end-1+e));
% UiL = ui(Nodeu(2+s:end-1+s,2-w:end-1-w));
% 
% Fu(:) = reshape(Fu,Ny,Nx-1) + dx/(dx*dy)*0.5*(...
%     (UiR+UiL).*V...
%   );
% %Fn Change from above is n=1
% w=0;e=0;s=0;n=1;
% %Average the two velocities located on each side of the CV
% V=(vj(Nodev(2-n:end-1+s,2:end-2))+vj(Nodev(2-n:end-1+s,3:end-1)))/2;
% UiR = ui(Nodeu(2-n:end-1-n,2+e:end-1+e));
% UiL = ui(Nodeu(2+s:end-1+s,2-w:end-1-w));
% 
% Fu(:) = reshape(Fu,Ny,Nx-1) - dx/(dx*dy)*0.5*(...
%     (UiR+UiL).*V...
%   );
% %Reformat flux for output
% Fu=Fu(idsu);


% %% FV
% %Very much a repetition of FU, but with index anamolies due to grid
% %stagerring :)
% %Fe
% w=1;e=0;s=0;n=0;
% V=(uj(Nodeu(2:end-2,2-w:end-1+e))+uj(Nodeu(3:end-1,2-w:end-1+e)))/2;
% UiR = vi(Nodev(2-n:end-1-n,2+e:end-1+e));
% UiL = vi(Nodev(2+s:end-1+s,2-w:end-1-w));
% 
% Fv(:) = reshape(Fv,Ny-1,Nx) + dy/(dx*dy)*0.5*(...
%     (UiR+UiL).*V...
%   );
% 
% %Fw
% w=0;e=1;s=0;n=0;
% V=(uj(Nodeu(2:end-2,2-w:end-1+e))+uj(Nodeu(3:end-1,2-w:end-1+e)))/2;
% UiR = vi(Nodev(2-n:end-1-n,2+e:end-1+e));
% UiL = vi(Nodev(2+s:end-1+s,2-w:end-1-w));
% 
% Fv(:) = reshape(Fv,Ny-1,Nx) - dy/(dx*dy)*0.5*(...
%     (UiR+UiL).*V...
%   );
% 
% %Fs
% w=0;e=0;s=1;n=0;
% %i'th mode
% UiR = vi(Nodev(2-n:end-1-n,2+e:end-1+e));
% UiL = vi(Nodev(2+s:end-1+s,2-w:end-1-w));
% 
% %j'the mode
% UjR = vj(Nodev(2-n:end-1-n,2+e:end-1+e));
% UjL = vj(Nodev(2+s:end-1+s,2-w:end-1-w));
% 
% %advect
% Fv(:) = reshape(Fv,Ny-1,Nx) + dx/(dx*dy)*0.25*(...
%      (UiR+UiL).*(UjR + UjL)...
%   );
% 
% %Fn
% w=0;e=0;s=0;n=1;
% %i'th mode
% UiR = vi(Nodev(2-n:end-1-n,2+e:end-1+e));
% UiL = vi(Nodev(2+s:end-1+s,2-w:end-1-w));
% 
% %j'the mode
% UjR = vj(Nodev(2-n:end-1-n,2+e:end-1+e));
% UjL = vj(Nodev(2+s:end-1+s,2-w:end-1-w));
% 
% %advect
% Fv(:) = reshape(Fv,Ny-1,Nx) - dx/(dx*dy)*0.25*(...
%     (UiR+UiL).*(UjR + UjL)... 
%   );
% 
% %reformat flux for output
% Fv=Fv(idsv);
