% This program plos the wigner represantation of a wavepacket of a particle
%  in half a box.

clear;

hbar=1;
m=1;

L=5;
Nx=2^8;                         % Number of grid points.
dx=2*L/Nx;
x=[0:Nx-1]'*dx;
Neigen=10;                      % Number of eigenstates included
P=zeros(Nx,Neigen);             % Eigenstates  
an=zeros(Neigen,1);             % Eigenstate coefficients
En=zeros(Neigen,1);             % Spectrum (Eigenvalues)


Nper=1;                         % Number of periods to propagate. 
tau=16*m*L^2/hbar/pi;           % time for one period 
T=tau*Nper;                     % total propagation time.
% Nt=300;
Nt=60;
dt=T/Nt;
t=[0:dt:T]';
C=zeros(Nt,1);

for n=1:Neigen;                 % preperation of eigenstates and coeff's
   P(:,n)=L^-.5*sin(n*pi*x/2/L);
   En(n)=(n*pi*hbar)^2/8/m/L^2;
   if n==2, an(n)=2^-.5;
   elseif mod(n,2)==0, an(n)=0;
   else an(n)=4*2^.5*(-1)^((n+1)/2)/(n+2)/(n-2)/pi;
   end;
     
end; 
psi0=sum(P*an,2);                % psi0;
in=1;
for tindx=1:Nt +1                 % Propagation loop
   t=(tindx-1)*dt;
   psin=P*(an.*exp(-(i*En*t)));     
   psi=sum(psin,2);
   % C(tindx)= psi0'*psi*dx;% autocorrelation function.
      W=wig2(psi',x);      % calculate the wigner 
      surf(W(15*Nx/32:17*Nx/32,:));
      view(0, 90);
      axis([0 250 2 16]);
      caxis([-0.0002 0.0012]);
      shading interp;
      drawnow;
       if (tindx==1)       % to delay of the 1st plot.
          for (ind=1:3)
              M(in)=getframe;
              in=in+1;
          end
      end
      %M(in)=getframe;%-------------------------------------------------------
      %in=in+1;
end;

%movie(M);%--------------------------------------------------------------