% this routine takes an array f and returns its spectral derivative.  it
% wants f to be periodic of length 2^n where f(1) = f(2^n+1)  And so 
% x(2^n) = 2*pi - dx, for example.  It assumes f is real-valued.
%
% to call it: >> fd = fdiff(f,L);
%
% where you've already defined f
%
function [fd] = fdiff(f,L);

N = length(f);

u = fft(f);

clear i
% The fft gives an array of length N.
% the first entry corresponds to wave number k=0.  The derivative of a
% constant function is zero and so we set v(1) to 0:
v(1) = 0;

for k = 1:(N-1)/2
  % first the positive wave numbers. The second entry of u
  % corresponds to wave number k=1, the third entry corresponds to wave
  % number k=2 and so on up to wave number floor((N-1)/2).
  % If N is odd then (N-1)/2 is an integer.  If N is
  % even then k will go up to floor( (N-1)/2 ).  Either way, we're going up
  % to the highest wave number we can resolve.  
  %
  % for positive wave numbers, multiply by i*k*(2*pi/L):
  v(1+k) = i*k*(2*pi/L)*u(1+k);
  % now for the negative wave numbers.  
  % the Nth entry of u corresponds to wave number k=-1, the N-1st entry of
  % u corresdponds to wave number k=-2 and so on.  So we multiply the N-k+1st 
  % entry of u by -i*k*(2*pi/L):
  v(N-k+1) = -i*k*(2*pi/L)*u(N-k+1);
end

% take the inverse transform.  I've cut off the imaginary parts since
% they're going to be at the level of round-off anyway.  (the true
% solution started out real-valued and will remain real.  If I wanted to
% be super-cautious, I'd keep the imaginary parts of the computed
% solution around to make sure they aren't growing on me.)
fd = real(ifft(v)); 

% test the code:
% >> N = 40; L = 1; dx = L/N; x = 0:dx:L-dx; f = exp(sin(2*pi*x)); fd = fdiff(f,L); fd(1)-2*pi
% ans =
%    2.0428e-14
% >> N = 41; L = 1; dx = L/N; x = 0:dx:L-dx; f = exp(sin(2*pi*x)); fd = fdiff(f,L); fd(1)-2*pi
% ans =
%   -4.4409e-15
  





