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
% Last Modified by Jing Lin in Apr. 2014


%% Fix boundary coordinates in X, and Y

dx = X(1,2) - X(1,1);
dy = Y(1,1) - Y(2,1);

% Attach the four external boundaries to X and Y
X = [X(:,1) - dx, X, X(:,end) + dx];
X = [X(1,:); X; X(end,:)];
Y = [Y(1,:) + dy; Y; Y(end,:) - dy];
Y = [Y(:,1), Y, Y(:,end)];

% Fix the coordinates of boundary cells adjacent to internal cells
Nbcs = length(bcsD) + length(bcsN);
if (lower(type) ~= 'u')
    % Extract the boundary cells with east adjacent internal cells
    ids = (Node<=Nbcs)&[(Node(:,2:end)>Nbcs), zeros(size(X,1),1)]; 
    X(ids) = X(ids) + dx/2;
    % Extract the boundary cells with west adjacent internal cells
    ids = (Node<=Nbcs)&[zeros(size(X,1),1), (Node(:,1:end-1)>Nbcs)];
    X(ids) = X(ids) - dx/2;
end;
if (lower(type) ~= 'v')
    % Extract the boundary cells with south adjacent internal cells
    ids = (Node<=Nbcs)&[(Node(2:end,:)>Nbcs); zeros(1, size(Y,2))]; 
    Y(ids) = Y(ids) - dy/2;
    % Extract the boundary cells with north adjacent internal cells
    ids = (Node<=Nbcs)&[zeros(1, size(Y,2)); (Node(1:end-1,:)>Nbcs)];
    Y(ids) = Y(ids) + dy/2;
end;


%% Extract mean coordinates for each boundary ID

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

