% The program was written by Yamit Maymon on 2002. (David Tannor's group).

%  This program plots a movie of a wavepacket 
%  that changes into an eigenfunction in time. The color 
%  of each point represents it's phase.
%  (In a scale of 0 to 2*pi).

clear;
xdim=7;
nop=128;               %number of points/2
dx=xdim/nop;
E=3.5;                 %energy level: 1.5, 3, 7.5  for 3 different movies.

x0=5/sqrt(2);
width=sqrt(2);

dt=0.01; %.001 
t=0;   % initial time.
w=1;
Nperiod=5; 
x=[-nop:nop-1]'*dx;
n=length(x);
T=((2*pi)/w)*Nperiod;                            % final time.
T=T-mod(T,dt);
 
pote=.5*x.*x;
epote=exp(-i*dt.*pote);
epote2=exp(.5*(-i)*dt.*pote);

 psinit=exp(-((x+x0)/width).^2);
 psinit=psinit/sqrt(psinit'*psinit*dx);
 psi=psinit;
 dk=pi/xdim;
 ak=fftshift([-nop:nop-1]'*dk);
 kine=.5*ak.*ak;
 ekine=exp(-i*dt.*kine);
 eigen=psinit;
 f=figure(1);
 set(f,'position',[0 50 400 350]);
 plot(x,abs(psi).^2*dx);
 mv_indx=1;                         %index for the movie.
 
      for t=0:dt:T 
         psi=psi.*epote2;
         psi=ifft(psi);
	     psi=psi.*ekine;
	     psi=fft(psi);
	     psi=psi.*epote2;
	     eigen=eigen+psi*(exp(i*E*t));      
         eigen_plot=eigen./(ones(n,1)*(sqrt(trapz(eigen.*conj(eigen))*dx)));%--normalizing the plot.
	     if mod(t/dt,10)==0,
            cla;                     % clearing the screen.
%	        myplot2(x,abs(eigen_plot),phase(eigen_plot'));
            myplot(x,real(eigen_plot),phase(eigen_plot'));
            axis([-8 8 -1 1]);
	        drawnow;
            t;
         	M(mv_indx)=getframe;   % for the movie.
            mv_indx=mv_indx+1;
       end
          
      end
     % movie2avi(M,'eig_three_and_half.avi','compression','cinepak');