function [y,yhat] = swell(x,M)  

% [0,1] periodic 'swell' function
% Returns gridpoint values y and Fourier coeffs yhat at gridpoints x
% uniformly spaced over [0,1).

  alpha = 1/3;
  a = 2*alpha/(1+alpha^2);
  y = 1./(1 - a*cos(2*pi*x));
  yhat = alpha.^abs(M)*(1 + alpha^2)/(1 - alpha^2);