warning off;

c = [1,10];
for j = 1:size(c,2)
cc = c(j)

for f_type = 1:4
    f = 0;
    I = [];
    x=[pi/6 pi/3; 0 2*pi];
    switch f_type
        case 1 
        % HESSE-WOMERSLEY TEST. (10^(-14) AT DEGREE 5)
            s='1./(cc*z).^2';
            f=inline(s,'x','y','z');
            ft=inline('1./(cc*sin(t)).^2.*(2+cos(t))','t','p');
        case 2 
        % EXPONENTIAL TEST. (10^(-16) AT DEGREE 30)
            s='exp(-x.^2-cc*y.^2-0.5*z.^2)';
            f=inline(s,'x','y','z');
            ft=inline('exp(-((2+cos(t)).*cos(p)).^2-cc*((2+cos(t)).*sin(p)).^2-0.5*(sin(t)).^2).*(2+cos(t))','t','p');
        % HUGE Error!
        case 3 
        %TRIGONOMETRIC TEST. (10^(-15) AT DEGREE 50)
            s='sin(-x.^2-cc*y.^2-0.5*z.^2)';
            f=inline(s,'x','y','z');
            ft=inline('sin(-((2+cos(t)).*cos(p)).^2-cc*((2+cos(t)).*sin(p)).^2-0.5*(sin(t)).^2).*(2+cos(t))','t','p');
        case 4 
        % HESSE-WOMERSLEY TEST. (SLOW CONVERGENCE. ODD BEH.)
            s='max(1-((x-3).^2+(y).^2+(z).^2),0)';
            f=inline(s,'x','y','z');
            ft=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');
    end

    Idbl=dblquad(ft,x(1,1),x(1,2),x(2,1),x(2,2))
    for i=5:5:100
        [tw,xw] = cubature_2_torus(x,i);
        fx = f(xw(:,1),xw(:,2),xw(:,3));
        i;
        res = xw(:,4)'*fx;
        err_res = abs(res - Idbl);
        I(i/5) = res;
        err(i/5) = err_res;
        printf('%d & %d & %1.2e & %1.2e\\\\ \n',cc,i,res,err_res);
    end

end
end
