	

	%%%%%%%%%%%%%%%%%%%%%%%%%
	%			%
	%    FILE NAME: wig_M	%
	%	 		%
	%%%%%%%%%%%%%%%%%%%%%%%%%



 clear
 
% A program to calculate the evolution of wave packets, and plot a movie of wigner functions
%  using the split operator propagator.
%
% parameters of the calculation

% MOLECULAR PARAMETERS

%mu=0.5*1836;              % the reduced mass
mu=2;
omega=2./mu;
alpha=mu*omega;

% X DIMENSION

xdim=12.0/sqrt(alpha);    % range of x dimension
nop=64;                   % no of points /2
xrange=nop/xdim;          %
dx=xdim/nop;
%x0=0.5/(dx*sqrt(alpha));  % shift of initial wavefunction
x0=0.5/(dx*sqrt(alpha));

x=[-1*nop:nop-1]/xrange;  % x range vector
n=length(x);

% T DIMENSION

dt=.001;                  % time step for propagation.
%niter=20;                % no of iterations
niter=70;
nstep=300;                % no of time steps per iteration
dt1=dt*nstep;             % time step of each iteration
c=-i;                     % constant copmlex -i
t=[0:niter-1]*dt1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% calculate k vector and vectors of kinetic and potential energy
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
        ak1=[0:nop-1]*pi/xdim;
        ak2=[-1*nop:-1]*pi/xdim;
        ak=[ak1 ak2];
        kine=1.*ak.*ak/(2.*mu);
        ekine=exp(c*dt.*kine);
        pote=0.5*mu*omega^2.*x.^2;     % excited states PES centered at x=0
        epote=exp(c*dt.*pote);
        epote2=exp(.5*c*dt.*pote);
        epotei2=exp(-.5*c*dt.*pote);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% initial function for propagation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%coherent wave function!
psinit=(alpha/pi)^(1/4)*exp(-0.5*alpha.*(x+x0).*(x+x0));
psisqr=psinit.*conj(psinit); norm=(psisqr*ones(n,1)*dx);
psinit=psinit/sqrt(norm);
psi=psinit;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%initialize movie making!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 %W=Wig(psi);%-------------------------------------------------
W=wig2(psi,x);
surf(W);
%axis([0 128 0 128 -5 15]); 
shading interp;
view(0,80);
%caxis([-5,12]);
%M=moviein(niter-1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% calculate psi(t) for iter=1,...,niter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        psi=psi.*epote2;         % first 1/2 step with potential
 
        for iter=2:niter         % niter "big" steps - i.e. output points
          for istep=1:nstep      % nstep "small" steps (alternating kine/pote)
             psi=ifft(psi);
             psi=psi.*ekine;
			 psi=fft(psi);
	  	     psi=psi.*epote;
          end


				psi=psi.*epotei2;      % 1/2 step backwards with potential
     
     
   		  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     	  % generate movie of wigner distributions%
     	  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
          
			%W=Wig(psi); %------------------------------
            W=wig2(psi,x);
            surf(W);
            axis([20 100 20 120 -5 5]); 
            shading interp;
            view(0,80);
  %         caxis([-5,12]);
			M(:,iter-1)=getframe;
            
         end       % --- end of loop over niter         
        % movie(M,5);%----------------------------
       %movie2avi(M,'mv_wg_coh.avi','compression','cinepak','fps',7);
         
