
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%				                %
	%     FILE NAME: i2_cars_M  	%
	%				                %
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%	



function i2_cars

% written by Dr. Jim Faeder for Quantum Molecular Dynamics, Weizmann, 2000.
% This program accompanies Ch. 13 of Introduction to Quantum Mechanics:
% A Time Dependent Perspective, by David Tannor (University Science Books, 2001).
% program for simulation of femtosecond CARS spectroscopy.  The potential
% curves correspond to the X and B states of I_2.  The program calculates
% the third order wavefunction, psi_123, produced by a sequence of three
% laser pulses corresponding to a pump, dump, pump (or probe) sequence,
% taking the wavefunction from X -> B -> X -> B.  The third order polarization
% is calculated as the overlap of psi_123 with the initial vibrational 
% state on the X surface, psi_0.  The program allows adjustment of the
% time delay between pump and dump, as well as the time delay between
% dump and probe.  The carrier frequencies of the three pulses can also
% be adjusted, although experimentally omega_1 is usually equal to omega_3.
% The default value of the pulse widths are 75 fs FWHM of the pulse
% INTENSITY.  The two different sample input files contrast the qualitative 
% difference in the polarization produced depending on whether the third 
% (probe) pulse is below or above the dissociation energy of the B state.  

% examples of input parameters:
%   third (probe) pulse below disociation limit: i2_cars(50,300,580,625,580,3);
%     "   (  "  )   "    above    "       "  i2_cars(50, 75, 586, 800, 586, 2);
% define various conversion constants
AU_TO_EV=27.2113961;
NM_TO_AU=45.56335;
AU_TO_CM=2.194746e5;
AU_TO_ANG=0.529177249;
AU_TO_FS=0.0241889;
AMU_TO_AU=1.8228885e3;
conv_r= AU_TO_ANG;
conv_en=AU_TO_CM;
conv_t= AU_TO_FS;

te =50;
Tau = 300;
w0_1 = 580;
w0_2 = 625;
w0_3 = 580;
v_init = 3;


% Turn on and off plotting of the wavepacket wavefunction
plot_wave=1;

% Read potentials and and eigenstates from file
pot_file='i2_256.mat';
sprintf('loading potential data from %s', pot_file)
pot= load(pot_file, 'ngrid', 'mass', 'rmin', 'rmax', ...
          'V1', 'V2', 'E1', 'E2', 'evecs1', 'evecs2'); 
ngrid= pot.ngrid;
rmin= pot.rmin;
rmax= pot.rmax;
mass= pot.mass;
V1= pot.V1;
V2= pot.V2;
E1= pot.E1;
E2= pot.E2;
evecs1= pot.evecs1;
evecs2= pot.evecs2;
clear pot;

% set up x grid
dr= (rmax - rmin)/(ngrid-1);
r= rmin:dr:rmax;

% quantum number of initial ground state wavefunction 
psi_init= (evecs1(:,v_init+1))';
psi_init= psi_init/sqrt(dot(psi_init,psi_init));
E_init= E1(v_init+1);

% pulse parameters
FWHM= 75.0/AU_TO_FS; % This is FWHM of pulse intensity, not amplitude
te =  te/AU_TO_FS; % Time delay between pump and Stokes pulses
Tau = Tau/AU_TO_FS;% Time delay between Stokes and probe pulses
w0_1= NM_TO_AU/w0_1; % Pulse 1 frequency
w0_2= NM_TO_AU/w0_2; % Pulse 2 frequency
w0_3= NM_TO_AU/w0_3; % Pulse 3 frequency
phi_2= 0;
phi_3= 0;
t0_1=  2*FWHM; % center of laser pulse 1
t0_2=  t0_1 + te;
t0_3=  t0_2 + Tau;
% Pulse intensities in atomic units
e0_1= 0.002;     
e0_2= 0.002;     
e0_3= 0.002;     
% Transition dipole in atomic units
mu12= 0.181;    

% integration parameters
dt      =     0.1 /AU_TO_FS;
t_sample=    1   /AU_TO_FS;
t_end   =   t0_3 + 300/AU_TO_FS;
update_freq= 10; % update_freq*t_sample is the interval at which plots updated
n_steps_sample= round(t_sample/dt);
n_sample= round(t_end/t_sample);
t_end= n_sample*t_sample;

% pulse function
t= 0:dt:t_end;
alpha= 2*log(2)/FWHM^2;
pulse1= e0_1*exp(-alpha*(t-t0_1).^2).*exp(-i*w0_1*(t-t0_1));
pulse2= e0_2*exp(-alpha*(t-t0_2).^2).*exp( i*w0_2*(t-t0_2))*exp(i*phi_2);
pulse3= e0_3*exp(-alpha*(t-t0_3).^2).*exp(-i*w0_3*(t-t0_3))*exp(i*phi_3);

cpu_start= cputime;
clock_start= clock;

% set up k grid
dk= 2*pi/(ngrid*dr);
k=[[0:ngrid/2-1],[-ngrid/2:-1]]*dk;

% set up kinetic energy and k-space propagator
kin= (0.5*(k.^2)/mass).';
sprintf('[Rmin,Rmax]:[%.1f, %.1f] Ang   ngrid: %d Tmax=%.0f wavenumbers', ...
          conv_r*rmin, conv_r*rmax, ngrid, max(kin)*conv_en)

% set up propagators
u_kin= exp(-i*dt*kin);
u_kin_half= exp(-i*0.5*dt*kin);
u_pot1= exp(-i*dt*(V1.'));
u_pot2= exp(-i*dt*(V2.'));

psi_0=    psi_init.';          % unperturbed wavefunction
psi_1=    zeros(ngrid, 1);     % First order wavefunction
psi_12=   zeros(ngrid, 1);     % second order wavefunction
psi_123=  zeros(ngrid, 1);     % third order wavefunction

% Plot potential curves and initial wavefunction
if (plot_wave)
   h=figure(1); %------------------------adding handle for the plot
   clf;
   set(h,'position',[50, 50,449,500])%----------------------------
   subplot(212);
   gs_plot=plot( r*conv_r, V1*conv_en);
   gs_min= conv_en*min(V1);
   gs_max= conv_en*V1(ngrid);
   axis manual;
   %axis([conv_r*min(r) conv_r*max(r) gs_min gs_max]);%% changed for the movie---------------------
   axis([conv_r*min(r) 5 gs_min gs_max]); %%dido-------------------------------------------
   hold on;

   subplot(211);
   es_plot=plot( r*conv_r, V2*conv_en);
   es_final= conv_en*(w0_1-w0_2+w0_3+E_init);
   es_final
   es_min=   conv_en*min(V2);
   es_max= es_final + 0.2*(es_final-es_min);
   es_max
   axis manual;
   %axis([conv_r*min(r) conv_r*max(r) es_min es_max]);%% changed for the movie-------------------------
   axis([conv_r*min(r) 5 es_min es_max]);%%dido-------------------------------------------------------
   hold on;

   scale_gs= 0.2*(gs_max-gs_min)/max(abs(psi_init));
   scale_es= 0.2*(es_max-es_min)/max(abs(psi_init));
   wave_scale= 5;

   subplot(212);
   psiplot_gs= ...
      plot( r*conv_r, scale_gs*abs(psi_init)+ E_init*conv_en, ... 
            r*conv_r, 0*psi_init+ (E_init+w0_1-w0_2)*conv_en, ...
            r*conv_r, 0*psi_init+ (E_init+w0_1-w0_2)*conv_en);
		  
   subplot(211);
   psiplot_es= plot( r*conv_r, 0*psi_init+ (E_init+w0_1)*conv_en,...
                     r*conv_r, 0*psi_init+ (E_init+w0_1)*conv_en,...
                     r*conv_r, 0*psi_init+ (E_init+w0_1-w0_2+w0_3)*conv_en,...
                     r*conv_r, 0*psi_init+ (E_init+w0_1-w0_2+w0_3)*conv_en);
end

in = 1;  %-----------------------------------------------------------

tot_steps=0;
t=0;
n_sample=0;
max_steps= round(t_end/dt);
start_123= round((t0_3 - 2*FWHM)/t_sample)*n_steps_sample;
end_123  = round((t0_3 + 2*FWHM)/t_sample)*n_steps_sample;
while (tot_steps<max_steps)
   % propagate wavefunctions 1/2 step forward
   U= u_kin_half;
   psi_1=    ifft(U.*fft(psi_1));
   psi_12=   ifft(U.*fft(psi_12));
   psi_123=  ifft(U.*fft(psi_123));

   % Propagate until next sample time 
   for step=1:n_steps_sample
      % potential step
      psi_0=    psi_0*exp(-i*dt*E_init);
      psi_1=    u_pot2.*psi_1;
      psi_12=   u_pot1.*psi_12;
      psi_123=  u_pot2.*psi_123;
      % kinetic step
      if (step==n_steps_sample)
         % apply only 1/2 kinetic step when wavefunction is sampled
         U= u_kin_half;
      else
         U= u_kin;
      end
      psi_1=    ifft(U.*fft(psi_1));
      psi_12=   ifft(U.*fft(psi_12));
      psi_123=  ifft(U.*fft(psi_123));

      % Add contribution due to field interactions
      psi_1=    psi_1    + i*dt*(pulse1(tot_steps+1)*mu12.*psi_0);
      psi_12=   psi_12   + i*dt*(pulse2(tot_steps+1)*mu12.*psi_1);
      psi_123=  psi_123  + i*dt*(pulse3(tot_steps+1)*mu12.*psi_12);
      t=t+dt;
      tot_steps=tot_steps+1;
   end
   n_sample=n_sample+1;

   pop1= dot(psi_1,psi_1);
   pop_1(n_sample)= pop1;
   if (pop1)
      r_1(n_sample)= real(dot(psi_1, (r.').*psi_1)/pop1);
   end

   pop12= dot(psi_12,psi_12);
   pop_12(n_sample)= pop12;
   if (pop12)
      r_12(n_sample)= real(dot(psi_12, (r.').*psi_12)/pop12);
   end

   if (tot_steps>=start_123)
      pop123= dot(psi_123,psi_123);
      pop_123(n_sample)= pop123;
      if (pop123)
         r_123(n_sample)= real(dot(psi_123, (r.').*psi_123)/pop123);
         P_3(n_sample)= 2*mu12*dot(psi_0,psi_123);
      end
   end

   if (tot_steps==end_123)
      % compute coefficients of psi_123 in eigenvector basis once field has 
      % become negligible
      psi_123_ev= (evecs2'*psi_123).*exp(i*E2'*t);   
   end

   if (plot_wave & mod(tot_steps/n_steps_sample, update_freq)==0)
      set(psiplot_gs(2), 'YData', wave_scale^2*scale_gs*abs(psi_12)+ ...
                                   (E_init+w0_1-w0_2)*conv_en);
      set(psiplot_gs(3), 'YData', wave_scale^2*scale_gs*real(psi_12)+ ...
                                   (E_init+w0_1-w0_2)*conv_en);

      set(psiplot_es(1), 'YData', wave_scale*scale_es*abs(psi_1)+ ...
                                   (E_init+w0_1)*conv_en);
      set(psiplot_es(2), 'YData', wave_scale*scale_es*real(psi_1)+ ...
                                   (E_init+w0_1)*conv_en);
      set(psiplot_es(3), 'YData', wave_scale^3*scale_es*abs(psi_123)+ ...
                                   (E_init+w0_1-w0_2+w0_3)*conv_en);
      set(psiplot_es(4), 'YData', wave_scale^3*scale_es*real(psi_123)+...
                                   (E_init+w0_1-w0_2+w0_3)*conv_en);
      s1= sprintf('(t-t_1)= %.0f fs', conv_t*(t-t0_1));
      s2= sprintf('  (t-t_2)= %.0f fs', conv_t*(t-t0_2));
      s3= sprintf('  (t-t_3)= %.0f fs', conv_t*(t-t0_3));
      title(strcat(s1,s2,s3));
      drawnow;%------------------------------------------
      
     M(in) = getframe(h); %-------------------------------------------------
     in = in+1;  %--------------------------------
   end
 
end
%movie2avi(M,'i2_crs1.avi','compression','cinepak');%----------------------------

%compute polarization from eigenvector coeffiecents analytically
p= mu12*(evecs2'*(psi_init.')).*psi_123_ev;
dE= E2-E1(v_init+1);

for n=(1:n_sample)
   P_3_ev(n)= 2*sum(p.*exp(-i*dE'*(n*t_sample)));
end

%Plot pulses and populations
figure(2);
clf;
subplot(411);
plot((0:dt:t_end)*conv_t, abs(pulse1+pulse2+pulse3).^2);
ylabel('Field intensity');
subplot(412);
plot((t_sample:t_sample:t_end)*conv_t, pop_1);
ylabel('B state population');
subplot(413);
plot((t_sample:t_sample:t_end)*conv_t, pop_12);
ylabel('Raman population');
subplot(414);
plot((t_sample:t_sample:t_end)*conv_t, pop_123);
ylabel('CARS population');



figure(3);
plot(((t_sample:t_sample:t_end)-t0_3)*conv_t, ...
       [abs(P_3); abs(P_3_ev)]);
title('Polarization, short time behavior/ numerical + analytical');
ylabel('Polarization envelope');
legend('numerical','analytical',0);



figure(4);
for n=(1:10000)
   P_3_ev(n)= 2*sum(p.*exp(-i*dE'*(n*10*t_sample)));
end
plot(((1:10000)*10*t_sample-t0_3)*conv_t, abs(P_3_ev));
title('Polarization, long time behavior/ analytical')
ylabel('Polarization envelope');

p_final= pop_123(length(pop_123));

cpu_end= cputime;
cpu_elapsed= cpu_end-cpu_start
time_elapsed= etime( clock, clock_start) 


return;


