function [I,x] = adaptgausspts(f, a, b, tol)
% Calculate integral of f from a to b with max error of tol (default 1e-6).
if ( nargin < 4 ), tol = 1e-6; end % Default tolerance
[x,c] = kronpts15(a,b); ff = f(x); I = sum(c.*ff);
[~,c] = gausspts7(a,b); G7 = sum(c.*ff(2:2:end-1));
E = (200*abs(I - G7)).^(1.5);
if ( abs(E) > tol )   % Fail. Subdivide
    c = (a + b) / 2;
%     I = adaptgauss(f, a, c, tol/2) + ...
%         adaptgauss(f, c, b, tol/2);
    [IL, xL] = adaptgausspts(f, a, c, tol/2);
    [IR, xR] = adaptgausspts(f, c, b, tol/2);
    I = IL + IR;
    x = unique([xL ; xR]);
end
end

function [x, c] = gausspts7(a,b)
x = [  -0.949107912342759
       -0.741531185599394
       -0.405845151377397
        0
        0.405845151377397
        0.741531185599394
        0.949107912342759];
c = [   0.129484966168870
        0.279705391489277
        0.381830050505119
        0.417959183673469
        0.381830050505119
        0.279705391489277
        0.129484966168870];
scl = (b-a)/2;
x = scl*(x+1) + a;
c = scl*c;
end

function [x, c] = kronpts15(a,b)
x = [  -0.991455371120813
       -0.949107912342759
       -0.864864423359769
       -0.741531185599394
       -0.586087235467691
       -0.405845151377397
       -0.207784955007898
        0
        0.207784955007898
        0.405845151377397
        0.586087235467691
        0.741531185599394
        0.864864423359769
        0.949107912342759
        0.991455371120813];

c   = [ 0.022935322010529
        0.063092092629979
        0.104790010322250
        0.140653259715525
        0.169004726639267
        0.190350578064785
        0.204432940075298
        0.209482141084728
        0.204432940075298
        0.190350578064785
        0.169004726639267
        0.140653259715525
        0.104790010322250
        0.063092092629979
        0.022935322010529];
scl = (b-a)/2;
x = scl*(x+1) + a;
c = scl*c;
end




