function [dpdx dpdy] = Grad2D(P,NodeP,dx,dy,idsu,idsv)
% function [dpdx dpdy] = Grad2D(P,NodeP,dx,dy)
% This function takes the gradient of P=[P(NodeP)] at the INTERIOR cell centres
% of u and v respectively. 
%
% INPUTS:
%       P:      Vector containing values of Pressure
%       NodeP:  Reference node numbering matrix for Pressure
%       dx:     delta x, grid spacing in x direction
%       dy:     delta y, grid spacing in y direction
%       idsu:   Location of masked u-velocity nodes in interior
%       idsv:   Location of masked v-velocity nodes in interior
%       
% OUTPUTS:
%       dpdx:   dp/dx, x-derivative of Pressure returned at interior
%               u-velocity points
%       dpdy:   dp/dy, y-derivative of Pressure returned at interior
%               v-velocity points
%
% Authors:   Matt Ueckermann, Pierre Lermusiaux, Pat Haley for MIT course 2.29

dpdx = (P(NodeP(2:end-1  ,3:end-1))-P(NodeP(2:end-1,2:end-2)))/dx;
dpdy = (P(NodeP(2:end-2,2:end-1  ))-P(NodeP(3:end-1,2:end-1  )))/dy;

dpdx = dpdx(idsu);
dpdy = dpdy(idsv);
