function [approx_val,error,err_max_vec,err_L2_vec,x,y] = ...
                                 EvalLstSq(fcn,n_deg,coeff,n_eval)

% Evaluate all of the least square approximating polynomials of degree less 
% than or equal to 'n_deg', over the unit disk, whose coefficients with 
% respect to the standard orthonormal basis are given in 'coeff'. Then 
% compare to the true function 'fcn' being approximated by the polynomials.
%
% Calculate both the maximum error and the L^2 error for each degree. 
%
% INPUT:
% fcn -     This is the function to be approximated over the unit disk.
% n_deg -   This is the degree of the final approximation. This program
%           looks at the error of all least squares approximations of
%           degrees 1 through n_deg.
% coeff -   These are the coefficients of the least squares approximation,
%           as produced by 'LeastSqCoeff'.
% n_eval -  This parameter is used in defining the evaluation points (x,y)
%           for the evaluation of the least squares approximation. The
%           number of evaluation points is approximately 2*n_eval^2. See
%           below in the program for the definition of (x,y).
%
% OUTPUT:   
% approx_val - This is the least squares approximation of degree n_deg,
%           evaluated at the points (x,y).
% error -   The error in approx_val at the points (x,y).
% err_max_vec - This is the maximum of the absolute value of the errors at 
%           the evaluation node points defined in the arrays given as 
%           (X,Y) and produced in 'EvalLstSq', as below. This is calculated
%           for the least squares approximation of degrees 1 through n_deg.
% err_L2_vec - These are an approximation of the L^2 estimate of the errors
%           as specified above.
% (x,y) -   These are the evaluation points for the least squares
%           approximation.

% Determine the points at which to evaluate the polynomial.
nq = n_eval; % Parameter to set up the evaluation points. 
hq = pi/nq; theta = hq*(0:2*nq);
hr = 1/nq; r_nodes = hr*(0:nq);
[R,Th] = meshgrid(r_nodes,theta);
x = R.*cos(Th); y = R.*sin(Th);

% Evaluate the orthonormal polynomials using the triple recursion relation.
[a_tot,c_tot,d_tot] = TripleRecurCoeff(n_deg-1);
poly_val = EvalOrthoPolys(x,y,n_deg,a_tot,c_tot,d_tot);

fcn_val = fcn(x,y);

% Evaluate the approximating polynomials of degrees 1 through n_deg.
approx_val = coeff(1)*poly_val(:,:,1);
err_max_vec = zeros(n_deg,1); err_L2_vec = zeros(n_deg,1);
for i=1:n_deg
    iord = (i*(i+1))/2 + 1;
    for j = 0:i
        k = iord + j;
        approx_val = approx_val + coeff(k)*poly_val(:,:,k);
    end
    error = fcn_val - approx_val;
    % Calculate the maximum error
    err_max_vec(i) = max(max(abs(error)));
    % Calculate the L2 error
    sum_error = sum(sum(R.*error.^2));
    err_L2_vec(i) = sqrt(hq*hr*sum_error)/sqrt(pi);
end

end % EvalLstSq

