warning off;

tols= [10^(-5) 10^(-8) 10^(-10) 10^(-15)];
x=[pi/6 pi/3; 0 pi/2];

ft{1}=inline('sin(-((2+cos(t)).*cos(p)).^2-10*((2+cos(t)).*sin(p)).^2-0.5*(sin(t)).^2).*(2+cos(t))','t','p');
ft{2}=inline('max(1-(((2+cos(t)).*cos(p)-3).^2+((2+cos(t)).*sin(p)).^2+(sin(t)).^2),0).*(2+cos(t))','t','p');
s{1}='sin(-x.^2-10*y.^2-0.5*z.^2)';
f{1}=inline(s{1},'x','y','z');
s{2}='max(1-((x-3).^2+(y).^2+(z).^2),0)';
f{2}=inline(s{2},'x','y','z');

for ii=1:1
    Idbl=dblquad(ft{ii},x(1,1),x(1,2),x(2,1),x(2,2))
    errm{1} ='higher';
    errm{2} ='anti';
    errm{3} ='kronrod';
    for err=1:3

        for t = 1:length(tols)
            tol = tols(t);

            [Iset,Eset]=set_adaptive_opt_alldeg(f{ii},x,'tol',tol,'err_method',err);
            printf('%s & %s & %1.1e & %1.2e & %1.2e & %1.2e & %1.2e \\\\ \n','set',errm{err},tol,Idbl,Iset,abs(Idbl-Iset),Eset);
            [Idoubly,Edoubly]=doubly_adaptive_optimized(f{ii},x,'tol',tol,'err_method',err);
            printf('%s & %s & %1.1e & %1.2e & %1.2e & %1.2e & %1.2e \\\\ \n','doubly',errm{err},tol,Idbl,Idoubly,abs(Idoubly-Idbl),Edoubly);
        end
    end
end
