function [JSf,nodes]=simpsonadatt(a,b,f,tol,hmin)
%SIMPSONADATT Integrazione di Simpson adattiva
%   [JSF]=SIMPSONADATT(A,B,F,TOL,HMIN) calcola
%   l'integrale fra A e B della funzione F a 
%   meno di una tolleranza TOL con la formula
%   di Simpson adattiva. HMIN e' il minimo passo
%   di integrazione accettato.
%   
%   [JSF,NODES]=SIMPSONADATT(A,B,F,TOL,HMIN) restituisce
%   la distribuzione degli estremi degli intervalli di
%   integrazione utilizzati.
%

%   Saleri Fausto 02-02-01

A=[a,b]; N=[]; S=[]; JSf = 0; ba = b - a; nodes=[];

while ~isempty(A)
     [deltaI,ISc]=caldeltai(A,f);
     if abs(deltaI) <= 15*tol*(A(2)-A(1))/ba;
       JSf = JSf + ISc;
       S = union(S,A); nodes = [nodes, A(1) (A(1)+A(2))*0.5 A(2)]; 
       S = [S(1), S(end)];
       A = N; N = [];
     elseif A(2)-A(1) < hmin
       JSf = JSf + ISc;
       S = union(S,A); S = [S(1), S(end)];
       A = N; N = [];
       disp('Warning. L''ampiezza di A e'' minore di hmin');
     else
       Am = (A(1)+A(2))*0.5;
       A = [A(1) Am]; N = [Am, b];
     end
end
nodes=unique(nodes);

return

function [deltaI,ISc]=caldeltai(A,f)
%CALDELTAI Calcolo dello stimatore dell'errore per
%   la formula di Simpson adattiva.
%
L=A(2)-A(1);
t=[0; 0.25; 0.5; 0.5; 0.75; 1];
x=L*t+A(1); L=L/6;
w=[1; 4; 1];
fx=eval(f);
IS=L*sum(fx([1 3 6]).*w); 
ISc=0.5*L*sum(fx.*[w;w]);
deltaI=IS-ISc;

return
