

	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%				                %
	%   FILE NAME: har_numeric_M	%
	%				                %
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function numeric
% written by Dr. Efrat Rosenman for Quantum Molecular Dynamics, Weizmann, 1998.
% 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 NUMERICALLY the autocorrelation function of a displaced Gaussian in 
% a harmonic potential.  The Fourier transform of the autocorrelation function is 
% used to calculate the electronic absorption spectrum.



clear


%
% A program to calculate the evolution of wave packets,
% in shifted H.O. using the split operator propagator.
%
% parameters of the calculation
xdim=7;                   % range of x dimension
nop=64;                   % no of points /2
dx=xdim/nop;

%--------------------------x0= 2/sqrt(2);
x0=3;
width = sqrt(2);  
           % displacement of initial gaussian (the shift of the upper level)  
dt=.001;                  % time step in each iteration
dt1=0.1;                  % time step for recording
time=30;                  % 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;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%             
        x=[-1*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;
        ak1=[0:nop-1]*dk;
        ak2=[-1*nop:-1]*dk;
        ak=[ak1 ak2];
%      
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

       % kine=1.*ak.*ak/4;
	    kine=1.*ak.*ak/2;
        ekine=exp(c*dt.*kine);
        pote=1*x.*x;
	    pote=.5*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;
%
   psisqr=conj(psinit).*psinit;
   norm=(psisqr*ones(n,1))*dx;
   psinit=psinit/sqrt(norm);
%
   psi=psinit;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%	
% the first step in c(t), before psi starts moving

   harm(1:n)=trapz(conj(psinit).*psi)*dx;
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% calculate the propagation for ntime time steps
%

   figure; %--------------------------------------------------
  
   for iter=2:niter
      psi=psi.*epote2;
      for istep=1:nstep
         psi=ifft(psi);
	     psi=psi.*ekine;
	     psi=fft(psi);
	     psi=psi.*epote;
      end
      psi=psi.*epotei2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  calculate  C_0j(t)  
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 harmonic(iter)=trapz(conj(psinit).*psi)*dx;
       drawnow
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  plot psi(t) 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
     subplot(2,2,1)%-----------------------------------
     plot(x,abs(psi),x,x.^2/50,'r')
     title('propagation of psi')
     axis ([-7 7 0 1.0])
     
     %f1(iter - 1) = getframe; %--------------------------------------------
     
   end 

	%movie2avi(f1,'numerc1.avi','compression','cinepak'); %---------------------------------------------- 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  plot C_0j(t) 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% comment - 
% formally w should go from -N/2 to +N/2, which correspond to
% freq from -wc(=-2*pi/(2*dt)) to +wc(=+2*pi/(2*dt)).
% but it's sufficient to go from 0 to wc since the negative half
% is symmetric. So we take this range (0 to pi/dt) and divide
% it by (nop-1) to get dw.
% for every w we have to calculate 2*real(sig(w))-1 because we
% had t only from 0 to tmax and we lack the left half (from -tmax
% to 0).
%
   t=[0.0:dt1:tmax];
   w=[0.0:dw:wmax];

%  the correlation fun. is now multiplied with a gaussian in order
%  to force it to decay to zero, to avoid "ringing"
 
   decay=exp(-0.008*t.*t);
   Corr=decay.*harmonic;

   
 
   subplot(2,2,2)
   plot(t,abs(harmonic),'b:')
   hold on
   plot(t,decay,'k:')
   plot(t,abs(Corr),'r-')
   title('correlation function')
   xlabel('t')
   ylabel('|C_0j(t)|')
   axis ([0 8*pi 0 1.0])
        
      
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  calculate the analytic expression for C(t)
%
   analyt=exp((-((x0)^2)/2*(1-exp(c*t)))+c*t/2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   Corr2=decay.*analyt;
%
    for n=1:niter
       sigma1(n)=trapz(t,Corr.*exp(-c*w(n)*t));
       sigma2(n)=trapz(t,Corr2.*exp(-c*w(n)*t));
       spectrum(n)=2*real(sigma1(n));
       spectrum2(n)=2*real(sigma2(n));
       
    end
%

 	
   subplot(2,2,3)
   plot(w,spectrum,'g-')
   title('the numeric spectrum of the decayed correlation function')
   axis ([0 15 -1 4])
   %axis ('on')
       	
	
  
    subplot(2,2,4)
    plot(w,spectrum2,'r-')
    title('Analytic spectrum')
    xlabel('w')
    axis ([0 15 -1 4])


return;




