function fig6_16

%  compute integral using adaptive simpson's rule
%  and plot points used for final computation
%  also: a comparison is made with MATLAB's quad fucntion
%
%  Inputs:
%        a = left endpoint of interval
%        b = right endpoint of interval
%        n = starting number of subintervals
%

global ifun ifun2

%%%% y = 1+(gg/pi)/((x-x0)^2+gg^2);
%%%% tolerances selected to get total error of 1e-6
a=0;
b=4;
exact=(atan(6*pi)+4*pi+atan(2*pi))/pi;
tol=5e-5;
tolQ=6e-6;  %%% error for quad command

%%%% y = cos(x.^3).^200;
%%%% tolerances selected to get total error of 1e-8
% a=0;
% b=3;
% exact=0.531594451912990;
% tol=3e-8;
% tolQ=2e-9;

%%%% y = exp(2*x)*sin(6*x);
%%%% tolerances selected to get total error of 1e-6
% a=0;
% b=2;
% exact=-8.225724190775976;
% tol=1e-4;
% tolQ=1e-5;


% L = max number of allowed levels
L=25;

%%% ifun = number of function evaluations
ifun=0;

tic
%  F = nf x 7 matrix, where rows contain info on subintervals that failed at previous level
%  A = rows are subintervals that are done, with last column containing the integral
%  nf = number of subintervals that failed at previous level
nf=1;
m=0.5*(a+b);
fm=f(m);
fa=f(a);  fb=f(b);
s=simp(a,m,b,fa,fm,fb);
F=[a m b fa fm fb s];
A=[];
Lmax=L;
%  this is the loop for the levels
for il=2:L
    FF=[];
    %%% use simpson on each remaining subinterval
    for is=1:nf
        mL=0.5*(F(is,1)+F(is,2));
        fmL=f(mL);
        mR=0.5*(F(is,2)+F(is,3));
        fmR=f(mR);
        s2L=simp(F(is,1),mL,F(is,2),F(is,4),fmL,F(is,5));
        s2R=simp(F(is,2),mR,F(is,3),F(is,5),fmR,F(is,6));
        s2=s2L+s2R;
        Err=abs(s2-F(is,7))/15;
        tolA=tol*(F(is,3)-F(is,1))/(b-a);
        %%% determine if the accuracy is attained
        if Err<tolA
            A=[A; F(is,1) F(is,3) (16*s2-F(is,7))/15];
        else
            FF=[FF; F(is,1) mL F(is,2) F(is,4) fmL F(is,5) s2L];
            FF=[FF; F(is,2) mR F(is,3) F(is,5) fmR F(is,6) s2R];
        end
    end
    nn=size(FF);
    %%% if no remaining subintervals, then leave the for-loop
    nf=nn(1);
    if nf==0
        Lmax=il;
        break
    end
    F=FF;
end

%%% if max number of levels exceeded, just use what was calcuated
if nf>0
    for is=1:nf
        A=[A; FF(is,1) FF(is,3) FF(is,7)];
    end
end

%%% add up values from each subinterval
ss=size(A)
sum=0;
for is=1:ss(1)
    sum=sum+A(is,3);
end
toc

% number of levels and subintevals used
levels_used=Lmax
subintervals=ss(1)

% value for integral, error, and total number of function evaluations
integral=sum
error=sum-exact
total_fns=ifun

% compare with MATLAB's quad routine
tic
ifun2=0;
quad_error=quad(@ff,a,b,tolQ)-exact
%quad(@ff,a,b,1e-8)
quad_fns=ifun2
toc

%%% plot function and subintervals
clf
% get(gcf)
set(gcf,'Position', [25 1115 658 230])
subaxis(1,1,1,1,'MT',0.01,'MB',0.18,'MR',-0.02,'ML',0.02,'P',0.04)
co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)

nx=200;
x=linspace(a-0.1,b+0.1,nx);
for ix=1:nx
    y(ix)=f(x(ix));
end

plot(x,y,'r','LineWidth',2)
hold on

for is=1:ss(1)
    xx=A(is,1); plot([xx xx],[0 f(xx)],'b','LineWidth',1.3)
    xx=0.5*(A(is,1)+A(is,2)); plot([xx xx],[0 f(xx)],'--b','LineWidth',1.1)
end
xx=b; plot([xx xx],[0 f(xx)],'b','LineWidth',1.3)

xlabel('x-axis')
ylabel('y-axis')
axis([x(1) x(nx) 0 3]);
set(gca,'FontSize',16,'FontWeight','bold')

%exportgraphics(gca,'/Users/mark/Desktop/aq8.eps')


function g=simp(a,m,b,fa,fm,fb)
g=(b-a)*(fa+4*fm+fb)/6;

function y=f(x)
global ifun
ifun=ifun+1;
%%%% 1 + Lorentzian Function
g=0.5*2/pi;  x0=3; gg=g/2;
y=1+(gg/pi)/((x-x0)^2+gg^2);
%y=cos(x^3)^200;
%y=50/(pi*(1+2500*x^2));
%y=exp(2*x)*sin(6*x);

% this is used for the MATLAB quad command
function y=ff(x)
global ifun2
ifun2 = ifun2 + length(x);
%%%% 1 + Lorentzian Function
g=0.5*2/pi;  x0=3; gg=g/2;
y=1+(gg/pi)./((x-x0).^2+gg^2);
%y=cos(x.^3).^200;
% y=50/(pi*(1+2500*x.^2));
%y=exp(2*x).*sin(6*x);















