
function [I,aerr,relerr]=cub_ellblend_adaptive(f,A,B,C,alpha,beta,tol,...
    method)

%--------------------------------------------------------------------------
% Object:
%--------------------------------------------------------------------------
% This routine is a basic adaptive code, to compute quadrature over
% bivariate regions that are linear blending of  two trigonometric arcs 
% with parametric equations
%           P(theta)=A1*cos(theta)+B1*sin(theta)+C1
%           Q(theta)=A2*cos(theta)+B2*sin(theta)+C2
% namely
% R = {(x,y)=t*P(theta)+(1-t)*Q(theta), t in [0,1], theta in [alpha,beta],
% 0<beta-alpha<=2*pi.
%
% As it is well-known, in specific and unfortunate situations, it may not
% provide correct results.
%--------------------------------------------------------------------------
% Input:
%--------------------------------------------------------------------------
% f        : function to integrate over the spherical triangle defined by
%            vertices;
%
% A,B,C: 2x2 matrices of the parametric arc coefficients:
%                 A1=A(1,:), B1=B(1,:), C1=C(1,:)
%                 A2=A(2,:), B2=B(2,:), C2=C(2,:)
%
% [alpha,beta]: angular interval, 0<beta-alpha<=2*pi
%
% tol      : vector of two components establishing the cubature error 
%            tolerance
%            tol(1): 'AbsTol', absolute error tolerance
%            tol(2): 'RelTol', relative error tolerance
%
%            if one uses "method=1" or "method=2", it attempts to satisfy
%                    |IH - I| <= max(AbsTol,RelTol*|IH|),
%            where "IH" is the numerical value computed by the algorithm;
%
%            if "tol" is scalar then "AbsTol=tol", "RelTol=tol".
%
% method   : 1: basic adaptive method
%            2: integral2 (Matlab built-in)
%            3: chebfun2
% 
% IMPORTANT: in case the user chooses "method=3", Chebfun toolbox must be
%            installed.
%--------------------------------------------------------------------------
% Output:
%--------------------------------------------------------------------------
% I        : approximation of the integral;
% aerr     : absolute error estimate ("method=1" only);
% rerr     : relative error estimate ("method=1" only).
%--------------------------------------------------------------------------
% Subroutines:
%--------------------------------------------------------------------------
% 1. cub_polygon_adaptive,
% 2. compute_jacobian,
%
% and their subroutines.
%--------------------------------------------------------------------------
% Dates:
%--------------------------------------------------------------------------
% Initial version: January 7, 2026, by A. Sommariva.
% Last update    : January 7, 2026, by A. Sommariva.
%--------------------------------------------------------------------------


if nargin < 7, tol=10^(-6); end
if nargin < 8, method=1; end

I=[];
aerr=[];

[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;

switch method
    case 1
        %fprintf('\n \t -> Using cub_polygon_adaptive \n');
        X=[a b b a]; Y=[c c d d]; vertices=[X' Y'];
        [I,IL]=cub_polygon_adaptive(vertices,F,tol);
        aerr=abs(I-IL); relerr=aerr/(abs(I)+(I == 0));

    case 2
        %fprintf('\n \t -> Using integral2 \n');
        I = integral2(F,a,b,alpha,beta,'AbsTol',tol(1),'RelTol',tol(end));
        aerr=NaN; relerr=NaN;

    case 3
       %fprintf('\n \t -> Using chebfun2 \n');
        Fch=chebfun2(F,[a b alpha beta]);
        I=sum2(Fch);
        aerr=NaN; relerr=NaN;
end










%--------------------------------------------------------------------------
% Subroutines
%--------------------------------------------------------------------------

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










function [IH,IL,flag,iters,PG,tri_type]=cub_polygon_adaptive(vertices,...
    f,tol,tri_type)

%--------------------------------------------------------------------------
% Object:
%--------------------------------------------------------------------------
% This routine is a basic adaptive code, to compute quadrature over
% polygonal regions of continuous bivariate functions.
%
% As it is well-known, in specific and unfortunate situations, it may not
% provide correct results.
%--------------------------------------------------------------------------
% Input:
%--------------------------------------------------------------------------
% vertices : 1. it is a M x 3 matrix, where the k-th row represents the
%            cartesian coordinates of the k-th vertex of the polygonal
%            region (counterclockwise order);
%            2. alternatively it is a polyshape object defining the region.
%
% f        : function to integrate over the spherical triangle defined by
%            vertices;
%
% tol      : vector of two components establishing the cubature error
%            tolerance
%            tol(1): 'AbsTol', absolute error tolerance
%            tol(2): 'RelTol', relative error tolerance
%
%            it attempts to satisfy |IH - I| <= max(AbsTol,RelTol*|IH|),
%            where "IH" is the numerical value computed by the algorithm;
%
%            if "tol" is scalar then "AbsTol=tol", "RelTol=tol".
%
% tri_type : in case a triangle must be subdivided,
%            0: triangulation via midpoint of each side  (default)
%            1: triangulation via midpoint of longer side
%            2: triangulation via centroid
%--------------------------------------------------------------------------
% Output:
%--------------------------------------------------------------------------
% IH       : high order approximation of the integral;
% IL       : low order approximation of the integral;
% flag     : 0: right termination, 1: some issues;
% iters    : iterations of the routine.
% PG       : domain in polyshape form.
% tri_type : in case a triangle must be subdivided,
%            0: triangulation via midpoint of each side  (default)
%            1: triangulation via midpoint of longer side
%            2: triangulation via centroid
%--------------------------------------------------------------------------
% Subroutines:
%--------------------------------------------------------------------------
% 1. add_triangle_data
% 2. generate_triangle_sons
% 3. cub_tri
% 4. lyness_jespersen
%--------------------------------------------------------------------------
% Dates:
%--------------------------------------------------------------------------
% Initial version: January 05, 2021, by A. Sommariva.
% Relevant update: July 14, 2021, by A. Sommariva.
% Last update: January 2, 2026, by A. Sommariva.
%--------------------------------------------------------------------------

% ................. Troubleshooting and settings ..........................

if nargin < 3, tol=[10^(-6) 10^(-6)]; end
if isscalar(tol), tol=[tol tol]; end
if nargin < 4, tri_type=0; end

AbsTol=tol(1);
RelTol=tol(2);

% max number of triangles in which the domain is partitioned
max_triangles=5000;
safety_factor=0.1;

% ........ Procedure start up, data stored in structures ............

%..........................................................................
% Main strategy:
% 1. Triangulate the spherical polygon, providing the vertices of the
% triangulation and its connectivity list (description of each spherical
% triangle via its vertices).
% 2. Compute low and high order integrals on each triangle.
% 3. Subdivide the triangle with worst error in triangles and compute
%    low and high order integrals in each of them.
% 4. Erase from the list the triangle with worst error.
% 5. Update the list with the new triangles.
% -> Iterate until a criterion is satisfied (absolute and relative
%    error estimates are below a threshold).
%..........................................................................

flag=0;

L1_vertices={}; L1_integrals_L=[]; L1_integrals_H=[]; L1_errors=[];

% If input is not a polyshape determines its polyshape "PG".
S=class(vertices);
if isa(S,'polyshape') == 0
    XV=vertices(:,1); YV=vertices(:,2); PG=polyshape(XV,YV);
else
    PG=vertices;
end

tri=triangulation(PG);
tri_vertices=tri.Points; tri_conn_list=tri.ConnectivityList;

for k=1:size(tri_conn_list,1)
    verticesL=tri_vertices((tri_conn_list(k,:))',:);
    [L1_vertices,L1_integrals_L,L1_integrals_H,L1_errors]=...
        add_triangle_data(verticesL,f,L1_vertices,L1_integrals_L,...
        L1_integrals_H,L1_errors);
end

IL=sum(L1_integrals_L); IH=sum(L1_integrals_H); AE=abs(IL-IH);

% ........................ successive refinement ..........................

iters=1;

while AE > safety_factor*max(AbsTol,RelTol*abs(IH))

    N=length(L1_integrals_L);

    % Too many triangles: exit with errors
    if N > max_triangles, flag=1; return; end

    % Determine triangle with worst error.
    [~,kmax]=max(L1_errors);
    vertices_kmax=L1_vertices{kmax};

    % Erase "kmax" triangle from LIST 1
    k_ok=setdiff(1:N,kmax);

    if not(isempty(k_ok))
        L1_vertices={L1_vertices{k_ok}};
        L1_integrals_L=L1_integrals_L(k_ok);
        L1_integrals_H=L1_integrals_H(k_ok);
        L1_errors=L1_errors(k_ok);
    else
        L1_vertices={}; L1_integrals_L=[]; L1_integrals_H=[]; L1_errors=[];
    end

    new_triangles_cell=generate_triangle_sons(vertices_kmax);

    % Replace worst triangle "data", with "those" of a triangulation.
    for j=1:length(new_triangles_cell)
        verticesL=new_triangles_cell{j};
        [L1_vertices,L1_integrals_L,L1_integrals_H,L1_errors]=...
            add_triangle_data(verticesL,f,L1_vertices,L1_integrals_L,...
            L1_integrals_H,L1_errors);
    end

    IL=sum(L1_integrals_L); IH=sum(L1_integrals_H);
    AE=abs(IL-IH);

    % fprintf('\n \t #: %5.0f AE:%1.1e',iters,AE);
    iters=iters+1;

end









function [L1_vertices,L1_integrals_L,L1_integrals_H,L1_errors]=...
    add_triangle_data(vertices,f,L1_vertices,L1_integrals_L,...
    L1_integrals_H,L1_errors)

%--------------------------------------------------------------------------
% Input:
% vertices : it is a 3 x 3 matrix, where the k-th row represents the
%            cartesian coordinates of the k-th vertex (counterclockwise);
% f        : function to integrate over the spherical triangle defined by
%            vertices;
% L1_vertices: list of vertices determining the triangulation of the
%            polygonal region (it is a cell);
% L1_integrals_L: list of integrals (low order) on the triangles defined by
%            L1_vertices (it is a vector);
% L1_integrals_H: list of integrals (high order) on the triangles defined
%            by L1_vertices (it is a vector);
% L1_errors: list of approximation of absolute errors of the list of
%            integrals on the triangles defined by L1_vertices (it is a
%            vector).
%--------------------------------------------------------------------------
% Output:
% L1_vertices: list of vertices determining the triangulation of the
%            polygonal region (it is a cell);
% L1_integrals_L: list of integrals (low order) on the triangles defined by
%            L1_vertices (it is a vector);
% L1_integrals_H: list of integrals (high order) on the triangles defined
%            by L1_vertices (it is a vector);
% L1_errors: list of approximation of absolute errors of the list of
%            integrals on the triangles defined by L1_vertices (it is a
%            vector).
%--------------------------------------------------------------------------

L1_vertices{end+1}=vertices;
[IL,IH]=cub_tri(f,vertices);
L1_integrals_L(end+1)=IL;
L1_integrals_H(end+1)=IH;
L1_errors(end+1)=abs(IL-IH);









function triangles_cell=generate_triangle_sons(vertices,tri_type)

%--------------------------------------------------------------------------
% Input:
% vertices : it is a 3 x 3 matrix, where the k-th row represents the
%            cartesian coordinates of the k-th vertex (counterclockwise);
% tri_type : 0: triangulation via midpoint of each side (default)
%            1: triangulation via midpoint of longer side
%            2: triangulation via centroid
%--------------------------------------------------------------------------
% Output:
% triangles_cell: coordinates of 4 triangles that determine the "natural"
%           triangulation, via midpoint of each side, of the triangle
%           determined by "vertices".
%--------------------------------------------------------------------------

if nargin < 2, tri_type=0; end

OA=(vertices(1,:)); OB=(vertices(2,:)); OC=(vertices(3,:));


switch tri_type


    case 0 % via sides midpoints
        % .... compute side midpoints ....
        OAB_mid=(OA+OB)/2; OAC_mid=(OA+OC)/2;  OBC_mid=(OB+OC)/2;

        % .... triangle data ....
        triangles_cell{1}=[OA;OAB_mid;OAC_mid];
        triangles_cell{2}=[OAB_mid; OB; OBC_mid];
        triangles_cell{3}=[OBC_mid; OC; OAC_mid];
        triangles_cell{4}=[OAB_mid; OBC_mid; OAC_mid];



    case 1 % via midpoint of longer side

        % ... compute midpoint of longer side ...

        Lsides(1)=norm(OB-OA); % AB
        Lsides(2)=norm(OC-OA); % AC
        Lsides(3)=norm(OC-OB); % BC

        [~,imax]=max(Lsides);

        % .... triangle data ....

        switch imax
            case 1
                OAB_mid=(OA+OB)/2;
                triangles_cell{1}=[OC; OA; OAB_mid];
                triangles_cell{2}=[OB; OC; OAB_mid];
            case 2
                OAC_mid=(OA+OC)/2;
                triangles_cell{1}=[OA; OB; OAC_mid];
                triangles_cell{2}=[OB; OC; OAC_mid];
            case 3
                OCB_mid=(OC+OB)/2;
                triangles_cell{1}=[OA; OB; OCB_mid];
                triangles_cell{2}=[OC; OA; OCB_mid];
        end



    case 2 % via centroid

        % ... compute centroid  ...

        OCC=(OA+OB+OC)/3;

        % .... triangle data ....

        triangles_cell{1}=[OA; OB; OCC];
        triangles_cell{2}=[OB; OC; OCC];
        triangles_cell{3}=[OC; OA; OCC];


end









function [IL,IH]=cub_tri(f,vertices)

%--------------------------------------------------------------------------
% Object:
% Computation of low order and higher order rules of integral of "f" on the
% triangle T determined by vertices.
%--------------------------------------------------------------------------
% Input:
% f: function to integrate;
% vertices: 3 x 3 matrix whose k-th row represents the cartesian
% coordinates of the k-th vertex.
%--------------------------------------------------------------------------
% Output:
% IL: low order approximation of the integral of "f" on the triangle T;
% IH: high order approximation of the integral of "f" on the triangle T.
%--------------------------------------------------------------------------

tri_area=polyarea(vertices(:,1),vertices(:,2));

% orders=[5 9];
% orders=[9 11];
orders=[11 18];

% ........................... low order rule ..............................

xyw_bar=lyness_jespersen(orders(1));

% nodes and weights
lambda=xyw_bar(:,1:3);
nodes=lambda*vertices;

% Low integral
w=tri_area*xyw_bar(:,4); fXY=feval(f,nodes(:,1),nodes(:,2)); IL=w'*fXY;

% ........................... high order rule ..............................

xyw_bar=lyness_jespersen(orders(2));

% nodes and weights
lambda=xyw_bar(:,1:3);
nodes=lambda*vertices;

% Low integral
w=tri_area*xyw_bar(:,4); fXY=feval(f,nodes(:,1),nodes(:,2)); IH=w'*fXY;









function xyw_bar=lyness_jespersen(deg)

%--------------------------------------------------------------------------
% Object:
% Cubature rules on triangles, in baricentric coordinates (sum of the
% weights equal to 1).
%--------------------------------------------------------------------------
% Input:
% deg: algebraic degree of exactness.
%--------------------------------------------------------------------------
% Output:
% xyw_bar: cubature rule in baricentric coordinates where the nodes are
%       (x,y,z) with x=xyw_bar(:,1), y=xyw_bar(:,2), z=xyw_bar(:,3) and the
%       weights are w=xyw_bar(:,4).
%--------------------------------------------------------------------------


switch deg
    case   5

        % ALG. DEG.:    5
        % PTS CARD.:    7
        % NEG. W.  :    0
        % OUT PTS. :    0
        % M.E.INF. : 1.6e-16

        xyw_bar=[
            3.333333333333333148296162562474e-01 3.333333333333333148296162562474e-01 3.333333333333333703407674875052e-01 2.250000000000000055511151231258e-01
            7.974269853530873142233303951798e-01 1.012865073234563428883348024101e-01 1.012865073234563428883348024101e-01 1.259391805448271672851490166067e-01
            1.012865073234563428883348024101e-01 7.974269853530873142233303951798e-01 1.012865073234563428883348024101e-01 1.259391805448271672851490166067e-01
            1.012865073234563428883348024101e-01 1.012865073234563428883348024101e-01 7.974269853530873142233303951798e-01 1.259391805448271672851490166067e-01
            5.971587178976978105282569231349e-02 4.701420641051151094735871538433e-01 4.701420641051150539624359225854e-01 1.323941527885061641978126090180e-01
            4.701420641051151094735871538433e-01 5.971587178976978105282569231349e-02 4.701420641051150539624359225854e-01 1.323941527885061641978126090180e-01
            4.701420641051151094735871538433e-01 4.701420641051151094735871538433e-01 5.971587178976978105282569231349e-02 1.323941527885061641978126090180e-01
            ];



    case   9

        % ALG. DEG.:    9
        % PTS CARD.:   19
        % NEG. W.  :    0
        % OUT PTS. :    0
        % M.E.INF. : 3.1e-16

        xyw_bar=[
            3.33333333333335257719909350271337e-01 3.33333333333331482961625624739099e-01 3.33333333333333259318465024989564e-01 9.71357962827988224985276133338630e-02
            4.89682519198736176946340492577292e-01 4.89682519198739063526204517984297e-01 2.06349616025247595274549894384108e-02 3.13347002271390229211078803928103e-02
            4.89682519198737009613608961444697e-01 2.06349616025245374828500644071028e-02 4.89682519198738397392389742890373e-01 3.13347002271387731209273397325887e-02
            2.06349616025248878969922117221358e-02 4.89682519198738230858936049116892e-01 4.89682519198736898591306498929043e-01 3.13347002271393143546518444964022e-02
            4.37089591492937745709213004374760e-01 4.37089591492935525263163754061679e-01 1.25820817014126673516472010305733e-01 7.78275410047742782770896496913338e-02
            4.37089591492937745709213004374760e-01 1.25820817014125674315749847664847e-01 4.37089591492936524463885916702566e-01 7.78275410047739035768188387010014e-02
            1.25820817014127700472769788575533e-01 4.37089591492935580774314985319506e-01 4.37089591492936690997339610476047e-01 7.78275410047746946107238841250364e-02
            1.88203535619033524017851277676527e-01 1.88203535619031969705616802457371e-01 6.23592928761934617298834382381756e-01 7.96477389272102626049942841746088e-02
            1.88203535619033940351485512110230e-01 6.23592928761932951964297444646945e-01 1.88203535619033135439792658871738e-01 7.96477389272104707718114013914601e-02
            6.23592928761936171611068857600912e-01 1.88203535619031553371982568023668e-01 1.88203535619032275016948574375419e-01 7.96477389272100266826015513288439e-02
            4.47295133944526565605848134055123e-02 4.47295133944527745217811798283947e-02 9.10540973211094617489891334116692e-01 2.55776756586980312524470804191878e-02
            4.47295133944526010494335821476852e-02 9.10540973211094728512193796632346e-01 4.47295133944526357439031016838271e-02 2.55776756586979271690385218107622e-02
            9.10540973211094395445286409085384e-01 4.47295133944528369718263149934501e-02 4.47295133944527675828872759211663e-02 2.55776756586981388053025909812277e-02
            3.68384120547366744613526634566369e-02 7.41198598784498008384957756788936e-01 2.21962989160765289398113964125514e-01 4.32835393772897300546098620088742e-02
            3.68384120547364316000660267036437e-02 2.21962989160765927776353123590525e-01 7.41198598784497564295747906726319e-01 4.32835393772896051545195916787634e-02
            7.41198598784497897362655294273281e-01 3.68384120547358973052354258470587e-02 2.21962989160766205332109279879660e-01 4.32835393772890292263255673788080e-02
            7.41198598784497564295747906726319e-01 2.21962989160766371865562973653141e-01 3.68384120547360638386891196205397e-02 4.32835393772892443320365885028878e-02
            2.21962989160765400420416426641168e-01 3.68384120547361401665220626000519e-02 7.41198598784498452474167606851552e-01 4.32835393772891610653097416161472e-02
            2.21962989160765011842357807836379e-01 7.41198598784498452474167606851552e-01 3.68384120547364801723233540542424e-02 4.32835393772895010711110330703377e-02
            ];

    case  11

        % ALG. DEG.:   11
        % PTS CARD.:   28
        % NEG. W.  :    0
        % OUT PTS. :    0
        % M.E.INF. : 4.2e-16

        xyw_bar=[
            3.33333333333333425851918718763045e-01 3.33333333333333370340767487505218e-01 3.33333333333333148296162562473910e-01 8.79773011622321937652557721776247e-02
            2.59891409282874893960091355893383e-02 2.59891409282874824571152316821099e-02 9.48021718143424951819042689749040e-01 8.74431155373607202352381762011646e-03
            2.59891409282874408237518082387396e-02 9.48021718143425062841345152264694e-01 2.59891409282875240904786551254801e-02 8.74431155373604947211862992162423e-03
            9.48021718143425173863647614780348e-01 2.59891409282874408237518082387396e-02 2.59891409282873853126005769809126e-02 8.74431155373604773739515394481714e-03
            9.42875026479226829856372660287889e-02 9.42875026479227384967884972866159e-02 8.11424994704154634028725467942422e-01 3.80815719939349636713465940829337e-02
            9.42875026479225580855469956986781e-02 8.11424994704154745051027930458076e-01 9.42875026479227384967884972866159e-02 3.80815719939349636713465940829337e-02
            8.11424994704154634028725467942422e-01 9.42875026479225997189104191420483e-02 9.42875026479227662523641129155294e-02 3.80815719939349428546648823612486e-02
            4.94636775017213814464867027709261e-01 4.94636775017213869976018258967088e-01 1.07264499655723155591147133236518e-02 1.88554480561312423625430767515354e-02
            4.94636775017213758953715796451434e-01 1.07264499655723415799668529757582e-02 4.94636775017213869976018258967088e-01 1.88554480561312527708839326123780e-02
            1.07264499655723502535842328597937e-02 4.94636775017213814464867027709261e-01 4.94636775017213814464867027709261e-01 1.88554480561312597097778365196064e-02
            2.07343382614511434480775164956867e-01 2.07343382614511434480775164956867e-01 5.85313234770977075527298438828439e-01 7.21596975447395400093952844144951e-02
            2.07343382614511378969623933699040e-01 5.85313234770977297571903363859747e-01 2.07343382614511351214048318070127e-01 7.21596975447395122538196687855816e-02
            5.85313234770977297571903363859747e-01 2.07343382614511351214048318070127e-01 2.07343382614511351214048318070127e-01 7.21596975447395122538196687855816e-02
            4.38907805700492092970677049379447e-01 4.38907805700492148481828280637274e-01 1.22184388599015814058645901241107e-01 6.93291387055358782065539458017156e-02
            4.38907805700492203992979511895101e-01 1.22184388599015814058645901241107e-01 4.38907805700491926437223355605965e-01 6.93291387055358782065539458017156e-02
            1.22184388599015911203160555942304e-01 4.38907805700492148481828280637274e-01 4.38907805700491926437223355605965e-01 6.93291387055359059621295614306291e-02
            7.85647334377654633569776739674958e-17 8.58870281282636427455656757956604e-01 1.41129718717363461522040779527742e-01 7.36238378330058337167818294233257e-03
            1.57669518850046390999681416911743e-16 1.41129718717363572544343242043396e-01 8.58870281282636316433354295440949e-01 7.36238378330061546406248851326382e-03
            8.58870281282636427455656757956604e-01 1.42062476402819843287610091505200e-16 1.41129718717363433766465163898829e-01 7.36238378330060505572163265242125e-03
            8.58870281282636427455656757956604e-01 1.41129718717363350499738317012088e-01 2.22044604925031308084726333618164e-16 7.36238378330062066823291644368510e-03
            1.41129718717363572544343242043396e-01 2.02916231519083255217512748366332e-16 8.58870281282636205411051832925295e-01 7.36238378330063020921203431612412e-03
            1.41129718717363378255313932641002e-01 8.58870281282636427455656757956604e-01 2.22044604925031308084726333618164e-16 7.36238378330061199461553655964963e-03
            4.48416775891304561496575331602799e-02 6.77937654882590168270439789921511e-01 2.77220667528279340885433157382067e-01 4.10563154292885590379569293872919e-02
            4.48416775891304908441270526964217e-02 2.77220667528279451907735619897721e-01 6.77937654882590057248137327405857e-01 4.10563154292885729157447372017486e-02
            6.77937654882590168270439789921511e-01 4.48416775891304561496575331602799e-02 2.77220667528279396396584388639894e-01 4.10563154292885590379569293872919e-02
            6.77937654882590168270439789921511e-01 2.77220667528279340885433157382067e-01 4.48416775891304908441270526964217e-02 4.10563154292885451601691215728351e-02
            2.77220667528279340885433157382067e-01 4.48416775891305255385965722325636e-02 6.77937654882590168270439789921511e-01 4.10563154292885798546386411089770e-02
            2.77220667528279285374281926124240e-01 6.77937654882590279292742252437165e-01 4.48416775891303798218245901807677e-02 4.10563154292885590379569293872919e-02
            ];

    case 18

        % ALG. DEG.:   18
        % PTS CARD.:   66
        % NEG. W.  :    0
        % OUT PTS. :    0
        % M.E.INF. : 6.1e-16

        xyw_bar=[3.697809640550978396711379048156e-03 7.435792350812416018346340251810e-02 9.219442668513249117268060217611e-01 2.637791126882527061731087414387e-03
            7.435792350812416018346340251810e-02 9.219442668513248007045035592455e-01 3.697809640550997478669614793034e-03 2.637791126882527061731087414387e-03
            9.219442668513248007045035592455e-01 3.697809640550978396711379048156e-03 7.435792350812420181682682596147e-02 2.637791126882527061731087414387e-03
            1.491439922799888545545954343652e-02 1.369091307336264828897753886849e-02 9.713946876986384593166690137878e-01 2.674906650996751195786815813449e-03
            1.369091307336264828897753886849e-02 9.713946876986384593166690137878e-01 1.491439922799886463877783171483e-02 2.674906650996751195786815813449e-03
            9.713946876986384593166690137878e-01 1.491439922799888545545954343652e-02 1.369091307336267604455315449741e-02 2.674906650996751195786815813449e-03
            1.697399962193029071891303694031e-03 2.089136906154610839880803041524e-01 7.893889094223458524623993071145e-01 3.193038273763521006642918109719e-03
            2.089136906154610839880803041524e-01 7.893889094223459634847017696302e-01 1.697399962192980282793541846331e-03 3.193038273763521006642918109719e-03
            7.893889094223459634847017696302e-01 1.697399962193029071891303694031e-03 2.089136906154610562325046885235e-01 3.193038273763521006642918109719e-03
            8.066208394662527370044813324057e-02 8.243532690291427780793220847499e-03 9.110943833630833488257394492393e-01 4.170698437610267927944107668736e-03
            8.243532690291427780793220847499e-03 9.110943833630833488257394492393e-01 8.066208394662521818929690198274e-02 4.170698437610267927944107668736e-03
            9.110943833630833488257394492393e-01 8.066208394662527370044813324057e-02 8.243532690291432984963648777921e-03 4.170698437610267927944107668736e-03
            1.987479014023108003783590902458e-01 6.976595096183165964354166987960e-03 7.942755035015059972280937472533e-01 5.259203237678910895780415302170e-03
            6.976595096183165964354166987960e-03 7.942755035015059972280937472533e-01 1.987479014023108003783590902458e-01 5.259203237678910895780415302170e-03
            7.942755035015059972280937472533e-01 1.987479014023108003783590902458e-01 6.976595096183202393547162500909e-03 5.259203237678910895780415302170e-03
            3.467027191472764058310929158324e-01 9.059771353479850741763357291347e-03 6.442375094992437434271437268762e-01 7.188445007987385830527937713441e-03
            9.059771353479850741763357291347e-03 6.442375094992436324048412643606e-01 3.467027191472764613422441470902e-01 7.188445007987385830527937713441e-03
            6.442375094992436324048412643606e-01 3.467027191472764058310929158324e-01 9.059771353479906252914588549174e-03 7.188445007987385830527937713441e-03
            4.863463338710427574795858163270e-02 5.579105297174669769733057478334e-02 8.955743136411490334936047474912e-01 8.958782692065399982928042277308e-03
            5.579105297174669769733057478334e-02 8.955743136411490334936047474912e-01 4.863463338710427574795858163270e-02 8.958782692065399982928042277308e-03
            8.955743136411490334936047474912e-01 4.863463338710427574795858163270e-02 5.579105297174669075843667087611e-02 8.958782692065399982928042277308e-03
            1.283186298494463049857383651897e-02 3.446925484779704285109858119540e-01 6.424755885370849739501863950863e-01 9.248725387125594776804859975528e-03
            3.446925484779704285109858119540e-01 6.424755885370849739501863950863e-01 1.283186298494465304997902421746e-02 9.248725387125594776804859975528e-03
            6.424755885370849739501863950863e-01 1.283186298494463049857383651897e-02 3.446925484779703729998345806962e-01 9.248725387125594776804859975528e-03
            2.935953899422051593948168601855e-02 1.399641103183917323349305661395e-01 8.306763506873877656033755556564e-01 1.007424624849140998295027316090e-02
            1.399641103183917323349305661395e-01 8.306763506873877656033755556564e-01 2.935953899422050206169387820410e-02 1.007424624849140998295027316090e-02
            8.306763506873877656033755556564e-01 2.935953899422051593948168601855e-02 1.399641103183917323349305661395e-01 1.007424624849140998295027316090e-02
            4.907883946937592889447898869548e-01 4.958225177896646296282767707453e-01 1.338908751657608142693334229989e-02 1.023007488644936922861372607940e-02
            4.958225177896646296282767707453e-01 1.338908751657613693808457355772e-02 4.907883946937592334336386556970e-01 1.023007488644936922861372607940e-02
            1.338908751657613693808457355772e-02 4.907883946937592889447898869548e-01 4.958225177896645741171255394875e-01 1.023007488644936922861372607940e-02
            1.386187840597887799276577425189e-01 4.470767938550187364787902311036e-02 8.166735365547093117299937148346e-01 1.273721618736533590932058501721e-02
            4.470767938550187364787902311036e-02 8.166735365547093117299937148346e-01 1.386187840597887799276577425189e-01 1.273721618736533590932058501721e-02
            8.166735365547093117299937148346e-01 1.386187840597887799276577425189e-01 4.470767938550190834234854264650e-02 1.273721618736533590932058501721e-02
            2.638976063522009973816295769211e-01 4.367209525175506923400092773591e-02 6.924302983960439750177329187864e-01 1.412698155144117943826653771566e-02
            4.367209525175506923400092773591e-02 6.924302983960439750177329187864e-01 2.638976063522009418704783456633e-01 1.412698155144117943826653771566e-02
            6.924302983960439750177329187864e-01 2.638976063522009973816295769211e-01 4.367209525175508311178873555036e-02 1.412698155144117943826653771566e-02
            4.559117449306109015605414924721e-02 2.452185735061957194691473205239e-01 7.091902520007431487414351067855e-01 1.533277441555348612201026270441e-02
            2.452185735061957194691473205239e-01 7.091902520007432597637375693012e-01 4.559117449306104852269072580384e-02 1.533277441555348612201026270441e-02
            7.091902520007432597637375693012e-01 4.559117449306109015605414924721e-02 2.452185735061956917135717048950e-01 1.533277441555348612201026270441e-02
            3.940103585370680616684069264011e-01 5.927693378164468340774106991375e-02 5.467127076812872132904885802418e-01 1.915098830446838179542012881029e-02
            5.927693378164468340774106991375e-02 5.467127076812873243127910427575e-01 3.940103585370680061572556951432e-01 1.915098830446838179542012881029e-02
            5.467127076812873243127910427575e-01 3.940103585370680616684069264011e-01 5.927693378164455850765079958364e-02 1.915098830446838179542012881029e-02
            1.068672069409647568116739080324e-01 1.283255269479803006849039093140e-01 7.648072661110549841367856060970e-01 1.958137585110737710625272711695e-02
            1.283255269479803006849039093140e-01 7.648072661110549841367856060970e-01 1.068672069409647429338861002179e-01 1.958137585110737710625272711695e-02
            7.648072661110549841367856060970e-01 1.068672069409647568116739080324e-01 1.283255269479802729293282936851e-01 1.958137585110737710625272711695e-02
            2.205752794254834114884289419933e-01 1.132758213016846898790745967744e-01 6.661488992728319402658598846756e-01 2.251921487278832692546792770827e-02
            1.132758213016846898790745967744e-01 6.661488992728319402658598846756e-01 2.205752794254833837328533263644e-01 2.251921487278832692546792770827e-02
            6.661488992728319402658598846756e-01 2.205752794254834114884289419933e-01 1.132758213016846760012867889600e-01 2.251921487278832692546792770827e-02
            1.236538350412987319337787539553e-01 2.427204372160924328305497965630e-01 6.336257277426088352356714494817e-01 2.374408805361872057404504232636e-02
            2.427204372160924328305497965630e-01 6.336257277426087242133689869661e-01 1.236538350412987874449299852131e-01 2.374408805361872057404504232636e-02
            6.336257277426087242133689869661e-01 1.236538350412987319337787539553e-01 2.427204372160924883417010278208e-01 2.374408805361872057404504232636e-02
            7.133487941009795940683346771038e-02 3.879084789424865342688519831427e-01 5.407566416474154369353755100747e-01 2.451258023475103020594190184056e-02
            3.879084789424865342688519831427e-01 5.407566416474154369353755100747e-01 7.133487941009808430692373804050e-02 2.451258023475103020594190184056e-02
            5.407566416474154369353755100747e-01 7.133487941009795940683346771038e-02 3.879084789424865897800032144005e-01 2.451258023475103020594190184056e-02
            2.241269103249716176939898559795e-01 2.132898000096189738972185523380e-01 5.625832896654093806532159760536e-01 2.468147360226879277744771457037e-02
            2.132898000096189738972185523380e-01 5.625832896654093806532159760536e-01 2.241269103249716732051410872373e-01 2.468147360226879277744771457037e-02
            5.625832896654093806532159760536e-01 2.241269103249716176939898559795e-01 2.132898000096190571639453992248e-01 2.468147360226879277744771457037e-02
            3.516190499556426929395058778027e-01 1.379375099764565126836401987021e-01 5.104434400679007666212783078663e-01 2.696766291995268335912605550675e-02
            1.379375099764565126836401987021e-01 5.104434400679008776435807703820e-01 3.516190499556426374283546465449e-01 2.696766291995268335912605550675e-02
            5.104434400679008776435807703820e-01 3.516190499556426929395058778027e-01 1.379375099764563739057621205575e-01 2.696766291995268335912605550675e-02
            3.607623405677324313067799721466e-01 4.651258199647673174936812756641e-01 1.741118394675003067106899834471e-01 3.234621087379146303630861325473e-02
            4.651258199647673174936812756641e-01 1.741118394675001956883875209314e-01 3.607623405677324868179312034044e-01 3.234621087379146303630861325473e-02
            1.741118394675001956883875209314e-01 3.607623405677324313067799721466e-01 4.651258199647674285159837381798e-01 3.234621087379146303630861325473e-02
            3.113001289825194084492920865159e-01 2.700736293886714101120105624432e-01 4.186262416288091259275461197831e-01 3.399685451717549683481323086198e-02
            2.700736293886714101120105624432e-01 4.186262416288091259275461197831e-01 3.113001289825194639604433177738e-01 3.399685451717549683481323086198e-02
            4.186262416288091259275461197831e-01 3.113001289825194084492920865159e-01 2.700736293886714101120105624432e-01 3.399685451717549683481323086198e-02];

end









