%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Simple Demo Script for Fourier decomposition of the Gaussian wavepacket 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.

%Gaussian
Amplitude = 1;
d=2;
funct = Amplitude*exp(-xgrid.*xgrid/d/d);

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

K = pi/dx;
dk = 2*K/Npts;
kgrid = linspace(-K,K-dk,Npts);
ftildeplot = fftshift(ftilde);
kplot = 5;
%%%%%%%plot successive approximations

for nn=1:Nmodes,
subplot(3,1,1)  
ms = 40;
plot(kgrid,abs(ftildeplot),'k.','MarkerSize',ms)
hold on
  if(abs(coscoeff(nn)*coscoeff(nn))>0.00001*maxpop)
     plot(kgrid(Npts/2 + nn),abs(ftildeplot(Npts/2 + nn)),'m.','MarkerSize',ms)
  else
     plot(kgrid(Npts/2 + nn),abs(ftildeplot(Npts/2 + nn)),'b.','MarkerSize',ms)  
  end
hold off
axis tight
xlim([0 kplot])
ylim([0 1.2*max(abs(ftildeplot))])
ylabel('~g(k)')
zlabel('k')
title(['cos wave no ',num2str(nn-1)])
  
subplot(3,1,2:3)  
plot(xgrid,funct,'k')
xlabel('x')
ylabel('y(x)')
xlim([-L L])
hold on
  coscoeff(nn) = (1 + (nn>1))*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:')
 ylim([-1.2 1.2]) 
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lw=4;
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



