
function [I,abserr]=cubature_disk(intfcn, tol, cub_type, starting_degree, varargin)


%------------------------------------------------------------------------------------
% 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.
%
% [cub_type] : QUADRATURE TYPE.
%            [cub_type=1]: FEJER 1.
%            [cub_type=2]: FEJER 2.
%            [cub_type=3]: CLENSHAW-CURTIS.
%            [cub_type=4]: GAUSS-LEGENDRE.
%            [cub_type=5]: TRAPEZOIDAL RULE.
%
% OUTPUT:
%
% [I]        : INTEGRAL ON THE ANNULAR SECTOR OF THE FUNCTION SPECIFIED BY "intfcn".
%
% {abserr]   : HISTORY OF THE ABSOLUTE ERROR.
%
% [warn]     : IF [warn=0] THE METHOD ENDS NORMALLY. 
%              IF [warn=1] THE METHOD ENDS BECAUSE TOO MANY POINTS ARE NEEDED. 
%
%------------------------------------------------------------------------------------


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, cub_type, varargin{:});
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, cub_type, varargin{:}); 
    
    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, degree_GL, cub_type, varargin)

%------------------------------------------------------------------------------------
% 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.
%
% [degree_GL]: 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(degree_GL,cub_type); % 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=degree_GL+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,cubature_type)

% SEE WALDVOGEL PAPER. ADDED NODES

% Weights of the Fejer2, Clenshaw-Curtis and Fejer1 quadrature by DFTs
% n>1. Nodes: x_k = cos(k*pi/n)

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

switch cubature_type
    
case 1 % FEJER 1.
v0=[2*exp(i*pi*K/n)./(1-4*K.^2); zeros(l+1,1)];
v1=v0(1:end-1)+conj(v0(end:-1:2)); 
weights=ifft(v1);
k=(1/2):(n-(1/2)); nodes=(cos(k*pi/n))'; 

case 2 % FEJER 2.
v0=[2./N./(N-2); 1/N(end); zeros(m,1)];
v2=-v0(1:end-1)-v0(end:-1:2); 
wf2=ifft(v2); weights=[wf2;0];
k=0:n; nodes=(cos(k*pi/n))'; 

case 3 % CLENSHAW CURTIS.
g0=-ones(n,1); g0(1+l)=g0(1+l)+n; g0(1+m)=g0(1+m)+n;
g=g0/(n^2-1+mod(n,2));
v0=[2./N./(N-2); 1/N(end); zeros(m,1)];
v2=-v0(1:end-1)-v0(end:-1:2);
wcc=ifft(v2+g); weights=[wcc;wcc(1,1)];
k=0:n; nodes=(cos(k*pi/n))'; 

case 4 % GAUSS LEGENDRE
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';

case 5 % TRAPEZOIDAL RULE.
h=2/n;
x=-1:h:1;
w=h*[0.5 ones(1,length(x)-2) 0.5];
nodes=x';
weights=w';
end

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