function [valFS] = evalFS (Fcoeff, N)
%
% evalFS evaluates a trigonometric polynomial in N equidistant points
% for given Fourier coefficients Fcoeff stored in the format 
% [a_M,...,a_1,a_0,b_1,...b_{M-1}]. 
% Regard that N should be even, else it is increased by 1.
%

% Determine length of Fcoeff. Make sure it's a column.
n = length (Fcoeff);
Fcoeff = reshape (Fcoeff, n, 1);

% Determine length of the coefficient vector and of the sine and
% cosine coefficients.
if mod (n,2)
    display ('number of coefficients is not even')
    valFS = 0;
    return;
end

M = n / 2;

% Check size of N and increase it to the minimal possible size.
if N < n
    N = n;
end

% If N was odd on the input, increase it by 1
if mod (N,2)
    N = N + 1;
end

% Compute the bias for the actual length of the Fourier series
% and the desired number of evaluation points.
bias = (N - n) / 2;

% Initialize the complex coefficient vector for the matlab builtin
% inverse Fourier transformation ifft.
beta = zeros (N, 1);

% Set the Fourier coefficents and regard the required zeros
% by the bias. This is the tricky part...
beta (2:M+1) =        .5 * (flipud (Fcoeff(1:M)) - i * [Fcoeff(M+2:n); 0]);
beta (M+1+2*bias:N) = .5 * (Fcoeff(1:M) + i * flipud ([Fcoeff(M+2:n); 0]));

% Set first Fourier coefficient and repair weight of a_M.                             
beta (1) =        Fcoeff(M+1) / 2;
beta (M+1+bias) = beta(M+1+bias) / 2;

% Evaluate polynomial by the ifft and correct the weights.
valFS = N * ifft (beta);

end