% function [area,x,h,n,n_max] = adaptive_simpsons(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 itnervals you'd have to use if it was done with
% uniform stepping.

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

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

% simpsons rule is h/3 (f(a) + 4 f(c) + f(b)) where c-a = b-c = h.

while (b_orig-a) > tol

  c = (a+b)/2;
  h = c-a;

  simp_coarse = (f(a) + 4*f(c) + f(b))*h/3;

  h = h/2;
  c_left = (a+c)/2;
  c_right = (b+c)/2;

  simp_fine = (f(a) + 4*f(c_left) + f(c))*h/3;
  simp_fine = simp_fine + (f(c) + 4*f(c_right) + f(b))*h/3;

  test = abs(simp_fine-simp_coarse)/15;

% 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;
    h = c-a;

    simp_coarse = (f(a) + 4*f(c) + f(b))*h/3;

    h = h/2;
    c_left = (a+c)/2;
    c_right = (b+c)/2;

    simp_fine = (f(a) + 4*f(c_left) + f(c))*h/3;
    simp_fine = simp_fine + (f(c) + 4*f(c_right) + f(b))*h/3;

    test = abs(simp_fine-simp_coarse)/15;
  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 + simp_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));