
function demo_product_gaussian_rule(example)

% Demo: see
% K.Atkinson, W.Han, Spherical Harmonics and Approximations on the Unit
% Sphere: An introduction, p.172.
if nargin < 1, example=3; end
nv=2:6;
I=[];

switch example
    case 1
        f=@(x,y,z) exp(x)+0*y+0*z;
        I=1.476801374576529e+01;
    case 2
        f=@(x,y,z) exp(-x.^2-0.1*y.^2-2*z.^2);
        I=5.028153009823267;
     case 3
        f=@(x,y,z) (0.1*x+y.^2+0.5*z.^3).^3;
        I=1.969928193450972;
end

if isempty(I),  F=spherefun(f); I=sum2(F); format long; I, end
AE=[]; RE=[];

for n=nv
    [nodes,w]=gaussian_product_rule(n);
    x=nodes(:,1); y=nodes(:,2); z=nodes(:,3);
    fnodes=feval(f,x,y,z);

    In=w'*fnodes;
    AE(end+1)=abs(In-I);
    RE(end+1)=abs(In-I)/(abs(I)+(I == 0));

    fprintf('\n \t n: %3.0f nodes: %5.0f AE: %1.2e RE: %1.2e',...
        n,length(x),AE(end),RE(end))

end

fprintf('\n \n');

% sphere(20);
% hold on;
% plot3(x,y,z,'o','MarkerEdgeColor','k',...
%                        'MarkerFaceColor','g',...
%                        'MarkerSize',10)
% hold off;