function [F] = SCAadvect_CDS_mex(Nx, Ny, dx, dy, rho, u, v, NodeP, Nodeu, Nodev, idsP,Nbcs,NbcsDP,NbcsDu,NbcsDv)
% function [F] = SCAadvect_CDS_mex(Nx, Ny, dx, dy, P, u, v, NodeP, Nodeu,Nodev,idsP)
% 
% This function calculates the FV fluxes in a staggered grid for tracer
% advection using CDS advection scheme.
%
% 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
%       P:      Vector containing the values of Tracer, including boundary data
%       u:      Vector containing the u values of velocity, including boundary
%                 data
%       v:      Vector containing the v values of velocity, including boudary
%                 data
%       NodeP:  Node numbering matrix for the scalar
%       Nodeu:  Node numbering matrix for the u velocity
%       Nodev:  Node numbering matrix for the v velocity
%       idsP:   Location of active Tracer nodes in interior
%       Nbcs:   Number of boundary conditions
%       NbcsDP: Number of Dirichlet Tracer boundary conditions
%       NbcsDu: Number of Dirichlet u-velocity boundary conditions
%       NbcsDv: Number of Dirichlet v-velocity boundary conditions
% OUTPUTS:
%       F:      Output vector: (Fe-Fw)*dy/dV+(Fs-Fn)*dx/dV
%
% Authors:   Matt Ueckermann, Pierre Lermusiaux, Pat Haley for MIT course 2.29

%% Call C++ function to do advection. Note the casting to int32's (important)
F = SCAadvect_CDS_mx(int32(Nx),int32(Ny),dx,dy,rho,u,v,...
    int32(NodeP),int32(Nodeu),int32(Nodev),...
    int32(Nbcs),int32(NbcsDP),int32(NbcsDu),int32(NbcsDv));

F = F(idsP);
%Mexed

% %% Initialize variables
% F  = zeros(Nx*Ny,1);
% %Using upwind fluxes at boundaries, and quadratic upwind fluxes 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)
% %Fe (East flux)
% w=1;e=0;s=0;n=0; %West, East, South, North (abbreviations)
% %Select values on Right and Left (or East and West) of CV boundary (some repetition here)
% PR = P(NodeP(2-n:end-1-n,2+e:end-1+e));
% PL = P(NodeP(2+s:end-1+s,2-w:end-1-w));
% %Select quadratic left and right nodes in interior
% PR(:,2:end-1) =...
%     3/8*P(NodeP(2+s:end-1+s,3-w:end-2-w))+...
%     6/8*P(NodeP(2-n:end-1-n,3+e:end-2+e))+...
%    -1/8*P(NodeP(2-n:end-1-n,4+e:end-1+e));
% PL(:,2:end-1) =...
%     3/8*P(NodeP(2-n:end-1-n,3+e:end-2+e))+...
%     6/8*P(NodeP(2+s:end-1+s,3-w:end-2-w))+...
%    -1/8*P(NodeP(2+s:end-1+s,2-w:end-3-w));
% %Use the upwind value of P
% F(:) = reshape(F,Ny,Nx) + 0.5*dy/(dx*dy)*(...
%          u(Nodeu(2-n:end-1+s,2-w:end-1+e))  .* ( PR + PL )...
%     -abs(u(Nodeu(2-n:end-1+s,2-w:end-1+e))) .* ( PR - PL )...
%     );
% 
% %Fw (West Flux)
% %Only change in code from above is e=1;
% %Also PR is now top (north) and PL bottom (south) of CV boundary
% w=0;e=1;s=0;n=0;
% PR = P(NodeP(2-n:end-1-n,2+e:end-1+e));
% PL = P(NodeP(2+s:end-1+s,2-w:end-1-w));
% PR(:,2:end-1) =...
%     3/8*P(NodeP(2+s:end-1+s,3-w:end-2-w))+...
%     6/8*P(NodeP(2-n:end-1-n,3+e:end-2+e))+...
%    -1/8*P(NodeP(2-n:end-1-n,4+e:end-1+e));
% PL(:,2:end-1) =...
%     3/8*P(NodeP(2-n:end-1-n,3+e:end-2+e))+...
%     6/8*P(NodeP(2+s:end-1+s,3-w:end-2-w))+...
%    -1/8*P(NodeP(2+s:end-1+s,2-w:end-3-w));
% F(:) = reshape(F,Ny,Nx) - 0.5*dy/(dx*dy)*(...
%          u(Nodeu(2-n:end-1+s,2-w:end-1+e))  .* ( PR + PL )...
%     -abs(u(Nodeu(2-n:end-1+s,2-w:end-1+e))) .* ( PR - PL )...
%     );
% 
% %Fs (South Flux)
% %Change in code here is using v instead of u, and s=1.
% w=0;e=0;s=1;n=0;
% PR = P(NodeP(2-n:end-1-n,2+e:end-1+e));
% PL = P(NodeP(2+s:end-1+s,2-w:end-1-w));
% PR(2:end-1,:) =...
%     3/8*P(NodeP(3+s:end-2+s,2-w:end-1-w))+...
%     6/8*P(NodeP(3-n:end-2-n,2+e:end-1+e))+...
%    -1/8*P(NodeP(2-n:end-3-n,2+e:end-1+e));
% PL(2:end-1,:) =...
%     3/8*P(NodeP(3-n:end-2-n,2+e:end-1+e))+...
%     6/8*P(NodeP(3+s:end-2+s,2-w:end-1-w))+...
%    -1/8*P(NodeP(4+s:end-1+s,2-w:end-1-w));
% F(:) = reshape(F,Ny,Nx) + 0.5*dx/(dx*dy)*(...
%          v(Nodev(2-n:end-1+s,2-e:end-1+w))  .* ( PR + PL )...
%     -abs(v(Nodev(2-n:end-1+s,2-e:end-1+w))) .* ( PR - PL )...
%     );
% 
% %Fn (North Flux)
% %change in code from above is n=1.
% w=0;e=0;s=0;n=1;
% PR = P(NodeP(2-n:end-1-n,2+e:end-1+e));
% PL = P(NodeP(2+s:end-1+s,2-w:end-1-w));
% PR(2:end-1,:) =...
%     3/8*P(NodeP(3+s:end-2+s,2-w:end-1-w))+...
%     6/8*P(NodeP(3-n:end-2-n,2+e:end-1+e))+...
%    -1/8*P(NodeP(2-n:end-3-n,2+e:end-1+e));
% PL(2:end-1,:) =...
%     3/8*P(NodeP(3-n:end-2-n,2+e:end-1+e))+...
%     6/8*P(NodeP(3+s:end-2+s,2-w:end-1-w))+...
%    -1/8*P(NodeP(4+s:end-1+s,2-w:end-1-w));
% F(:) = reshape(F,Ny,Nx) - 0.5*dx/(dx*dy)*(...
%          v(Nodev(2-n:end-1+s,2-e:end-1+w))  .* ( PR + PL )...
%     -abs(v(Nodev(2-n:end-1+s,2-e:end-1+w))) .* ( PR - PL )...
%     );
% F = F(idsP);
