function FourierCoeff = LeastSqCoeff(fcn,n_deg,n_quad)

% Produce the coefficients for the least squares approximation
% of 'fcn' by a polynomial of degree less than or equal to 'n_deg'. 
% The orthonormal basis is that computed in 'EvalOrthoPolys'. The 
% Fourier coefficients are calculated using numerical integration
% over the unit disk, and 'n_quad' is the quadrature parameter.
% That means we use Gaussian quadrature in the radial direction with
% n_quad+1 nodes, and we use the trapezoidal rule in the angular 
% direction with 2*n_quad+1 subdivisions.
%
% INPUT:
% fcn -     This is the function that is to be approximated over 
%           the unit disk. 
% n_deg -   This is the degree of the polynomial approximation.
% n_quad -  This is the integration parameter for the numerical
%           integration of the Fourier coefficients. The quadrature
%           scheme is exact when applied to integrate any polynomial 
%           of degree less than 2*n_quad+1.
%
% OUTPUT:
% FourierCoeff - The approximate Fourier coefficients for all terms
%           of degree less than or equal to n_deg. 

% Produce the nodes and weights for the numerical integration.
nq = n_quad; hq = 2*pi/(2*nq+1);
theta = hq*(1:2*nq+1);
[r_nodes,r_wgts] = gl_nodes_wts(0,1,nq+1);
r_wgts = r_nodes.*r_wgts;
[R,Th] = meshgrid(r_nodes,theta);
X = R.*cos(Th); Y = R.*sin(Th);
[nr,nc] = size(X);
wgts = zeros(nr,nc);
for i=1:nr
    wgts(i,:) = hq*r_wgts';
end

% Produce the function values.
fcn_val = fcn(X,Y); fcn_wgt = fcn_val.*wgts;

% Produce the orthonormal polynomials of degrees 0 through n_deg.
[a_tot,c_tot,d_tot] = TripleRecurCoeff(n_deg-1);
poly_val = EvalOrthoPolys(X,Y,n_deg,a_tot,c_tot,d_tot);

% Calculate the Fourier coefficients using quadrature.
order = ((n_deg+1)*(n_deg+2))/2;
cf = zeros(order,1);
for j=1:order
    cf(j) = sum(sum(fcn_wgt.*poly_val(:,:,j)));
end

FourierCoeff = cf;

end % LeastSqCoeff
