function [aj,bj,cj,kj] = fft314(xn,fn);
%
%  fft314.m -- djm -- 28 jan 2008
%
%  [cj,kj] = fft314(xn,fn) -- inputs:  xn, fn (col vect)
%                          outputs:  cj, kj (col vect)
%
%  !!! DO NOT INCLUDE PERIODIC POINT !!!

%  system size
N = size(xn,1);
if (2*floor(N/2) ~= N)
    disp('ERROR:  # of x-points must be EVEN')
    cj = [];  kj = [];
    return
else
    N2 = round(N/2);
end

%  left endpoint & interval half-width
x0 = xn(1);  Lx = xn(N/2+1)-xn(1);

%  wavenumber vector
kj = (pi/Lx)*(-N2:N2-1)';

%  complex fourier coefficients
cj = exp(i*x0*kj).*fftshift(fft(fn))/N;

%  *** zero the "nyquist" c_j ***
cj(1) = 0;

%  real fourier coefficients
aj =  2*real(cj(N2+1:end));
bj = -2*imag(cj(N2+1:end));
