function [x ,eigen]=har_eigfuncFT(x0,width)
% Adapted by Shlomo Sklarz for Quantum Molecular Dynamics, Weizmann, 2001.
% This program accompanies Ch. 6 and Ch. 14 of Introduction to Quantum Mechanics:
% A Time Dependent Perspective, by David Tannor (University Science Books, 2001).
%
% The program calculates the eigenfunctions of a harmonic potential by numerically fourier transforming a 
% propagting Wave packet. 
%
%  Input vars: 
%   a) x0- displacement of initial gaussian wave packet. [x0=5/sqrt(2)]
%   b) width- of initial gaussian wave packet.          [width=sqrt(2)] 
%  Output vars:
%   a) x- the spatial dimention vector.
%   a) eigen- a matrix, the columns of which represent the eigenfunctions 
%            shifted upwords by an ammount corresponding to their respective energies.
%
%   Algorithm:
%   Wave Packet propagation using Split operator method.
figure;

% parameters of the calculation
xdim=7;                   % range of x dimension
nop=64;                   % no of points /2
dx=xdim/nop;
if nargin==0
   x0=5/sqrt(2);          % displacement of initial gaussian (the shift of the upper level)  
   width=sqrt(2);	  % width of initial gaussian;
end;
dt=.001;                  % time step in each iteration
dt1=0.05;                  % time step for recording
time=5*2*pi;              % the time for the whole process
niter=time/dt1;           % no of iteration (recording)
nstep=dt1/dt;             % no of time steps per iteration
tmax=(niter-1)*dt1;
dw=pi/tmax;
wmax=(niter-1)*dw;
c=-i;                     % constant copmlex -i
t=0;
w=1;			  % energy spacing between harmonic levels.
N=10;			  % Number of eigenfunctions to calculate.
E=w*([0:N-1]+.5);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
             
x=[-nop:nop-1]'*dx;       % x range vector
n=length(x);              % No of components of x

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% calculate k vector and vectors of kinetic and potential energy.
% the k vector should go from -pi/dx to +pi/dx
% krange=2*pi/dx
% dk=krange/(2*nop)=pi/xdim
% (the number of points in the x & k vectors is 2*nop)
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
        dk=pi/xdim;
        ak=fftshift([-nop:nop-1]'*dk);
%      
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        kine=.5*ak.*ak;
        ekine=exp(c*dt.*kine);

        pote=.5*w^2*x.*x;
        epote=exp(c*dt.*pote);
        epote2=exp(.5*c*dt.*pote);
        epotei2=exp(-.5*c*dt.*pote);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        psinit=exp(-((x+x0)/width).^2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  normalize psint;
%
   psinit=psinit/sqrt(psinit'*psinit*dx);

   psi=psinit;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%	

   eigen=psinit*ones(1,N);
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% calculate the propagation for ntime time steps
%
   t=0;
   for iter=2:niter
      for istep=1:nstep
         psi=psi.*epote2;
         psi=ifft(psi);
	 psi=psi.*ekine;
	 psi=fft(psi);
	 psi=psi.*epote2;
	 eigen=eigen+psi*(exp(-c*E*t));
	 t=t+dt;
      end

      drawnow

     %  plot psi(t) 
     % subplot(2,2,1)
     h=plot(x,abs(psi).^2,x,x.^2/50,'r');
     title('propagation of psi');
     axis ([-7 7 0 3])
     %f(iter-1) = getframe; %---------------------------------------------
   end 
   
   
 %movie2avi(f,'eigfuncFT1.avi','compression','cinepak');

%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Normalize and plot final eigenstates
%%%%%%%%%%%%%%%%%%%%%%%%%%

eigen=eigen./(ones(n,1)*(sqrt(trapz(eigen.*conj(eigen))*dx)))+ones(n,1)*E;

   plot (x,eigen);
   hold on;
   plot(x,.5*w*x.^2,'b');
   title('Harmonic Eigenfunctions Obtained by Fourier Tr. of W.P.');
   xlabel('x');  
  
   
return;