function [a, b] = fourierCoeffs(f, k, r_end = pi)
% function [a, b] = fourierCoeffs(f, k, r_end = pi)
%
%   This routine computes the classical Fourier series coefficients
%   for f up to order k on the interval [-r_end, r_end].
%
%   INPUTS:
%     f      a function handle, pointing to the function to approximate
%     k      the number of sine terms to include
%     r_end  the right endpoint for the spatial domain; defaults to pi
%
%   OUTPUTS:
%      a     the coefficients of the cosine terms
%      b     the coefficients of the sine terms
  r_end = abs(r_end);
  l_end = -r_end;

  a = quad(f, l_end, r_end);
  b = [];
  for m = 1:k
    a = [a; quad(@(x) f(x).*cos(m*pi*x/r_end), l_end, r_end)];
    b = [b; quad(@(x) f(x).*sin(m*pi*x/r_end), l_end, r_end)];
  end

  a /= r_end;
  b /= r_end;
  if (l_end == 0)
    a *= 2; b *= 2;
  end
end

