function [Dpts, Opts, Npts] = get_bc_coord(type, X, Y, Node, bcsD, bcsO, bcsN)
%function [DXY OXY NXY] = get_bc_coord(type, X, Y, Node, bcsD, bcsO, bcsN)
%
% This function retrieves the coordinates associated with different
% boundary types.
%
% INPUTS:
%       type:   This is required because of grid staggering. Use:
%               type = 'u' for the u-velocity grid
%               type = 'v' for the v-velocity grid
%               type = 'P' for the Pressure or tracer (density) grid
%       X:      The x-coordinate grid (from meshdomain.m)
%       Y:      The y-coordinate grid (from meshdomain.m)
%       Node:   The Node-numbering matrix
%       bcD:    Boundary values for Dirichlet boundary conditions
%       bcO:    Boundary values for Open boundary conditions
%       bcN:    Boundary values for Neumann boundary conditions
%
% OUTPUTS:
%       XYD:    Coordinate matrix for Dirichlet boundaries
%       XYO:    Coordinate matrix for Open boundaries
%       XYN:    Coordinate matrix for Neumann boundaries
%
% Note, XYD(:, 1) contains the x-coordinates, and XYD(:, 2) contains the y-coordinates.
%
% See also: meshdomain.m
%
% Written by:   Matt Ueckermann

%% First we have to pad X, and Y
dx = abs(X(1,2) - X(1,1));
dy = abs(Y(2,1) - Y(1,1));
xoff = dx/2;
yoff = dy/2;
if lower(type) == 'u'
    xoff = dx;
end
if lower(type) == 'v'
    yoff = dy;
end


%%
X = [X(:, 1) - xoff, X, X(:, end) + xoff];
X = [X(1, :); X; X(end, :)];

Y = [Y(1, :) + yoff; Y; Y(end, :) - yoff];
Y = [Y(:, 1), Y, Y(:, end)];

ND = length(bcsD) - length(bcsO);
NO = length(bcsO);
NN = length(bcsN);

Dpts = zeros(ND, 2);
Opts = zeros(NO, 2);
Npts = zeros(NN, 2);

for i=1:ND
    ids = Node == i;
    Dpts(i, :) = [mean(X(ids(:))), mean(Y(ids(:)))];
end
for i=ND+1:ND+NO
    ids = Node == i;
    Opts(i - ND, :) = [mean(X(ids(:))), mean(Y(ids(:)))];
end
for i=ND+NO+1:ND+NO+NN
    ids = Node == i;
    Npts(i - ND - NO, :) = [mean(X(ids(:))), mean(Y(ids(:)))];
end

