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

% Code 8: aliasing effects
% A. Meseguer (2007)

  clear all 
  close all 

% Spectral parameters  
  N = 8 ;                     % 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 simple function
 f = cos(5*x) ; 

% Lets do its Fourier transform (from X to F)
 a = EX2F*f ;

% Now we calculate f' with Fourier Transform

  invec = i*[-N:N]' ;   % a vector containing the i*n factors
  b = invec.*a ;        % product component x component (.*)

% If that is correct, the inverse Fourier of b should look like
% -5*sin(5x)...lets check

  fprime = EF2X*b ;
  plot(-pi:.01:pi,-5*sin(5*(-pi:.01:pi)),'-b') ; hold on % Exact in blue
  plot(x,fprime,'.r') % Approx. in red
 daspect([.25 1 1])

% Now we calculate the Fourier coeffs of f*f' 
  
  b2 = EX2F*(f.*fprime) ; % We do the physical product point to point (.*)

% If this is correct, its inverse Fourier should look like
% -(5/2)*sin(10*x)

  ftimesfprime = EF2X*b2 ;
  x2 = [-pi:.01:pi] ;

figure(2)
  plot(x2,-(5/2)*sin(10*x2),'-b') ; hold on % Exact in blue
  plot(x,ftimesfprime,'.r') % Approx. in red

% But it should be zero!!!!!!!

% Actually, if we look at b2 coeffs, we observe that n=+/-7 are
% filled:

%b2 =
%
%   0.0000 - 0.0000i
%  -0.0000 - 1.2500i
%  -0.0000 + 0.0000i
%  -0.0000 + 0.0000i
%   0.0000 + 0.0000i
%  -0.0000 - 0.0000i
%   0.0000 + 0.0000i
%   0.0000 - 0.0000i
%   0.0000 - 0.0000i
%   0.0000 + 0.0000i
%   0.0000 - 0.0000i
%  -0.0000 + 0.0000i
%   0.0000 - 0.0000i
%  -0.0000 - 0.0000i
%  -0.0000 - 0.0000i
%  -0.0000 + 1.2500i
%   0.0000 + 0.0000i

% So: it is actually the function: 
% -1.25*i*exp(-i*7*x) + 1.25*i*exp(i*7*x)
   plot(x2,-1.25*i*exp(-i*7*x2) + 1.25*i*exp(i*7*x2),'-m') ; 
% Aliased function plotted in magenta
 daspect([.25 1 1])
