function [x ,eigen]=har_eigfuncRlx(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 using 
% the 'imaginary time' or  'relaxation' method.
% This method can only extract the lowest eigenfunction existing in the
% initial wave packet therefor in real situations one can only get the 
% eigenfunctions one by one beggining from the lowest. We first find 
% the lowest eigenfunction by relaxation. Then we extract this component 
% out of the initial Wave Packet. The new initial wave packet can now be 
% relaxed again to obtain the next lowest eigenfunction and so on.  
%
%  In the following program we simultainiously obtain all eigenfunctions, 
% phi_n, by assumming for each that all lower eigenfunctions are 
% already known. Using the analytic eigenfunctions we construct psi(:,n), 
% the n'th column of which contains the initial wave packet with all 
% phi_m m<n extracted out. These are now simultainiously propagated to 
% obtain each of the n eigenfunctions.  
%
%  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;
w=1;			  % energy spacing between harmonic levels.
N=10;			  % Number of eigenfunctions to calculate.
E=w*([0:N-1]+.5);
dt=.001;                  % time step in each iteration
dt1=0.1;                  % time step for recording
time=10*2*pi/w;            % 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=-0.04;                   % constant copmlex -i
t=0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
             
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)*ones(1,N);

        pote=.5*w*x.*x*ones(1,N)-ones(n,1)*E;
        epote=exp(c*dt.*pote);
        epote2=exp(.5*c*dt.*pote);

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

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

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


   psi=psi*ones(1,N);

   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Compute the analytical harmonic eigenfunctions:
%   phi_n=(w/pi)^.25/sqrt(n!2^n)*H_n(y)exp(-y^2/2) where y=sqrt(w)x.
%
%  Note: Hermite recurrence formula:
%  		H(n+1)=2xH(n)-2nH(n-1)
%  		phi are the analytical harmonic eigenfunctions
%  We extract out of the initial WP psi(:,n) all the components of 
%  eigenfunctions lower than n. n now becomes the lowest energy eigenfunction in 
%  psi(:,n), and can be found found by the relaxation method. 


   y=sqrt(w)*x;
   H=zeros(n,N);
   H1=1; H2=0;
   for idx=0:8,
      if idx==0, H(:,idx+1)=1; 
      else H(:,idx+1)=2*y.*H1-2*(idx-1)*H2; end;
      phi=(w/pi)^.25*H(:,idx+1).*exp(-y.^2/2)/sqrt(2^idx*factorial(idx));

      curr=idx+2;
      % a is the overlap of the initial WP with the idx eigenfunction. 
      a=trapz(x,psi(:,curr).*conj(phi));
      % extract the n-1'th component out of psi(:,curr-1).
      psi(:,curr)=psi(:,curr-1)-a*phi;
      
      H2=H1; H1=H(:,idx+1);
      phin(:,idx+1)=phi+E(idx+1);
   end;

   
   plot(x,.5*w*x.^2,'b',x,phin);
   hold on;
   title('a) Harmonic Eigenfunctions- analytic');
   axis([-7 7 0 10]);
 
   
   figure;%--------------------------------------------------------
   plot(x,.5*w*x.^2,'b',x,psi+ones(n,1)*E);
   title('b) Initial Wave Packet\newline with lower eigenfunctions extracted.');
   axis([-7 7 0 10]);
   
   
   figure;%--------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% calculate the propagation for ntime time steps
%
   for iter=2:niter
      for istep=1:nstep
         psi=psi.*epote2;
         psi=ifft(psi);
	 psi=psi.*ekine;
	 psi=fft(psi);
	 psi=psi.*epote2;
      end

      drawnow

     %  plot psi(t) 
     plot(x,real(psi))
     title('c) propagation of psi')
     axis([-8,8,-0.4,0.8]);%---------------------------------------------
     
     f(iter-1) = getframe;%------------------------------------------------
   end 

  % movie2avi(f,'eigfuncRlx3.avi','compression','cinepak');%-----------------
   
%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Normalize and plot final eigenstates
%%%%%%%%%%%%%%%%%%%%%%%%%%
eigen=psi;


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

   figure;%---------------------------------------------------------------
   plot (x,eigen);
   hold on;
   plot(x,.5*w*x.^2,'b');
   axis([-7 7 0 10]);
   title('d) Harmonic Eigenfunctions- numeric\newline Relaxation method');
   xlabel('x');  
   
return;