clear;
%
% This program generates wave functions using the semiclassical method
% described by M.J.Davis and E.J.Heller in ref.(1) -
% J.Chem. Phys. 75(8) 3916-3924 (1981)
% 
% The parameters in the program generate the wave function as ploted in
% FIG.8 of ref.(1). The numeration of the equations and figures in the footnotes refer
% to ref.(1). A different set of parameters for state 41 in FIG.2 is given
% at the final footnote.
%
%
omegax=1.1; % The angular frequency in the x-direction. See equation (3.1).
omegay=1.0; % The angular frequency in the y-direction. See equation (3.1).
lam=-.042;  % The lamda parameter in equation (3.1).
e=62.480;   % The energy of the system. (Since the classical trajectory initiates from 
            %  x=0 and y=0 this energy equals the initial kinetic energy)
fo=.836;    % The parameter f_x in equation (2.2). Notice m=1 hence p_x(0)=v_x(0) and p_y(0)=v_y(0).
v_x0=sqrt(2*e*fo);      % Initial velosity in the x-direction.  
v_y0=sqrt(2*e*(1-fo));  % Initial velosity in the y-direction. 
num=0;
poe(1)=0;
par=(1.036626+.954787)/2; % The quantity (omega_x'+omega_y')/2 which is the second term in equantion
                          % (2.10). omega_x' and omega_y' are derived using
                          % equantions (2.11a) and (2.11b) where M_x=38,
                          % M_y=35 and tao_c=230.325 which is the time
                          % it takes the trajectory to nearly close itself. 
%
options=odeset('MaxStep',.1,'InitialStep',.1);
[t,y]=ode45(@fad11,[0,230.325],[0 v_x0 0 v_y0 0],options); % The calculation of the classical trajectory 
                                                           % and gamma_t in equantion (2.10).See function fad11.m. 
siz=size(t,1);
n=280;                    % Roughly the number of tics in the x and y direction in which the wave function is evaluated.
x_min=min(y(:,1));        % The following quantities refer to the boundaries of the classical trajectory.
x_max=max(y(:,1));
y_min=min(y(:,3));
y_max=max(y(:,3));
dy=(y_max-y_min);
dx=(x_max-x_min);
d=max(dy,dx);
dd=d/n;
n1=round((d+4)/dd);       % The final number of tics.
mat=zeros(n1);
%
tic
 for ii=1:n1
  xx(ii)=x_min+dd*ii-2;   % The x and y coordinats for the evaluation of the wave function. 
  yy(ii)=y_min+dd*ii-2;
 end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for it=2:siz              % This loop calcutates the sum in equation (2.14). 
                          % The data is saved for the indices ii and jj (which correspond to x-y coordinates)
                          % in matrix mat(jj,ii). 
tt=(it-1)*.025            % The time of propagation. .025 is the time step.
   a=exp(i*y(it,5));
     for ii=1:n1
         b=exp(-omegax/2*(xx(ii)-y(it,1))^2+i*y(it,2)*(xx(ii)-y(it,1)));
         for jj=1:n1
             c=exp(-omegay/2*(yy(jj)-y(it,3))^2+i*y(it,4)*(yy(jj)-y(it,3)));
             mat(jj,ii)=mat(jj,ii)+a*b*c; 
         end 
     end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%     
if tt>num*10           % Each 10 time units the emerging wave function is saved in fish2
    num=num+1
    fish2(:,:,num)=mat;
	t1(num)=t(it)
    save data
end
%
end
toc
%
num=num+1              % The final formation.
fish2(:,:,num)=mat;
t1(num)=t(it);
save data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This part plots the series of the energing wave function.
%
for n=1:num
figure
matt=fish2(:,:,n);
for ii=1:n1;
  for jj=1:n1;
    matt1(ii,jj)=abs(real(matt(ii,jj))); % The absolute value of the real part of the wave function. 
  end
end
%
surf(xx,yy,matt1);view(0,90);shading interp;axis off;axis square
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% By changing the initial parameters other modes can be created.
% 
or example the 41st state in FIG.2 in ref.(1) can be obtained with the 
% aid of TABLE I in ref.(1) and by obtaining tao_c from the classical
% trajectory. To obtain this state the following changes are required in
% this program and in fad11.m :
% lam=-.11, e=9.215, fo=.840, tao_c=230.55, par=(1.0356+.9539)/2;
%
% THE END.
























