function GaussQuadRect = SetGaussQuadRect(ch2x,ch2y)
%% Setting up a Gaussian quadrature for rectangles 
% Note: ch2x<=10;  ch2y<=10 
% James Liu, Rachel Cali, Zhuoran Wang, ColoState; 2007/01--2017/02 

%% Format: barycentric coordinates and weights 
% For rectangles:  alphax, alphay, weight 

% Choice (ch2x,ch2y) for a rectangle with opposite vertices (x1,y1),(x2,y2)
% Usage: 
%   ptx = alphax*x1 + (1-alphax)*x2
%   pty = alphay*y1 + (1-alphay)*y2 

GaussQuadLineX = SetGaussQuadLine(ch2x);
GaussQuadLineY = SetGaussQuadLine(ch2y);

NumQuadPts = ch2x * ch2y;
GaussQuadRect = zeros(NumQuadPts,3);

%% Note: Order of quadrature points: 
% from left to right, from bottom to top 
k = 0;
for j=1:ch2y
  for i=1:ch2x
    k = k + 1;
    GaussQuadRect(k,1) = GaussQuadLineX(i,1);
    GaussQuadRect(k,2) = GaussQuadLineY(j,1);
    GaussQuadRect(k,3) = GaussQuadLineX(i,3) * GaussQuadLineY(j,3);
  end
end

return;