function [I,abserr]=cubature_loganl(intfcn, tol, starting_degree, k)

%------------------------------------------------------------------------------------
% OBJECT:
%
% COMPUTATION OF INTEGRAL OVER AN ANNULAR SECTOR BY A GAUSS LEGENDRE RULE ON THE RADIAL 
% PART AND ON THE ANGLE COMPONENT. BASIC ADAPTIVE PROCEDURE.
%
% INPUT:
%
% [intfcn]   : IT DEFINES THE FUNCTION TO INTEGRATE. IT WORKS AS "dblquad".
%
% [tol]      : TOLERANCE. RELATIVE ERROR.
%
% [k]        : GRADO DI PRECISIONE
%
% [n]        : NUMERO DEI POLINOMI DI LOGAN-SHEPP
%
% OUTPUT:
%
% [I]        : INTEGRAL ON THE ANNULAR SECTOR OF THE FUNCTION SPECIFIED BY "intfcn".
%
% {abserr]   : HISTORY OF THE ABSOLUTE ERROR.
%
%------------------------------------------------------------------------------------


thetamin=0; thetamax=2*pi; rmin=0; rmax=1;

f = fcnchk(intfcn); err=[]; warn=0; max_level=9; 

degree_GL_loc=starting_degree; 

I_loc(1)=integration_annulus_loc(f, thetamin, thetamax, rmin, rmax, degree_GL_loc);
abserr=[];

for index=2:max_level
    
    degree_GL_loc=2*degree_GL_loc; 
    I_loc(index)=integration_annulus_loc(f, thetamin, thetamax, rmin, rmax, degree_GL_loc); 
    
    den=abs(I_loc(index)); 
    abserr_loc=abs(I_loc(index)-I_loc(index-1));
    abserr=[abserr; abserr_loc];
    
    fprintf('\n \t [Degree]: %4.0f [Value]: %5.15f [Abs. Error]: %1.1e',degree_GL_loc, I_loc(index),abserr_loc)
    
    if abserr_loc < tol
        I=I_loc(end); 
        return;
    end
    
end

I=I_loc(end); warn=1;



            


%----------------------------------------------------------------
% "step_error".
%----------------------------------------------------------------

function err=step_error(I_loc,err)

if (I_loc(end) == 0)
    err(end+1)=abs( I_loc(end)-I_loc(end-1) );
else
    err(end+1)=abs( I_loc(end)-I_loc(end-1) )/abs(I_loc(end));
end

function I=integration_annulus_loc(f, thetamin, thetamax, rmin, rmax, k)

%------------------------------------------------------------------------------------
% OBJECT:
%
% COMPUTATION OF INTEGRAL ON AN ANNULAR SECTOR BY A GAUSS LEGENDRE RULE ON THE RADIAL 
% PART AND ON THE ANGLE COMPONENT.
%
% INPUT:
%
% [f]        : IT DEFINES THE FUNCTION TO INTEGRATE. 
%
% [thetamin] : MINIMUM ANGLE OF THE ANNULAR SECTOR.
%
% [thetamax] : MAXIMUM ANGLE OF THE ANNULAR SECTOR.
%
% [rmin]     : MINIMUM RADIUS OF THE ANNULAR SECTOR.
%
% [rmax]     : MAXIMUM RADIUS OF THE ANNULAR SECTOR.
%
% [n]: DEGREE OF CUBE RULE IN RADIAL AND ANGLE COMPONENT.
%
% OUTPUT:
%
% [I]        : INTEGRAL ON THE ANNULAR SECTOR OF THE FUNCTION SPECIFIED BY "intfcn".
%
%------------------------------------------------------------------------------------

[x,w]=quadrature_rules_1D(k); % 1D CUB RULE (SCALED IN [0,1])

r=rmin+(rmax-rmin)*x;
r_weights=(rmax-rmin)*w;

theta=thetamin+(thetamax-thetamin)*x;
theta_weights=(thetamax-thetamin)*w;

Lr=k+1;
Ltheta=Lr;

% THIS PART IS COMPLICATED TO BE READ. DATA ARE STORED SO TO 
% AVOID "FOR" CYCLES. 

rep_theta=repmat(theta,Lr,1); 
rep_theta=rep_theta(:);
rep_cos_theta=cos(rep_theta);
rep_sin_theta=sin(rep_theta);

rep_r=repmat(r',Ltheta,1);
rep_r=rep_r(:);

x=rep_r.*rep_cos_theta;
y=rep_r.*rep_sin_theta;

% f_pts=feval(f,x,y,varargin{:}); 
f_pts=feval(f,x,y);
 
integrand=rep_r.*f_pts;       % INTEGRAND IN POLAR COORDINATES.
integrand=reshape(integrand,Lr,Ltheta);

I_inner=r_weights'*integrand; % INNER INTEGRAL IN "[0,1]", I.E.
                              % RADIAL COMPONENT.

I=I_inner*theta_weights;  


%----------------------------------------------------------------------
% Some quadrature_rules for the interval "[0,1]".
%----------------------------------------------------------------------

function [nodes,weights]=quadrature_rules_1D(n)


N=[1:2:n-1]'; l=length(N); m=n-l; K=[0:m-1]';

beta=0.5./sqrt(1-(2*(1:n)).^(-2));
T=diag(beta,1)+diag(beta,-1);
[V,D]=eig(T);
x=diag(D); [x,ii]=sort(x); x=x';
w=2*V(1,ii).^2; 
nodes=x';
weights=w';

nodes=0.5*(nodes+1);
weights=real(weights)*0.5;