% THE WAVE SPLITS TO 2 AFTER HITTING THE BARRIER.
%


clear all
%
% A program to calculate the evolution of a wave packet hitting  the Eckart barrier
% using the split operator propagator
% written by Yossi Elran for Quantum Molecular Dynamics, Weizmann, 1998.
% This program accompanies Ch. 7 of Introduction to Quantum Mechanics:
% A Time Dependent Perspective, by David Tannor (University Science Books, 2001).
%
% Parameters of the calculation
%
xdim=60;                   % range of x dimension 
%nop=2^10;                 % no. of points  
nop=2^7;
xrange=nop/xdim;
x0=10;                     % the starting point (50) 
dt=.0005;                  % time step in each iteration 
dt1=0.5;                   % time step for recording 

time=20;                   % the time for the whole process (45)

%time 1;
niter=time/dt1;             % no. of iteration 
nstep=dt1/dt;               % no. of time steps per iteration 
tmax=(niter-1)*dt1;         % maximum time of propagation
dw=pi/tmax;                 % see NUMERICAL RECIPES about fft
wmax=(niter-1)*dw;
c=-i;                       % constant complex -i
t=0;

k0=2; %7;	     	        % initial momentum -------------changed====================
a=1;	  	                % Eckart barrier parameter

%  to create 3 different movies in the web, v0 and v1 got three different
%  values. The two other options are:
%  v0= 35/5         v1= 20/5
%  v0= 35/22        v1= 20/22

v0=35/9;		            % barrier height (initially I chose 100)
v1=20/9;		            % potential on right side of barrier (initially 0.3 for assymmetric, 0.0 for symmetric)


q0=30;		                % barrier coordinate
dx=xdim/nop;
alpha=0.2;         	        % initial Gaussian width (0.2)  
enstep=0.62;                % Numerical energy step for creating energy vector 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
	x=[0:nop-1]/xrange;      % x range vector
	n=length(x);             % No of components of x

%
% calculate k vector and vectors of kinetic and Eckart potential energy
%
	dk=2*pi/xdim;
        ak=fftshift([-nop/2:nop/2-1])*dk;

	kine=1.*ak.*ak/2;
	ekine=exp(c*dt.*kine);
	xi=-exp((x-q0)/a);
	pote=(-v1.*xi./(1-xi)-v0.*xi./((1-xi).*(1-xi)));
	epote=exp(c*dt.*pote);
	epote2=exp(.5*c*dt.*pote);
	epotei2=exp(-.5*c*dt.*pote);

%	Initial normalized Gaussian 

     psinit=exp(-alpha.*(x-x0).*(x-x0)-c.*k0.*(x-x0));

    	psinit=psinit/sqrt(trapz(conj(psinit).*psinit)*dx);
	psi=psinit;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Initial propagateion separating reflected and transmitted parts
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

     corrfun(1)=trapz(conj(psinit).*psi)*dx;

%
% Take half time step with potential
%
      psi=psi.*epote2;
%
% Calculate propagation for ntime time steps 
%    	


 W=wig2(psi,x);
 surf(W);
 axis([0 120 0 120]);
 caxis([-0.0007    0.0022]);
 shading interp;
 view([0,0,90]);
 
      for iter=2:niter   
	   for istep=1:nstep
		 psi=fft(psi);
		 psi=psi.*ekine;
		 psi=ifft(psi);
		 psi=psi.*epote;
	   end
	psi=psi.*epotei2;
%
% Calculate Correlation function C(t)
%

	%corrfun(iter)=trapz(conj(psinit).*psi)*dx;%---------------------------------
	
%
%plot(x,abs(psi).^2,x,pote);axis([0,60,0, .25]);drawnow; 
            W=wig2(psi,x);%------------------------------
            surf(W);
            axis([0 120 0 120]);
            caxis([-0.0007    0.0022]);
            shading interp;
            view([0,0,90]);
            drawnow;
		 
     end