function CofRT0 = Hdiv_RectRT0_NmlzBas_ProjVecFxn(VecFxn,RectMesh,GAUSSQUAD)
%% "Hdiv" local RT[0] on a rectangular mesh: Normalized basis 
% Computing the projection (into RT[0]) coeff. of a given vector function/field
% James Liu, ColoState; 2012/07--2017/02 

%% Mesh info 
NumEms = RectMesh.NumEms;
k1 = RectMesh.elem(:,1);  k2 = RectMesh.elem(:,3);
x1 = RectMesh.node(k1,1);  y1 = RectMesh.node(k1,2); 
x2 = RectMesh.node(k2,1);  y2 = RectMesh.node(k2,2);
Deltax = x2 - x1;
Deltay = y2 - y1;
area = Deltax *Deltay;
xc = 0.5*(x1+x2);
yc = 0.5*(y1+y2);

%% Computing the inner products / integrals 
RHS = zeros(NumEms,4);
NumQuadPts = size(GAUSSQUAD.RECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
qp = zeros(NumEms,2);
for k=1:NumQuadPts
  qp(:,1) = GAUSSQUAD.RECT(k,1)*x1 + GAUSSQUADRECT1(k,1)*x2;
  qp(:,2) = GAUSSQUAD.RECT(k,2)*y1 + GAUSSQUADRECT1(k,2)*y2;
  vec = VecFxn(qp);
  X = qp(:,1) - xc;  
  Y = qp(:,2) - yc;
  RHS(:,1) = RHS(:,1) + GAUSSQUAD.RECT(k,3)*vec(:,1);
  RHS(:,2) = RHS(:,2) + GAUSSQUAD.RECT(k,3)*vec(:,2);
  RHS(:,3) = RHS(:,3) + GAUSSQUAD.RECT(k,3)*(vec(:,1).*X);
  RHS(:,4) = RHS(:,4) + GAUSSQUAD.RECT(k,3)*(vec(:,2).*Y);
end
RHS = RHS.*area;

%% Computing the Gram matrix for local RT[0] nmlz.bas.
% GMRT0 has size NumEms*4 
GMRT0 = Hdiv_RectRT0_NmlzBas_GramMat(RectMesh);

%% Comuting the projection coefficients in local RT[0] 
CofRT0 = RHS./GMRT0;

return;