function [cf, poly_val, X, Y, Z, weight, fcn_val] = ...
    LeastSqCoeff2(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 'EvalPolys3D'. The 
% Fourier coefficients are calculated using numerical integration
% over the unit ball, and 'n_quad' is the quadrature parameter. That
% means we use Gauss-Jacobi quadrature in the radial direction (rho) with
% n_quad nodes, we use the trapezoidal rule in the azimuthal angular 
% (theta) direction with 2*n_quad+1 subdivisions, and we use Gauss-Legendre 
% quadrature in the polar angular direction (phi) with n_quad subdivisions.
% 0<= rho <= 1, 0<= theta < 2*Pi, 0<= phi < Pi.  See details in the 
% article.
%
% INPUT:
% fcn -     This is the function that is to be approximated over 
%           the unit ball. 
% 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:
% cf -       The approximate Fourier coefficients for all terms of
%            degree less than or equal to n_deg.
% poly_val - The function values of the orthonormal basis computed in
%            'EvalPolys3D'.
% (X,Y,Z) -  The points at which the orthonormal polynomials are evaluated.
% weight -   Combined weight calculated by the program 'GaussQuadrature'.
% fcn_val -  The function values of 'fcn' evaluated at (X,Y,Z).

% Produce the nodes and weights for the numerical integration.
nq = n_quad;
hq = pi/nq;
the = hq*(1:2*nq);
[r_nodes,r_wgts] = GaussQuadrature(0,2,nq);
[p_nodes,p_wgts] = GaussQuadrature(0,0,nq);
r_wgts = r_wgts/8;
p_wgts = p_wgts*hq;

[rweight,~,pweight]= meshgrid(r_wgts,the,p_wgts);
weight = rweight.*pweight;

r_nodes = (r_nodes+1)/2; p_nodes = acos(p_nodes);
[Rho,Theta, Phi] = meshgrid(r_nodes,the, p_nodes);
SinPhi = sin(Phi);

X = Rho.*SinPhi.*cos(Theta); Y = Rho.*SinPhi.*sin(Theta); Z=Rho.*cos(Phi);

% Produce the function values.
fcn_val = fcn(X,Y,Z); fcn_wgts = fcn_val.*weight;

% Produce the orthonormal polynomials of degrees 0 through n_deg.
poly_val = EvalPolys3D(X,Y,Z,n_deg);

% Calculate the Fourier coefficients using quadrature.
order = (n_deg+3)*(n_deg+2)*(n_deg+1)/6;
cf = zeros(order,1);

for n1=1: order
    cf(n1)=sum(sum(sum(fcn_wgts.*poly_val(:,:,:,n1))));
end

end % LeastSqCoeff2
