function har_analytic(x0, width)
% 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 ANALYTICALLY the autocorrelation function of a displaced Gaussian in a 
% harmonic potential.  The Fourier transform of the correlation function is used to 
% calculate the electronic absorption spectrum.
figure;
clear


x0=2;             % the shift of the initial wf
dt1=0.1;                  % time step for recording
time=30;                  % the time for the whole process
niter=time/dt1;           % no of iteration (recording)
tmax=(niter-1)*dt1;
dw=pi/tmax;
wmax=(niter-1)*dw;
c=-i;                     % constant copmlex -i

             
   t=[0.0:dt1:tmax];
   w=[0.0:dw:wmax];

%  calculate the analytic expression for C(t)

   analyt=exp((-x0*x0*(1-exp(c*t)))+c*t/2);
 
%  the correlation fun. is multiplied with a gaussian in order
%  to force it to decay to zero, to avoid "ringing"
   decay=exp(-0.008*t.*t);
   Corr=decay.*analyt;

   subplot(2,2,1)
   plot(t,abs(analyt))
	title('the correlation function')
	xlabel('t')
	ylabel('|C(t)|')
        axis ([0 8*pi 0 1.0])
   
   subplot(2,2,2)
   plot(t,abs(analyt),':')
   hold on
   plot(t,decay,'--')
   plot(t,abs(Corr))
        title('correlation function + decay')
	xlabel('t')
	ylabel('|C(t)|')
 
    for n=1:niter
%        sigma1(n)=0.0;
%        sigma2(n)=0.0;     % S.F.T. for the analytic c(t)
%        for k=1:niter
%            sigma1(n)=sigma1(n)+(Corr(k)*exp(-c*w(n)*t(k)));
%            sigma2(n)=sigma2(n)+(analyt(k)*exp(-c*w(n)*t(k)));
%        end
%        sig1(n)=2*real(sigma1(n))-1.0;
%        sig2(n)=2*real(sigma2(n))-1.0;
%        spectrum(n)=dt1*sig1(n);
%        spectrum2(n)=dt1*sig2(n);

        sigma1(n)=trapz(t,Corr.*exp(-c*w(n)*t));
        sigma2(n)=trapz(t,analyt.*exp(-c*w(n)*t));
        spectrum(n)=2*real(sigma1(n));
        spectrum2(n)=2*real(sigma2(n));
    end
%
   subplot(2,2,3)
   plot(w,spectrum2)
        title('the spectrum')
        xlabel('\omega')
        ylabel('\sigma(\omega)')
        axis ([0 10 -1 11])

   subplot(2,2,4)
   plot(w,spectrum)
        title('the spectrum (no "ringing")')
        xlabel('\omega')
        ylabel('\sigma(\omega)')
        axis ([0 10 -1 11])
        
return