

	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%			                    %
	%   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.



figure; %--------------------------------
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
xdim=9;
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;
%fac=1;           % for the 3rd movie.
%fac=0.3;         % for the 1st movie.
fac=5;            % for the 2nd movie.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%             
        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(-fac*((x+x0)/width).^2);               %change fac

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  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('position',[50 50 350 250]); %--------------------------------------------------
  
   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,-abs(psi),'b',x,real(psi),x,x.^2/50,'r')
     title('propagation of psi')
     %axis ([-7 7 0 1.0])
     axis ([-xdim xdim -1 1.0])
    %f(iter - 1) = getframe; %--------------------------------------------
     
   end 

	%movie2avi(f,'numerc_5.avi','compression','cinepak'); %---------------------------------------------- 


