%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Simple Demo Script for Fourier decomposition using MATLAB
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%  define a spatial interval from -L to L discretized into Npts points.
L = 12;
Npts = 1024;
dx = 2*L/Npts;
xgrid = linspace(-L,L-dx,Npts);

%%%%%  define into how many Fourier modes we want to decompose this later
Nmodes = 100;

%%%%% define function to be decomposed. You may change this.

% sawtooth example

Amplitude = 1;
d=6;
funct = Amplitude*(abs(d/2 - abs(mod(xgrid,d)))-d/4);

%%%%% Fourier transform function
ftilde = ifft(fftshift(funct));  
coscoeff = 2*ftilde;
maxpop = max(abs(coscoeff));
fouriersum = zeros(size(xgrid));

%%%%%%%plot successive approximations

for nn=1:Nmodes,
  
plot(xgrid,funct,'k')
xlabel('x')
ylabel('y(x)')
xlim([-L L])
title(['cos wave no ',num2str(nn)])
hold on
  coscoeff(nn) = 2*ftilde(nn);
  nextterm = cos(pi/L*(nn-1)*xgrid);
  fouriersum = fouriersum + coscoeff(nn)*nextterm;
  if(abs(coscoeff(nn)*coscoeff(nn))>0.00001*maxpop)
    plot(xgrid,nextterm,'m')
  else  
    plot(xgrid,nextterm,'b')
  end  
    
  plot(xgrid,fouriersum,'r:')
  
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lw=6;
fs=30;

lineobj = findobj('type', 'Line');
set(lineobj,'LineWidth', lw);

patchobj = findobj('type', 'Patch');
set(patchobj,'LineWidth', lw);

textobj = findobj('type', 'axes');
set(textobj,'LineWidth', lw);
set(textobj, 'fontsize', fs);
set(textobj, 'FontName', 'Times');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  hold off
  
  
  dump = input(['press ENTER to advance'])
end



