function [AP, rhsP, APBcs] = FixDiffPObcs(AP, bcsOP, Nbcs, rhsP, NodeP, dx, dy, APBcs)
%function [AP rhsP APBcs] = FixDiffPObcs(AP,bcsOP,Nbcs,rhsP,NodeP,dx,dy,APBcs)
% This function fixes the A matrix and RHS vector for the discretized 2D
% Laplace equation using Finite Volume method on a staggerd uv grid for open 
% pressure boundary conditions.
%
% INPUTS:
%       AP:     Original A matrix for Pressure
%       bcsOP:  Vector containing the BC numbers associated with the OBC
%       Nbcs:    Total number of BCS.
%       rhsP:   The Right-hand-side vector from BC contributions for P
%       NodeP:  Is the node numbering matrix. size(NodeP) = [Ny+2,Nx+2]
%               NOTE: Node needs to include boundary conditions.
%       dx:     Size of CV in x direction
%       dy:     Size of CV in y direction
%       APBcs:   The matrix used to create the Right-hand-side vector for
%                  Pressure
% OUTPUTS:
%       AP:     The corrected A matrix for Pressure
%       rhsP:   The Right-hand-side vector from BC contributions for Pressure
%       APBcs:   The matrix used to create the Right-hand-side vector for Pressure
%
% Authors:   Matt Ueckermann, Pierre Lermusiaux, Pat Haley for MIT course 2.29

if nargin<8, APBcs=zeros(length(AP),length(Nbcs)); end

[Ny, Nx] = size(NodeP);

ids = find((NodeP >= bcsOP(1))&(NodeP <= bcsOP(end)));
% Find Ids of interior cells with a northern boundary
n = ids + 1;
n(n>Ny*Nx) = [];
n(NodeP(n)<=Nbcs) = [];
% Find Ids of points with a southern boundary
s = ids - 1;
s(s<1) = [];
s((NodeP(s)<=Nbcs)) = [];
% Find Ids of points with a western boundary
w = ids + Ny;
w(w>Ny*Nx) = [];
w(NodeP(w)<=Nbcs) = [];
% Find Ids of points with a eastern boundary
e = ids - Ny;
e(e<1)=[];
e(NodeP(e)<=Nbcs) = [];

% Re-Constructing the Inner part of the 'A' Matrix for north boundaries
for j = 1 : length(n)
    % Start with the middle node, subtract out old value, add new value
    Rowij = NodeP(n(j)) - Nbcs;
    AP(Rowij, Rowij) = AP(Rowij, Rowij) - 2/dy^2;
    % Now fix the south part in the computational cell
    colno = NodeP(n(j)+1) - Nbcs;
    AP(Rowij, colno) = AP(Rowij, colno) + 1/dy^2;
    rhsP(Rowij) = 0;
    APBcs(Rowij, NodeP(n(j)-1)) = 0;
end;
% Re-Constructing the Inner part of the 'A' Matrix for south boundaries
for j = 1 : length(s)
    % Start with the middle node, subtract out old value, add new value
    Rowij = NodeP(s(j)) - Nbcs;
    AP(Rowij, Rowij) = AP(Rowij, Rowij) - 2/dy^2;
    % Now fix the north part in the computational cell
    colno = NodeP(s(j)-1)-Nbcs;
    AP(Rowij, colno) = AP(Rowij, colno) + 1/dy^2;
    rhsP(Rowij) = 0;
    APBcs(Rowij, NodeP(s(j)+1)) = 0;
end;
% Re-Constructing the Inner part of the 'A' Matrix for west boundaries
for j = 1 : length(w)
    % Start with the middle node, subtract out old value, add new value
    Rowij = NodeP(w(j)) - Nbcs;
    AP(Rowij, Rowij) = AP(Rowij, Rowij) - 2/dx^2;
    % Now fix the east part in the computational cell
    colno = NodeP(w(j)+Ny) - Nbcs;
    AP(Rowij, colno) = AP(Rowij, colno) + 1/dx^2;
    rhsP(Rowij) = 0;
    APBcs(Rowij, NodeP(w(j)-Ny)) = 0;
end;
% Re-Constructing the Inner part of the 'A' Matrix for east boundaries
for j = 1 : length(e)
    % Start with the west node, subtract out old value, add new value
    Rowij = NodeP(e(j)) - Nbcs;
    AP(Rowij, Rowij) = AP(Rowij, Rowij) - 2/dx^2;
    % Now fix the north part in the computational cell
    colno = NodeP(e(j)-Ny) - Nbcs;
    AP(Rowij, colno) = AP(Rowij, colno) + 1/dx^2;
    rhsP(Rowij) = 0;
    APBcs(Rowij, NodeP(e(j)+Ny)) = 0;
end;

% Remove the singularity
AP(end,:) = 0; AP(end, end) = 1;
