function adaptivequad(tol)
% demo adaptive quadrature with Clenshaw-Curtis panels. Barnett 5/20/13
% needs clencurt.m

% real analytic func but with very near singularity
a = 0.1; % singularity distance parameter (+-ia)
f = @(x) a./(a^2+x.^2); fint = 2*atan(1/a)

% Or, can even do unbounded (improper) integral:
%a = -0.3;  % location of sqrt singularity
%f = @(x) (x>a)./sqrt(x-a); fint = 2*sqrt(1-a)

%a = 100; f = @(x) cos(a*x); fint = (sin(a)-sin(-a))/a % oscillations

%f = @(x) sin(1/(x-1));  fint = 0;   % fails - runs forever

[I fevals] = adaptiveCCquad(f,-1,1, 1e-8)

fprintf('err = %.3g\n',I-fint)

end

function [I fevals] = adaptiveCCquad(fun,a,b,tol,s)
% recursive adaptive quadrature with Clenshaw-Curtis panels.
% s is a structure that contains nodes and weights for n=9,17.

n = 9; n2 = 17;              % CC orders for coarse and fine nodes
if nargin<5        % fill s structure with node and weights so don't recompute
  [s.x s.w]  = clencurt(n-1);
  [s.x2 s.w2]  = clencurt(n2-1);
end
mid = (a+b)/2; siz = (b-a)/2;
x2 = mid + siz*s.x2;         % scale fine nodes for this panel
w = siz*s.w; w2 = siz*s.w2;  % scale weights

f2 = fun(x2);  f = f2(1:2:end);      % get f at the 2 node sets (one is subset)
I = dot(f,w); I2 = dot(f2,w2);
if abs(I2-I) < tol             % criterion for splitting
  %plot(x2, w2/siz, 'r.'); hold on;   % show the nodes used
  plot(x2, f2, '.'); hold on;      % show the func at the nodes used
  plot([a a],[0 f2(1)],'-'); plot([b b],[0 f2(end)],'-'); 
  I = I2; fevals = n2;                   % use the most accurate
  return
else                         % split and recurse...
  [Il fevalsl] = adaptiveCCquad(fun,a,mid,tol,s);
  [Ir fevalsr] = adaptiveCCquad(fun,mid,b,tol,s);
  I = Il + Ir;
  fevals = fevalsl + fevalsr;  % keep track of how many func evals used
end
  
end
