warning off;

% Demo of doubly_adaptive_optimized routine
% showcasing calls and different parameters

% First, define some test functions
for f_type = 1:1
    f = 0;
    I = [];
%    x=[pi/8 4*pi/5; -pi/6 2*pi/5];
    x=[0 pi; 0 pi];
    switch f_type
        case 1 
        % HESSE-WOMERSLEY TEST. (10^(-14) AT DEGREE 5)
            s='1./((x-2.5).^2+(y-2.5).^2+(z-6.5).^2).^1';
            f=inline(s,'x','y','z');
        case 2 
        % HESSE-WOMERSLEY TEST. (eps AT DEGREE 5)
            s='25*max(x+0*y+0*z,0).^3';
            f=inline(s,'x','y','z');
        case 3 
        % HESSE-WOMERSLEY TEST. (SLOW CONVERGENCE. ODD BEH.)
            s='(max(1/4-((x-1/sqrt(18)).^2+(y-1/sqrt(18)).^2+(z-4/sqrt(18)).^2),0)).^1';
            f=inline(s,'x','y','z');
        case 4 
        % POLYNOMIAL TEST. (10^(-16) AT DEGREE 20)
            s='3*x.^40+4*y.^40+100*z.^39';
            f=inline(s,'x','y','z');
        % BIG error
        case 5 
        % EXPONENTIAL TEST. (10^(-16) AT DEGREE 30)
            s='exp(-x.^2-100*y.^2-0.5*z.^2)';
            f=inline(s,'x','y','z');
        % HUGE Error!
        case 6 
        %TRIGONOMETRIC TEST. (10^(-15) AT DEGREE 50)
            s='sin(-x.^2-10*y.^2-0.5*z.^2)';
            f=inline(s,'x','y','z');
        case 7 
        % MODIFICATION HESSE-WOMERSLEY TEST. (10^(-15) AT DEGREE 15)
            s='1./((x-1).^2+(y-0).^2+(z-0).^2).^1';
            f=inline(s,'x','y','z');
        case 8 
        % HESSE-WOMERSLEY TEST. (10^(-14) AT DEGREE 5)
            s='25*max(x+0*y+0*z,0).^1';
            f=inline(s,'x','y','z');
        case 9 
        % MODIFICATION HESSE-WOMERSLEY TEST. (SLOW CONVERGENCE. GOOD BEH.)
            s='(max(1/4-((x-2/sqrt(5)).^2+(y-2/sqrt(5)).^2+(z-1/sqrt(5)).^2),0)).^2';
            f=inline(s,'x','y','z');
        case 10 
        % MODIFICATION HESSE-WOMERSLEY TEST. (MID CONVERGENCE. GOOD BEH.)
            s='(max(1/4-((x-1/sqrt(5)).^2+(y-2/sqrt(5)).^2+(z-2/sqrt(5)).^2),0)).^3';
            f=inline(s,'x','y','z');
    end

    % Basic call, all default parameters
    I = doubly_adaptive_optimized(f,x);

    % We can specify a tolerance (lax in this case)
    I = [I; doubly_adaptive_optimized(f,x,'tol',10^(-5))];

    % Or the error estimation method (Kronrod in this case)
    I = [I; doubly_adaptive_optimized(f,x,'err_method',3)];

    % Or the hint parameter 
    I = [I; doubly_adaptive_optimized(f,x,'f',1)];

    % Or the maximum order of the formulas
    I = [I; doubly_adaptive_optimized(f,x,'maxord',50)];

    % Or the maximum number of domain subdivisions
    %I = [I; doubly_adaptive_optimized(f,x,'maxdiv',10)];

    % And of course combine these parameters
    I = [I; doubly_adaptive_optimized(f,x,'maxdiv',5,'maxord',20,'f',0.1,'tol',10^(-5),'err_method',2)];

    % Notice whether we get different behaviours
    I

    % Integral on the sphere
    Is = doubly_adaptive_optimized(f,x,'tol',10^(-5),'domain',2)

    % Integrand can be given in polar or cartesian coordinates
    s='sin(t)';
    f=inline(s,'t','p');
    Ip = doubly_adaptive_optimized(f,x,'coord',1)

    s='z';
    f=inline(s,'x','y','z');
    Ic = doubly_adaptive_optimized(f,x,'coord',0)

    % And compare the performance to the standard algorithm
    I2 = doubly_adaptive_torus(f,x,10^(-10))
    err12 = abs(I(end) - I2)/abs(I2)

end
