function f = foursum(x,fh0,fhpos,fhneg)
% f = foursum(x,fh0,fhpos,fhneg)
% evaluate Fourier sum with coefficients given by
%   fh0:    fhat_0
%   fhpos:  vector [fhat_1,...,fhat_n]
%   fhneg:  vector [fhat_{-1},...,fhat_{-n}]
% at the points in vector x
%
% For a real-valued function we have fhneg=conj(fhpos)
% and this argument can be omitted:
% f = fourser(x,fh0,fhpos)   % uses fhneg=conj(fhpos)
%
% Example: fhat_0=0, fhat_k=i/k for nonzero k
%   nv = 1:100;              % vector of n values
%   x = 0:.001:3;            % points for plotting
%   plot( x, foursum(x,0,1i./nv) )

% This is the straightforward way to evaluate a Fourier series.
% We will later use a more efficient way to do this with fft.

if nargin==3
 fhneg = conj(fhpos);
end

f = fh0*ones(size(x));
for k=1:length(fhpos)
  f = f + fhpos(k)*exp(k*2i*pi*x);
end
for k=1:length(fhneg)
  f = f + fhneg(k)*exp(-k*2i*pi*x);
end
