function fx=logan_shepp(pts,j,n)

% INPUT:
% pts: IT IS A MATRIX N x 2, SUCH THAT pts_k=(x_k,y_k).
% j: IT IS A COLUMN VECTOR M x 1 WHOSE COMPONENTS ARE INTEGERS,
%    LESS OR EQUAL TO n (SEE THIRD INPUT OF logan_shepp).
% n: LOGAN-SHEPP DEGREE.

% OUTPUT:
% fx: EVALUATION OF j-th LOGAN-SHEPP POLYNOMIAL OF ORDER n
%     AT pts. IT IS A MATRIX, WITH N ROWS AND M COLUMNS.

normalized_ls=0; % 0: normalized, 1: not normalized.

x=pts(:,1); y=pts(:,2);

M=size(j,1);
N=size(pts,1);

x=repmat(x,M,1); y=repmat(y,M,1);
j=repmat(j,1,N); j=j'; j=j(:);

args= x.*cos( j*pi/(n+1) ) + y.*sin( j*pi/(n+1) );

if normalized_ls == 0
    fx=(1/sqrt(pi))*chebyshev_second_type(args,n);
else
    fx=chebyshev_second_type(args,n);
end

fx=reshape(fx,N,M);

function fx=chebyshev_second_type(x,n)

method=2;

switch method 
case 1
    theta=acos(x);
    s=sin(theta);
    s=s+(1-abs(sign(s)))*eps;
    fx=sin((n+1)*theta)./s;
case 2
    
    % WE COMPUTE RECURSIVELY U_n=sin((n+1)*theta)./s.
    un_2=ones(size(x));
    un_1=2*x;
    
    switch n
    case 0
        fx=un_2;
    case 1
        fx=un_1;
    otherwise
        
        for index=2:n
            un=2*x.*un_1-un_2;
            un_2=un_1;
            un_1=un;
        end
        fx=un;
    end
    
end
