function [err_vec_max, err_vec_L2] = ...
    Test_EvalLstSq_Multiple(indx_f,n_deg,n_quad)

% Test the program EvalLstSq_Multiple. As an example, run the test case
%       [err_vec_max, err_vec_L2] = Test_EvalLstSq_Multiple(10,30,30);
% This will produce the orthonormal polynomial approximation of degree 30
% of the function #10 over the unit ball. 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 30 is used.
%
% INPUT:
% indx_f -  This is the index of the function 'fn' that is to be 
%           approximated over the unit ball. 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.

% Obtain the Fourier coefficients
[cf, poly_val, X, Y, Z, weight, fcn_val] = ...
    LeastSqCoeff2(@fn,n_deg,n_quad);

% Evaluate the least squares approximation and produce the errors.
[nr,nc,nz]=size(X);
approx_val = zeros(nr,nc,nz);
err_vec_max = zeros(n_deg,1);
err_vec_L2 = zeros(n_deg,1);

approx_val = cf(1).*poly_val(:,:,:,1);

for i=1:n_deg
    s1 = NthFunction(i-1,0,i-1)+1;
    s2 = NthFunction(i,0,i);
    for j= s1:s2
        approx_val = approx_val + cf(j).*poly_val(:,:,:,j);
    end
    err = fcn_val - approx_val;
    err_vec_max(i) = max(max(max(abs(err))));
    err_vec_L2(i) = sqrt(sum(sum(sum(err.*err.*weight))));
end
SqrtVolume = sqrt(4*pi/3);
err_vec_L2 = err_vec_L2/SqrtVolume;

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

figure
semilogy(deg,err_vec_max,'b-', deg, err_vec_L2, 'g--')
xlabel('degree n')
ylabel('error')
legend('max error', 'L^2 error','Location','NE')

% =====================================================================

function fval = fn(x,y,z)

% 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;   
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(10*x.*y.^2)).*(1+x)./(1+x.^2+y.^2);
    
case 10
    fval = ((1+x+y)./(1+x.^2 + y.^2 +z.^2)).*cos(6*x.*y.^2.*z);
 
end

end % fn

end % Test_EvalLstSq_Multiple


