function [fek,R1,R2,leb] = approxfek_sph(deg,mymesh,leja,controlmesh);


% extracts approximate Fekete or Leja interpolation points from a mesh on
% the sphere and estimates their Lebesgue constant

% by Alvise Sommariva and Marco Vianello, University of Padova 
% November 22, 2012

% see, e.g., the papers: 
%
% A. Sommariva and M. Vianello, Computing approximate Fekete points by QR 
% factorizations of Vandermonde matrices, Comput. Math. Appl. 57 (2009)
% preprint: http://www.math.unipd.it/~marcov/pdf/afek.pdf
% 
% L. Bos, S. De Marchi, A. Sommariva and M. Vianello, Computing 
% multivariate Fekete and Leja points by numerical linear algebra, 
% SIAM J. Numer. Anal. 48 (2010)
% preprint: http://www.math.unipd.it/~marcov/pdf/leja.pdf


% INPUT: 

% deg = interpolation degree

% mymesh = 2-column array of mesh points coordinates, from which the
% approximate Fekete or Leja points will be excracted; if mymesh is 
% already a set of interpolation points, their Lebesgue constant will be 
% estimated on controlmesh
% the coordinates on the sphere must be cylindrical coordinates, that is 
% [phi,z] so that
% x_1=cos(phi)*sqrt(1-z^2)
% x_2=sin(phi)*sqrt(1-z^2)
% x_3=z

% leja=0: computes approximate Fekete points
% leja=1: computes approximate Leja points

% controlmesh = 2-column array of mesh point coordinates
% (possibly finer mesh to evaluate the Lebesgue constant of the points) 
% if controlmesh == [ ], then it will be set controlmesh = mymesh

% OUTPUT: 

% fek = 2-column array of approximate Fekete or Leja points coordinates

% P = transformation matrix for the total-degree product Chebyshev basis 
% of the rectangle, namely 
% [new basis] = [old basis] * P
% [new Vandermonde] = [old Vandermonde] * P

% leb = approximate Lebesgue constant
% a call without this output argument will directly avoid to compute it 


% FUNCTION BODY

% approximate Fekete or Leja points extraction from mymesh 



% computing the rectangular Chebyshev-Vandermonde matrix on the mesh 
W=van_spharm_s2(deg,mymesh(:,2),mymesh(:,1));

dim=0;
for i=0:1:deg
    dim=dim+(2*i+3-2)*factorial(i+3-3)/(factorial(i)*factorial(3-2));
end

% change of basis via 2 successive QR factorizations of the 
% Vandermonde matrices; to keep the original basis set ort=0
ort=1; 
B=W;
if ort==1
[Q,R1]=qr(B,0);
B=B/R1;
[Q,R2]=qr(B,0);
B=B/R2;
else
R1=eye(dim); R2=R1;
end;

% P is the basis transformation matrix: 
% now B = W * P is the Chebyshev-Vandermonde matrix on mymesh 
% in the new basis and is nearly-orthogonal 
% (unless W is too severely ill conditioned) 

% extracting the interpolation points and computing the weights 

if leja == 0
% approximate Fekete points and corresponding weights
% by QR factorization with column pivoting of B'=(W*P)'
w=B'\ones(dim,1); 
ind=find(abs(w)>0);
fek=mymesh(ind,:);
else
% approximate Leja points and corresponding weights  
% by LU factorization with row pivoting of B=W*P
[L,U,perm]=lu(B,'vector');
fek=mymesh(perm(1:dim),:);
end;

% estimating the Lebesgue constant

% square Chebyshev-Vandermonde matrix at the interpolation points
V=van_spharm_s2(deg,mymesh(:,2),mymesh(:,1));

% setting controlmesh if not defined and computing the rectangular 
% Chebyshev-Vandermonde matrix on controlmesh
if nargin==3
controlmesh=mymesh;
controlV=W;
else
controlV=van_spharm_s2(deg,controlmesh(:,2),controlmesh(:,1));
end;

if nargout == 4
% numerical maximization of the Lebesgue function on controlmesh, 
% stabilizing by the basis transformation matrix 
leb=norm(((V/R1)/R2)'\((controlV/R1)/R2)',1);
end;

end