% function [area,x,h,n,n_max] = adaptive_trap(a,b,tol) 
% 
% this does an adaptive simpsons rule.  It takes endpoints a and b,
% and the tolerance tol.  It gives back an area, the mesh-points, x, the
% interval lengths, h, and the number of intervals n.  It also gives
% the number of intervals you'd have to use if it was done with
% uniform stepping.

function [area,x,h,i,n_max] = adaptive_trap(a,b,tol)

b_orig = b;
a_orig = a;
i = 1;
area = 0;

% trapezoidal rule is h (f(a)/2 + f(b)/2) where b-a = h.

while (b_orig-a) > tol

  c = (a+b)/2;
  trap_coarse = (f(a)/2 + f(b)/2)*(b-a);
  trap_fine = (f(a)/2 + f(c) + f(b)/2)*(c-a);
  test = abs(trap_fine-trap_coarse)/3;

% if test > tol/2 then want to refine
  while test > tol/2
% refining means that I move the right hand point b to the left.  So
% I take the new right-hand point to be c.
    b = c;
    c = (b+a)/2;

    trap_coarse = (f(a)/2 + f(b)/2)*(b-a);
    trap_fine = (f(a)/2 + f(c) + f(b)/2)*(c-a);
    test = abs(trap_fine-trap_coarse)/3;
  end

% having gotten out of that while loop means that no more refinement
% was needed.  This means that I update the area, and I move the
% left hand point a to the right, and keep going

  area = area + trap_fine;
  x(i) = a;
  i = i+1;
% make the new left-hand point the old right-hand point
  a = b;
% make the new right-hand point the original right-hand point
  b = b_orig;
end

% n is the number of mesh-points
n = i-1;
for i=1:n-1
 h(i) = x(i+1)-x(i);
end

n_max = ceil((b_orig-a_orig)/min(h));