function [err_max_vec,err_L2_vec] = ...
                     Test_EvalLstSq(indx_f,n_deg,n_quad,n_eval)

% Test the program EvalLstSq. As an example, run the test case
%       [err_max_vec,err_L2_vec] = Test_EvalLstSq(9,30,40,50);
% This will produce the orthonormal polynomial approximation of degree 30
% of the function #9 over the unit disk. In addition, it produces
% The maximum errors and the L^2 errors of the least squares approximations 
% of degrees 1 through 30. To produce the Fourier coefficients, numerical 
% integration with the parameter 40 is used. 
%
% Graph the approximating polynomial and its error for the degree 'n_deg' 
% case. Also graph degree vs. error on a semi-log scale. Display these
% errors as well.
%
% INPUT:
% indx_f -  This is the index of the function 'fn' that is to be 
%           approximated over the unit disk. See the function 'fn'  
%           given below.
% 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 for all polynomials of degree less than 2*n_quad+1.
% 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.
%
% OUTPUT:
% 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 the L2 error of the errors as specified above.

% Obtain the Fourier coefficients 
cf = LeastSqCoeff(@fn,n_deg,n_quad);

% Evaluate the least squares approximation and produce the errors.
[approx_val,error,err_max_vec,err_L2_vec,X,Y] = ...
                                    EvalLstSq(@fn,n_deg,cf,n_eval);
figure
mesh(X,Y,approx_val)
title(['degree = ',num2str(n_deg)])
xlabel('x'); ylabel('y');

figure
mesh(X,Y,error)
title(['Error:      n#quad = ',num2str(n_quad)])
xlabel('x'); ylabel('y');

%disp('[degree, Maximum errors, L2 errors]')
deg = (1:n_deg)';
%format long e
%disp([deg,err_max_vec,err_L2_vec])

figure
semilogy(deg,err_max_vec,deg,err_L2_vec,'--')
xlabel('degree n')
ylabel('error')
title(['Errors: n#quad = ',num2str(n_quad),...
                                   '    n#eval = ',num2str(n_eval)])
legend('max error','L^2 error')                               

function fval = fn(x,y)

% This is the function for which we are calculating the least squares
% approximation of degree 'n_deg'.

switch indx_f
case 1
    fval = ones(size(x));
case 2
    fval = x;
case 3
    fval = x.*y.^5;   
case 4
    fval = cos(x - 2*y);
case 5
    fval = exp(-(1-x).^2).*cos(x+y.^3);
case 6
    fval = abs(x) + abs(y);
case 7
    fval = (x.^2+y.^2).*(abs(x)+abs(y));
case 8
    a = .4; b = .8; c = pi;
    fval = exp(-(a*(1-x).^2+b*y.*x)).*cos(c*x.*y.^2);
case 9
    fval = (cos(6*x.*y.^2)).*(1+x)./(1+x.^2+y.^2);
end

end % fn

end % Test_EvalLstSq


