function [Fu Fv] = UVadvect_CDS(Nx, Ny, dx, dy, u, v, Nodeu, Nodev, idsu, idsv)
% function [Fu Fv] = UVadvect_CDS(Nx, Ny, dx, dy, u, v, Nodeu, Nodev, idsu, idsv)
% This function calculates the FV fluxes in a staggered grid for velocity
% 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
%       u:      Vector containing the u values of velocity, including boundary
%                 data
%       v:      Vector containing the v values of velocity, including boudary
%                 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
% 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, Pierre Lermusiaux, Pat Haley for MIT course 2.29


%% 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 from n(orth) to s(south) index 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)
UR = u(Nodeu(2-n:end-1-n,2+e:end-1+e));
UL = u(Nodeu(2+s:end-1+s,2-w:end-1-w));

%Select central value in flux equation
Fu(:) = reshape(Fu,Ny,Nx-1) + dy/(dx*dy)*0.25*(...
    (UR+UL).^2 ...
  );

%Fw
%Only change from above is e=1;
w=0;e=1;s=0;n=0;
UR = u(Nodeu(2-n:end-1-n,2+e:end-1+e));
UL = u(Nodeu(2+s:end-1+s,2-w:end-1-w));

%Select central value in flux equation
Fu(:) = reshape(Fu,Ny,Nx-1) - dy/(dx*dy)*0.25*(...
    (UR+UL).^2 ...
  );

%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=(v(Nodev(2-n:end-1+s,2:end-2))+v(Nodev(2-n:end-1+s,3:end-1)))/2;
%Find quadratic Right and left velocities as before.
UR = u(Nodeu(2-n:end-1-n,2+e:end-1+e));
UL = u(Nodeu(2+s:end-1+s,2-w:end-1-w));

Fu(:) = reshape(Fu,Ny,Nx-1) + dx/(dx*dy)*0.5*(...
    (UR+UL).*V ...
  );

%Fn Change from above is n=1
w=0;e=0;s=0;n=1;
V=(v(Nodev(2-n:end-1+s,2:end-2))+v(Nodev(2-n:end-1+s,3:end-1)))/2;
UR = u(Nodeu(2-n:end-1-n,2+e:end-1+e));
UL = u(Nodeu(2+s:end-1+s,2-w:end-1-w));
Fu(:) = reshape(Fu,Ny,Nx-1) - dx/(dx*dy)*0.5*(...
    (UR+UL).*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=(u(Nodeu(2:end-2,2-w:end-1+e))+u(Nodeu(3:end-1,2-w:end-1+e)))/2;
UR = v(Nodev(2-n:end-1-n,2+e:end-1+e));
UL = v(Nodev(2+s:end-1+s,2-w:end-1-w));

Fv(:) = reshape(Fv,Ny-1,Nx) + dy/(dx*dy)*0.5*(...
    (UR+UL).*V ...
  );
%Fw
w=0;e=1;s=0;n=0;
V=(u(Nodeu(2:end-2,2-w:end-1+e))+u(Nodeu(3:end-1,2-w:end-1+e)))/2;
UR = v(Nodev(2-n:end-1-n,2+e:end-1+e));
UL = v(Nodev(2+s:end-1+s,2-w:end-1-w));

Fv(:) = reshape(Fv,Ny-1,Nx) - dy/(dx*dy)*0.5*(...
    (UR+UL).*V ...
  );
%Fs
w=0;e=0;s=1;n=0;
UR = v(Nodev(2-n:end-1-n,2+e:end-1+e));
UL = v(Nodev(2+s:end-1+s,2-w:end-1-w));

Fv(:) = reshape(Fv,Ny-1,Nx) + dx/(dx*dy)*0.25*(...
    (UR+UL).^2 ...
  );
%Fn
w=0;e=0;s=0;n=1;
UR = v(Nodev(2-n:end-1-n,2+e:end-1+e));
UL = v(Nodev(2+s:end-1+s,2-w:end-1-w));

Fv(:) = reshape(Fv,Ny-1,Nx) - dx/(dx*dy)*0.25*(...
    (UR+UL).^2 ...
  );

%reformat flux for output
Fv=Fv(idsv);
