function [Fu Fv] = UVadvect(Nx, Ny, dx, dy, u, v, Nodeu, Nodev, idsu, idsv)
% function [Fu Fv] = UVadvect(Nx, Ny, dx, dy, u, v, Nodeu, Nodev, idsu, idsv)
% This function calculates the FV fluxes in a staggered grid for scalar
% 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

%% Prelinimaries
%% Initialize variables
Fu = zeros((Nx-1)* Ny   ,1);
Fv = zeros( Nx   *(Ny-1),1);
%Using upwind fluxes
%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
%Fe
e=1;w=0;s=0;n=0;
Fu(:) = reshape(Fu,Ny,Nx-1) + dy/(dx*dy)*0.5*(...
    u(Nodeu(2-n:end-1-n,2+w:end-1+w)).^2 - abs(u(Nodeu(2-n:end-1-n,2+w:end-1+w))).*u(Nodeu(2-n:end-1-n,2+w:end-1+w))...
  + u(Nodeu(2+s:end-1+s,2-e:end-1-e)).^2 + abs(u(Nodeu(2+s:end-1+s,2-e:end-1-e))).*u(Nodeu(2+s:end-1+s,2-e:end-1-e))...
  );
%Fw
e=0;w=1;s=0;n=0;
Fu(:) = reshape(Fu,Ny,Nx-1) - dy/(dx*dy)*0.5*(...
    u(Nodeu(2-n:end-1-n,2+w:end-1+w)).^2 - abs(u(Nodeu(2-n:end-1-n,2+w:end-1+w))).*u(Nodeu(2-n:end-1-n,2+w:end-1+w))...
  + u(Nodeu(2+s:end-1+s,2-e:end-1-e)).^2 + abs(u(Nodeu(2+s:end-1+s,2-e:end-1-e))).*u(Nodeu(2+s:end-1+s,2-e:end-1-e))...
  );
%Fs
e=0;w=0;s=1;n=0;
V=(v(Nodev(2-n:end-1+s,2:end-2))+v(Nodev(2-n:end-1+s,3:end-1)))/2;
Fu(:) = reshape(Fu,Ny,Nx-1) + dx/(dx*dy)*0.5*(...
    u(Nodeu(2-n:end-1-n,2+w:end-1+w)).*V - abs(V).*u(Nodeu(2-n:end-1-n,2+w:end-1+w))...
  + u(Nodeu(2+s:end-1+s,2-e:end-1-e)).*V + abs(V).*u(Nodeu(2+s:end-1+s,2-e:end-1-e))...
  );
%Fn
e=0;w=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;
Fu(:) = reshape(Fu,Ny,Nx-1) - dx/(dx*dy)*0.5*(...
    u(Nodeu(2-n:end-1-n,2+w:end-1+w)).*V - abs(V).*u(Nodeu(2-n:end-1-n,2+w:end-1+w))...
  + u(Nodeu(2+s:end-1+s,2-e:end-1-e)).*V + abs(V).*u(Nodeu(2+s:end-1+s,2-e:end-1-e))...
  );
Fu=Fu(idsu);
%FV
%Fe
e=1;w=0;s=0;n=0;
V=(u(Nodeu(2:end-2,2-e:end-1+w))+u(Nodeu(3:end-1,2-e:end-1+w)))/2;
Fv(:) = reshape(Fv,Ny-1,Nx) + dy/(dx*dy)*0.5*(...
    v(Nodev(2-n:end-1-n,2+w:end-1+w)).*V - abs(V).*v(Nodev(2-n:end-1-n,2+w:end-1+w))...
  + v(Nodev(2+s:end-1+s,2-e:end-1-e)).*V + abs(V).*v(Nodev(2+s:end-1+s,2-e:end-1-e))...
  );
%Fw
e=0;w=1;s=0;n=0;
V=(u(Nodeu(2:end-2,2-e:end-1+w))+u(Nodeu(3:end-1,2-e:end-1+w)))/2;
Fv(:) = reshape(Fv,Ny-1,Nx) - dy/(dx*dy)*0.5*(...
    v(Nodev(2-n:end-1-n,2+w:end-1+w)).*V - abs(V).*v(Nodev(2-n:end-1-n,2+w:end-1+w))...
  + v(Nodev(2+s:end-1+s,2-e:end-1-e)).*V + abs(V).*v(Nodev(2+s:end-1+s,2-e:end-1-e))...
  );
%Fs
e=0;w=0;s=1;n=0;
Fv(:) = reshape(Fv,Ny-1,Nx) + dx/(dx*dy)*0.5*(...
    v(Nodev(2-n:end-1-n,2+w:end-1+w)).^2 - abs(v(Nodev(2-n:end-1-n,2+w:end-1+w))).*v(Nodev(2-n:end-1-n,2+w:end-1+w))...
  + v(Nodev(2+s:end-1+s,2-e:end-1-e)).^2 + abs(v(Nodev(2+s:end-1+s,2-e:end-1-e))).*v(Nodev(2+s:end-1+s,2-e:end-1-e))...
  );
%Fn
e=0;w=0;s=0;n=1;
Fv(:) = reshape(Fv,Ny-1,Nx) - dx/(dx*dy)*0.5*(...
    v(Nodev(2-n:end-1-n,2+w:end-1+w)).^2 - abs(v(Nodev(2-n:end-1-n,2+w:end-1+w))).*v(Nodev(2-n:end-1-n,2+w:end-1+w))...
  + v(Nodev(2+s:end-1+s,2-e:end-1-e)).^2 + abs(v(Nodev(2+s:end-1+s,2-e:end-1-e))).*v(Nodev(2+s:end-1+s,2-e:end-1-e))...
  );
Fv=Fv(idsv);
