% Here is an example script of using PMPack to generate quadrature rules. 

% The Legendre parameter corresponds to a uniform distribution over the
% interval [-1,1]. 
s_uniform = parameter('Legendre');

% You can generate a 1D Gaussian quadrature rule with n points.
n = 10;
[points, weights] = gaussian_quadrature(s_uniform, n);

% To approximate an integral with these points...
f = @(x) cos(2*pi*x);
Int_f = f(points)'*weights;

% To generate a tensor product quadrature rule in two variables...
s = [parameter('Legendre'); parameter('Legendre')]; 
n = [10; 5];
[points2d, weights] = gaussian_quadrature(s, n);

% To integrate a bivariate function with these points...
f = @(x,y) cos(2*pi*x) + sin(3*pi*y);
Int_f = 0;
for i=1:length(points)
    Int_f = Int_f + f(points2d(i,1),points2d(i,2))*weights(i);
end

% The code above can be vectorized like the 1D example for efficiency.
% You can generate a parameter corresponding to the Chebyshev measure 
% using...
s = parameter('Chebyshev');

% The point sets of this measure are nested. In other words, the points of
% the 2^n-1 point quadrature rule are a subset of the points in the 
% 2^(n+1)-1 point rule. 

