function [Af,RHS,Arhs] = mk_DiffusionOperator(Nx, Ny, dx, dy, bcsD, bcsN, Node, stagr,...
    Rm_Singularity)
% function [Af,RHS,Arhs] = mk_DiffusionOperator(Nx, Ny, dx, dy, bcsD, bcsN, Node, stagr)
% This function creates the diffusion operator and its and RHS 
%
% INPUTS:
%       Nx:     Number of interior tracer nodes in the x-direction
%       Ny:     Number of interior tracer nodes in the y-direction
%       dx:     Step size in x direction
%       dy:     Step size in y direction
%       bcsD:   Vector containing the values of dirichlet boundary
%                   conditions. 
%       bcsN:   Vector containing the values of Neumann boundary
%                   conditions multiplied by the correct sign and length
%                   eg. bcsN = gN*dx, where du/dx = gN on the boundary
%       Node:   Is the node numbering matrix. size(Node) = [Ny+2,Nx+2]
%               NOTE: Node needs to include boundary conditions. In the 
%                       case for periodic boundaries, the periodicity needs
%                       to be included. See NodePad.m
%       stagr:  String: Which grid to use, 'P', 'u', or 'v'. 
% OUTPUTS:
%       Af:     The A matrix
%       RHS:    The Right-hand-side vector from BC contributions
%       Arhs:   The matrix used to create the Right-hand-side vector
%
% See also: NodePad.m
%
% Authors:   Matt Ueckermann, Pierre Lermusiaux, Pat Haley for MIT course 2.29

%% Preliminaries

% For staggered u and v grids, the number of interior unknowns require modification.
switch stagr
    case 'u', Nx = Nx-1;
    case 'v', Ny = Ny-1;
end;

% Initialize the A matrix. Notice we are initializing it as a sparse
% matrix in order to reduce memory requirements. For a large number of
% nodes, this sparse matrix allocation would be essential, otherwise we may
% run out of memory.
Nbcs = length(bcsD) + length(bcsN);
TotalNodes = max(Node(:)) - Nbcs;

% Initialize the RHS vector (in Ax=b, this is the 'b' vector)
RHS = zeros(TotalNodes,1);


%% Build matrix A

Rowij = Node(2 : Ny + 1, 2 : Nx + 1) - Nbcs;
idsx = 2 : (Nx + 1);
idsy = 2 : (Ny + 1);
%Create the indices for the entries
I = [Rowij, Rowij, Rowij, Rowij, Rowij];
J = [Node(idsy, idsx), Node(idsy, idsx - 1), Node(idsy, idsx + 1),...
    Node(idsy - 1, idsx), Node(idsy + 1, idsx)];
J(I <= 0) = [];
I(I <= 0) = [];

%Create the A matrix
One = ones(length(I(:)) / 5,1);
A = sparse(I(:), J(:),...
    [(2 / dx ^ 2 + 2 / dy^2) * One; 
    - 1 / dx ^ 2 * One; - 1 / dx ^ 2 * One; 
    - 1 / dy ^ 2 * One; - 1 / dy ^ 2 * One],...
    TotalNodes, TotalNodes + Nbcs, 5 * TotalNodes + Nbcs * TotalNodes);


%% Set up the Dirichlet and Neumann boundary conditions

%A quick formatting check
if ~isempty(bcsD)
    %Then we add the boundary condition contributions to the RHS vector
    RHS = RHS - A(:, 1:Nbcs-length(bcsN))*bcsD;
end;
% Neumann boundary conditions require a little more work. For Neumann on
% the west face, we have u_w = u_p + dx*bcN. We subtitute this in the
% governing equation, and we get a RHS contribution, AND an additional
% entry on the diagonal for u_p:
if ~isempty(bcsN)
    A(:, Nbcs+1:end) = A(:, Nbcs+1:end) + diag(sum(A(:, (length(bcsD)+1):Nbcs),2));
    
    % This step requires correction when dx != dy
    A(:, (length(bcsD)+1):Nbcs) = A(:, (length(bcsD)+1):Nbcs)*dx;
    
    %Then we add the boundary condition contributions to the RHS vector
    RHS = RHS - A(:, (length(bcsD)+1):Nbcs) * bcsN;
    %NOTE, we could combine the above two steps (i.e
    % RHS = RHS - A( : ,1:Nbcs) * [bcsD, bcsN]';
end;

% Remove Singularity
if (nargin == 9)&&(stagr == 'P')&&(Rm_Singularity)&&(isempty(bcsD))
    A(end, :) = 0; A(end, end) = 1; RHS(end) = 0;
end;

%% Finally, distribute entries of A to outputs

Arhs = A(:, 1:Nbcs);
Af = A(:, Nbcs+1:end);
