function [I,Ct] = adaptiveAvg(aa,bb,tol,ct,ld,sm,al,la,hh,wl,ws)
% Gauss-Legendre adaptive quadrature 
% using 4 & 6 points.
% See: Algorithm 8.1 on p.358 Heath
% fcn = (string) name of the function to be integrated.
% aa  = the lower limit.
% bb  = the upper limit.
% tol = the tolerance (largest possible absolute error).
% ct  = counter to keep track of the total number of times 
%       adaptiveQuadG2 calls itself. ct must be 0 initially.
% ld = center wavelength.
% sm = sigma of incident spectrum.
% al = number of sigma to include in integration.
% la = lattice constant.
% hh = etching depth.
% wl = width of the larger cavity
% ws = width of the smaller cavity.
% Functions needed: 
%   twoCavitiesAuSymmMDA4I.m twoCavitiesAuSymmMDA6I.m
%   4- & 6-pt Gauss-Legendre quadrature to compute avg I.
% K. Ming Leung, 08/20/03

I1 = twoCavitiesAuSymmMDA4I(aa,bb,ld,sm,al,la,hh,wl,ws);
I2 = twoCavitiesAuSymmMDA6I(aa,bb,ld,sm,al,la,hh,wl,ws);

m = aa + 0.5*(bb-aa);
ct = ct+1;  % increase counter Ct by 1.
% if m <= a | m >= b
%     error('Tolerance may not be met');
%     I = I2;
% end
if abs(I2-I1) < tol % works better in Matlab
    I = I2;
else
    [T1,Ct] = adaptiveAvg(aa,m,tol,ct,ld,sm,al,la,hh,wl,ws);
    [T2,Ct]=  adaptiveAvg(m,bb,tol,ct,ld,sm,al,la,hh,wl,ws);
    I = T1+T2;
end

