
function [I,aerr]=cub_ellblend_adaptive(f,domain_structure,tol)

if nargin < 3, tol=10^(-14); end

I=[];
aerr=[];

dom_str=domain_structure.domain;

A=domain_structure.A;
B=domain_structure.B;
C=domain_structure.C;
alpha=domain_structure.alpha;
beta=domain_structure.beta;

[abs_jac,UX,UY]=compute_jacobian(A,B,C,alpha,beta);

F=@(t,theta) f(UX(t,theta),UY(t,theta)).*abs_jac(t,theta);

a=0; b=1; c=alpha; d=beta; 

method=2;

tic;
switch method
    case 1
        X=[a b b a]; Y=[c c d d]; vertices=[X' Y'];
        [I,IL,flag,iters]=cub_polygon_adaptive(vertices,F,tol);
        aerr=abs(I-IL);
    case 2
        Fch=chebfun2(F,[a b alpha beta]);
        I=sum2(Fch);
        aerr=NaN;
end
toc;






function [abs_jac,UX,UY]=compute_jacobian(A,B,C,alpha,beta)

% "Algebraic cubature by linear blending of elliptical arcs"
% Gaspare Da Fies, Alvise Sommariva and Marco Vianello
% -> see section: "Linear blending of elliptical arcs";

% A: matrix 2 x 2: A1=A(1,:), A2=A(2,:).
% B: matrix 2 x 2: B1=B(1,:), B2=B(2,:).
% C: matrix 2 x 2: C1=C(1,:), C2=C(2,:).

A1=A(1,:); A2=A(2,:);
B1=B(1,:); B2=B(2,:);
C1=C(1,:); C2=C(2,:);

a11=A1(1); a12=A1(2);
a21=A2(1); a22=A2(2);

b11=B1(1); b12=B1(2);
b21=B2(1); b22=B2(2);

c11=C1(1); c12=C1(2);
c21=C2(1); c22=C2(2);

u0=(a11-a21)*(b12-b22)+(a12-a22)*(b21-b11);
u1=(b12-b22)*(c11-c21)+(b21-b11)*(c12-c22);
u2=(a11-a21)*(c12-c22)+(a12-a22)*(c21-c11);

v0=b21*(a22-a12)+b22*(a11-a21);
v1=b21*(c22-c12)+b22*(c11-c21);
v2=a21*(c12-c22)+a22*(c21-c11);
v3=a12*a21-a11*a22+b11*b22-b12*b21;
v4=a12*b21-a11*b22+a21*b12-a22*b11;

u=@(theta) u0+u1*cos(theta)+u2*sin(theta);
v=@(theta) v0+v1*cos(theta)+v2*sin(theta)+v3*cos(theta).*sin(theta)+v4*(sin(theta)).^2;

abm=(alpha+beta)/2;
c=sign(0.5*u(abm)+v(abm));

abs_jac=@(t,theta) c*(t.*u(theta)+v(theta));

if nargout >= 2
    PX=@(theta) a11*cos(theta)+b11*sin(theta)+c11;
    PY=@(theta) a12*cos(theta)+b12*sin(theta)+c12;
    QX=@(theta) a21*cos(theta)+b21*sin(theta)+c21;
    QY=@(theta) a22*cos(theta)+b22*sin(theta)+c22;

    UX=@(t,theta) t.*PX(theta)+(1-t).*QX(theta);
    UY=@(t,theta) t.*PY(theta)+(1-t).*QY(theta);
end




