
function [abs_jac,U]=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
    P=@(theta) A1'*cos(theta)+B1'*sin(theta)+C1';
    Q=@(theta) A2'*cos(theta)+B2'*sin(theta)+C2';
    U=@(t,theta) t*P(theta)+(1-t)*Q(theta);
end
