function [NodeP NodeU NodeV idsP idsu idsv] = NodePad(Mask, bnd, Masku, Maskv, periodicity)
% function [NodeP NodeU NodeV idsP idsu idsv] =...
%                           NodePad(Mask, bnd, Masku, Maskv, periodicity)
% This function adds 4 boundaries around the outside of the domain, and
% creates the Node numbering reference matrices.
%
% INPUTS:
%       Mask:   The tracer mask matrix (also the Pressure mask matrix)
%       bnd:    Vector containing true or false values of which boundaries
%               to add. bound = [left right bottom top]
%       Masku:  The u velocity mask matrix (optional)
%       Maskv:  The v velocity mask matrix (optional)
%       periodicity:  Switch to check for periodicty (optional)
%               NOTE: Can also specify this as third entry.
%   NOTE: All internal boundaries should be number from 1 to #internal BCs
% OUTPUTS:
%       NodeP:  Reference node numbering matrix for Pressure
%       NodeU:  Reference node numbering matrix for u-velocity 
%       NodeV:  Reference node numbering matrix for v-velocity 
%       idsP:   Location of active Pressure nodes in interior
%       idsu:   Location of active u-velocity nodes in interior
%       idsv:   Location of active v-velocity nodes in interior
%
% See also: maskuv.m
%
% Written by:   Matt Ueckermann

%Get mask size
[Ny Nx] = size(Mask);
%Get number of different internal bcs
bndnum  = max(Mask(:));
Nbcs    = max(Mask(:))+sum(bnd);
%Check inputs, and create u, v, masks if necessary
if nargin <3 
    [Masku, Maskv] = maskuv(Mask);
end
if nargin <5
    periodicity=false;
end;
if nargin == 3
   if length(Masku)<2
       periodicity = Masku;
       [Masku, Maskv] = maskuv(Mask);
   end
end
if isempty(Masku) || isempty(Maskv)
    [Masku, Maskv] = maskuv(Mask);
end
%Find internal active nodes id's.
idsP = find(Mask==0);
NodeP=Mask;
%Number only the internal active nodes, starting from #bcs+1
NodeP( idsP ) = [1:length(idsP)]+Nbcs;
%Do this same for U and V masks, if number of output arguments are greater
%than 1.
if nargout > 1
    idsu = find(Masku==0);
    NodeU=Masku;
    NodeU( idsu ) = [1:length(idsu)]+Nbcs;
    idsv = find(Maskv==0);
    NodeV=Maskv;
    NodeV( idsv ) = [1:length(idsv)]+Nbcs;
end
%Add boundaries
%Left
if bnd(1)
    bndnum = bndnum+1;
    NodeP = [bndnum*ones(1,Ny) ; NodeP']';
    if nargout >1
        NodeU  = [bndnum*ones(1,Ny  ) ; NodeU']';
        NodeV  = [bndnum*ones(1,Ny-1) ; NodeV']';
    end
    Nx = Nx + 1;
end%Right
if bnd(2)
    bndnum = bndnum+1;
    NodeP = [NodeP' ; bndnum*ones(1,Ny)]';
    if nargout >1
        NodeU  = [NodeU' ; bndnum*ones(1,Ny  )]';
        NodeV  = [NodeV' ; bndnum*ones(1,Ny-1)]';
    end
    Nx = Nx + 1;
end
%Bottom
if bnd(3)
    bndnum = bndnum+1;
    NodeP = [NodeP; bndnum*ones(1,Nx)];
    if nargout >1
        NodeU  = [NodeU; bndnum*ones(1,Nx-1)];
        NodeV  = [NodeV; bndnum*ones(1,Nx  )];
    end
end
%Top
if bnd(4)
    bndnum = bndnum+1;
    NodeP  = [bndnum*ones(1,Nx  ) ; NodeP];
    if nargout >1
        NodeU  = [bndnum*ones(1,Nx-1) ; NodeU];
        NodeV  = [bndnum*ones(1,Nx  ) ; NodeV];
    end
end
Nx = Nx - sum(bnd(1:2));
%% Check for periodicity
% Add a boundary to handle periodicity correctly.
if periodicity
    if  size(NodeP,1) == Ny
        NodeP = [NodeP(end,:) ; NodeP ; NodeP(1,:)];
        NodeU = [NodeU(end,:) ; NodeU ; NodeU(1,:)];
        NodeV = [NodeV(end,:) ; NodeV ; NodeV(1,:)];
    end
    if size(NodeP,2) == Nx
        NodeP = [NodeP(:,end)' ; NodeP' ; NodeP(:,1)']';
        NodeU = [NodeU(:,end)' ; NodeU' ; NodeU(:,1)']';
        NodeV = [NodeV(:,end)' ; NodeV' ; NodeV(:,1)']';
    end
end
