function P = rectprojmat( M, N ) 
% Calculate the MxN rectangular projection matrix for 
% spectral colloction.
% 
% Alex Townsend 

% chebpts of 1st kind:
y = sin(pi*((-M+1:2:M-1)/(2*M))).';
% chebpts of 2nd kind:
m = N - 1; x = sin(pi*(-m:2:m)/(2*m)).';
    
% projection matrix can be formed using barycentric formula: 
P = bsxfun(@minus, y, x.');

% Barycentric weights: 
w = ones(1, N);
w(2:2:end) = -1;
w([1, N]) = 0.5*w([1, N]);

% w(k)/(y(j)-x(k))
for k = 1:N
    P(:,k) = w(k)./P(:,k);          
end

% Normalisation.
c = 1./sum(P, 2);                   
for j = 1:M
    P(j,:) = P(j,:)*c(j);
end

end