% Programming and Numerical Methods
% MSc in Computational and Applied Physics

% Code 9: Gibbs phenomenon
% A. Meseguer (2007)


% Sawtooth function:
% f(x) = (1/pi)*(x+pi)     ; if x<0
% f(x) = (1/pi)*(x+pi)-2   ; if x>0

  clear all ; close all 

% We do the exploration for number of modes N=8,16,24,32,40,48
  iplot = 0 ; % Plot window index set to 0

  for N = [8:8:48]

% Spectral parameters  
%  N = 16 ;                     % Higher mode: exp(iNx) [a_{-N}...a_{N}]
  NX = 2*N+1 ;                 % No. grid points
  x = (-pi+(0:NX-1)*2*pi/NX)'; % Physical grid

% Fourier Transforms DFT - Matrices
  EX2F=exp(-i*[-N:N]'*x')/(2*N+1); % From X to F (X2F No-Dealiased)
  EF2X=exp(i*x*[-N:N]);            % From F to X (F2X No-Dealiased)
  DF = diag(i*[-N:N]) ;            % Fourier Derivative Matrix (in*a_{n})

% A discontinuous function
  for ix = 1:length(x)
   if x(ix) < 0 
     f(ix,1) = (x(ix)+pi)/pi ; 
   else
     f(ix,1) = (x(ix)+pi)/pi -2 ; 
   end
  end
 
% Lets do its Fourier transform (from X to F)
  a = EX2F*f ;


% We go back on a finer grid just for quality plotting
  xx = linspace(-pi,pi,400)' ;
  EF2XX = exp(i*xx*[-N:N]) ;
  invFa = EF2XX*a ;
  for ix = 1:length(xx)
   if xx(ix) < 0 
     ff(ix,1) = (xx(ix)+pi)/pi ; 
   else
     ff(ix,1) = (xx(ix)+pi)/pi -2 ; 
   end
  end

  figure(1)
  iplot = iplot + 1 ;
  subplot(6,1,iplot)
  plot(xx,real(ff),'-b','linewidth',1.5); hold on 
  plot(xx,real(invFa),'-r')
  title(['N = ', int2str(N)])


end

