function V = SZEFw(deg,wam,pts,method)
% -------------------------------------------------------------------------
% computes a discrete version of Siciak Zaharjuta Extremal Function
% by a Weakly Admissible Mesh in R^2
% by Federico Piazzon and Marco Vianello, October 2015
% -------------------------------------------------------------------------
% INPUT
% deg: polynomial degree
% wam: 2-column array of mesh points for degree deg
% pts: 2-column array of target points
% method: string, 'kernel','weighhtedkernel', 'bergman', 'weightedbergman'
% -------------------------------------------------------------------------
% OUTPUT
% V: 1-column array of values of the extremal function at pts
% -------------------------------------------------------------------------

% FUNCTION BODY

% rectangle containing the mesh wam
rect=[min(wam(:,1)) max(wam(:,1)) min(wam(:,2)) max(wam(:,2))];
M=size(wam,1);
N=(deg+1)*(deg+2)/2;
% Chebyshev-Vandermonde matrix at the mesh points
V=chebvandr(deg,wam,rect);
% 2-step orthogonalization with uniform weights
w0=sqrt(M).^(-1).*ones(1,M);
W0=sparse(diag(w0));
[Q0,R1]=qr(W0*V,0);
[Q1,R2]=qr(Q0,0);
Q1=sparse(diag(w0.^(-1)))*Q1;
ort=Q1'*W0.^2*Q1;
fprintf('orthogonality check %1.16f\n',norm(ort-eye(size(ort,1))))
switch method
    case 'kernel'
        Vop =(chebvandr(deg,pts,rect)/R1/R2);
        V = 1/deg .* log(sum(abs(Q1*Vop'))./M);
    case 'bergman'
        Vop =(chebvandr(deg,pts,rect)/R1/R2);
        V = 1/(2*deg).* log((sum(abs(Vop).^2,2)));
    case 'weightedkernel'
        % orthonormalization w.r.t Bergman weight
        w=(sum(abs(Q1).^2,2)./N);
        W=sparse(diag(w.^(1/2)));
        [Q,R3] = qr(W*Q1,0);
        DOP=(chebvandr(deg,pts,rect)/R1/R2/R3);
        V=1/(deg).*log(sum(abs(W*Q*DOP')));
        
    case 'weightedbergman'
          w=(sum(abs(Q1).^2,2)./(N*M));
          W=sparse(diag(w.^(1/2)));
          [Q,R3] = qr(W*Q1,0);
          %[Q, R4] =qr(Q,0);%opzionale?
          Q=diag(w.^(-1))*Q;
          DOP=(((chebvandr(deg,pts,rect)/R1)/R2)/R3);
          %DOP=(((chebvandr(deg,pts,rect)/R1)/R2)/R3)/R4;
          
          V= 1/(2*deg).* log((sum(abs(DOP).^2,2))'); 
          
      
          
end
end

%----------- SUBROUTINES --------------------------------------------------

function V = chebvandr(deg,gmesh,rect)
%--------------------------------------------------------------------------
% computes the bivariate Chebyshev-Vandermonde matrix (with the Chebyshev
% basis defined on the rectangle rect) at the target points gmesh
% INPUT:
% deg = polynomial degree
% gmesh = 2-column array of mesh point coordinates
% rect = 4-component vector such that the rectangle
% [rect(1),rect(2)] x [rect(3),rect(4)] contains the mesh
% -------------------------------------------------------------------------
% OUTPUT:
% V = Chebyshev-Vandermonde matrix at gmesh
%--------------------------------------------------------------------------
% FUNCTION BODY
% rectangle containing the mesh
if isempty(rect)
    rect=[min(gmesh(:,1)) max(gmesh(:,1)) min(gmesh(:,2)) max(gmesh(:,2))];
end;
% couples with length less or equal to deg
% graded lexicographical order
j=(0:1:deg);
[j1,j2]=meshgrid(j);
dim=(deg+1)*(deg+2)/2;
couples=zeros(dim,2);
for s=0:deg
    good=find(j1(:)+j2(:)==s);
    couples(1+s*(s+1)/2:(s+1)*(s+2)/2,:)=[j1(good) j2(good)];
end;
% mapping the mesh in the square [-1,1]^2
a=rect(1);b=rect(2);c=rect(3);d=rect(4);
map=[(2*gmesh(:,1)-b-a)/(b-a) (2*gmesh(:,2)-d-c)/(d-c)];
% Chebyshev-Vandermonde matrix on the mesh
for j=1:length(couples(:,1))
    V(:,j)=chebpoly(couples(j,1),map(:,1)).*chebpoly(couples(j,2),map(:,2));
end
end


function tn = chebpoly(deg,x)
% -------------------------------------------------------------------------
% computes the Chebyshev polynomials on the real line by recurrence
t0=ones(length(x),1);
t1=x;
if deg==0
    tn=t0;
end;
if deg>0
    for j=2:deg
        t2=2*x.*t1-t0;
        t0=t1;
        t1=t2;
    end;
    tn=t1;
end;
end

