function [x, c] = fejer(N, dom)
%   FEJER(N) returns N Chebyshev points X in (-1,1).
%
%   GAUSS(N, [A, B]) returns N Chebyshev points X in (C,B).
%
%   [X, C] = FEJER(N)  or [X, W] = FEJER(N, [A, B]) returns also a row vector W
%   of weights for Fejer quadrature.

% Nick Hale - Mar 2017

% Theta:
t = (2*(1:N)-1)*pi/(2*N);

% Chebyshev points:
x = -cos(t)';

% Chebyshev moments (RHS)
rhs = 2./(1-(0:N-1).^2)'; rhs(2:2:N) = 0;

% DCT matrix: (note - one can do better by using the FFT)
V = cos(t'*(0:N-1));
c = V'\rhs;

% Scale to arbitrary interval:
if ( nargin == 2)
    a = dom(1);
    b = dom(2);
    x = .5*(b-a)*x + .5*(b+a); 
    c = .5*(b-a)*c;
end

end
