clear
%
% 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=250;                   % range of x dimension 
nop=2^10;                   % no. of points  
xrange=nop/xdim;
x0=90;    %45                    % the starting point (50) 
dt=.0005;                   % time step in each iteration 
dt1=0.05;                   % time step for recording 
time=9;        %12            % 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=7;	     	         % initial momentum
a=1;	  	             % Eckart barrier parameter
v0=75;		             % barrier height (initially I chose 100)
%v1=20;		             % potential on right side of barrier (initially 0.3 for assymmetric, 0.0 for symmetric)
v1=10;
q0=125;		%70          % 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 
%    	

   
     for iter=2:niter   
	   for istep=1:nstep
		 psi=fft(psi);
		 psi=psi.*ekine;
		 psi=ifft(psi);
		 psi=psi.*epote;
	   end
%
% Take half time step backwards with potential
%			
	psi=psi.*epotei2;
%
% Calculate Correlation function C(t)
%

	corrfun(iter)=trapz(conj(psinit).*psi)*dx;
	drawnow
   
%
% plot psi(t)
%	
%	
	xtick=[-12:3:12];
	ytick=[0:1]*.2;
	h1=figure(1);
    set(h1,'position',[0, 50,449,399])%----------------------------
	subplot(2,1,1)
	plot(x,pote./25,'r',x,abs(psi))
		%title('x propagation of \psi');%------------------------
        title('Position Space');%--------------------------------
		%xlabel('x');     %---------------interrupts the movie
		ylabel('|\psi|')
		axis([0 xdim 0 1.0])
        
        
        
	subplot(2,1,2)
	plot(fftshift(ak),abs(fftshift(fft(psi)).*dx./sqrt(2*pi)))
		%title('k propagation of \psi') %------------------------
		title('Momentum Space');%--------------------------------
        xlabel('k')
		ylabel('|\psi|')
        M1(iter-1)=getframe(h1);%-----------------------------------------
     end
     
    % movie2avi(f1,'traa1.avi','compression','cinepak');
     
%	
% Plot C(t)
% 
	t=[0.0:dt1:tmax];
	w=[0.0:dw:wmax];
	figure(2)
	plot(t,abs(corrfun),'r-')
		title('whole WP correlation function')
		xlabel('t')
		ylabel('|C(t)|')
		axis([0 2 0 1.0])

   %named- traa2.eps (using: File--> export)
   %print -djpeg traa2.jpeg;     
        
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Propagation of transmitted ad reflected parts seperately
%  ALL THE NEXT STEPS In principal ARE PERFORMED TWICE - 
%           ONCE FOR THE REFLECTED WP AND ONCE FOR THE TRANSMITTES WP
%  HOWEVER (!) In practice for efficiency we propagate psi including 
%         both parts of the wavepacket (transmitted and reflected) but
%         separate them out when calculating the auto correlation.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	psirinit=psi.*(x<q0);
        psitinit=psi.*(x>q0);
%
%
% The first step in C(t) befor psi starts moving
%
        corrfunr=zeros(1,niter);
        corrfunt=zeros(1,niter);

	corrfunr(1)=trapz(conj(psirinit).*psi)*dx;
	corrfunt(1)=trapz(conj(psitinit).*psi)*dx;

%
%
% Take half time step with potential
%

    	psi=psi.*epote2;

%
% Calculate propagation for ntime time steps
%    	
     for iter=2:niter
	for istep=1:nstep
      	   psi=fft(psi);
	   psi=psi.*ekine;
	   psi=ifft(psi);
	   psi=psi.*epote;
	end
%
% Take half time step backwards with potential
%			
	psi=psi.*epotei2;
%
% Calculate C(t)
%
       corrfunr(iter)=trapz(conj(psirinit).*psi)*dx;
       corrfunt(iter)=trapz(conj(psitinit).*psi)*dx;
       drawnow

%
% plot psi(t)
%	
 psit=psi.*(x>q0);
 psir=psi.*(x<q0);
	
	h3 = figure(3);%-------------------------------------------
    set(h3,'position',[0, 50,449,399])%----------------------------
	subplot(2,1,1)
	plot(x,abs(psir),'r:', ...
             x,abs(psirinit),'k-', ...
             x,abs(psit),'b-.', ...
             x,abs(psitinit),'k--', ...
             x,pote/25,'r')
		legend('reflected WP','original reflected WP','transmitted WP','original transmitted WP',2)
		%title('propagation of \psi_r and \psi_t');
		title('Position Space');%--------------------------------
       % xlabel('x') %---------------------------interrupts the movie
		ylabel('|\psi|')
%		axis([xr(1) xt(nop) 0 1.0])
        %f3(iter -1) = getframe;%-----------------------------------------
	subplot(2,1,2)
	plot(fftshift(ak),fftshift(abs(fft(psir)).*dx./sqrt(2*pi)),'r:', ...
             fftshift(ak),fftshift(abs(fft(psit)).*dx./sqrt(2*pi)),'b-.')
		legend('reflected WP','transmitted WP')
		%title('k propagation of \psi_r and \psi_t')
        title('Momentum Space');%--------------------------------     
		axis([(-nop/2).*2.*pi./xdim (nop/2).*2.*pi./xdim 0 1])
		xlabel('k')
		ylabel('|\psi|')
        
   % M3(iter -1) = getframe(h3);%-----------------------------------------
     end
     
     %   movie2avi(M3,'traa3.avi','compression','cinepak');
     
%     
% Plot C(t)
% 
	t=[0.0:dt1:tmax];
	w=[0.0:dw:wmax];
	figure(4)
	subplot(2,1,1)
	plot(t,abs(corrfunr),'r:')
		title('reflected correlation function')
		xlabel('t')
		ylabel('|C_r(t)|')
		axis([0 time 0 1.0])
	subplot(2,1,2)
	plot(t,abs(corrfunt),'b-.')
		title('transmitted correlation function')
		xlabel('t')
		ylabel('|C_t(t)|')
		axis([0 time 0 1.0])
     
%named- traa4.eps         
%print -djpeg -r72 'traa4.jpeg'; %-------------------------  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% calculating and plotting the spectrums.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 	for mn=1:niter
 		signa1(mn)=0.0; % S.F.T for the numeric c(t)
 		signa1(mn)=trapz(corrfun.*exp(-c*w(mn)*t));
		spectrum(mn)=dt1* 2*real(signa1(mn)-signa1(1)); 
	end
%
	figure(5)
	%subplot(2,2,1)
	plot(w,spectrum,'g-')
		title('the numeric spectrum of the whole WP')
		axis ('on')
		xlabel('energy (w)')
		ylabel('sigma(w)')
		axis ([0 60 0 1.0])
% named- traa5.eps
        
%
 	for mn=1:niter
 		sigma1r(mn)=0.0; % S.F.T for the reflected numeric c(t)
 		sigma1t(mn)=0.0; % S.F.T for the transmitted numeric c(t)
 		sigma1r(mn)=trapz(corrfunr.*exp(-c*w(mn)*t));
 	       	sigma1t(mn)=trapz(corrfunt.*exp(-c*w(mn)*t));
		
  		spectrumr(mn)=dt1*2*real(sigma1r(mn)-sigma1r(1));
		spectrumt(mn)=dt1*2*real(sigma1t(mn)-sigma1t(1));	
	end
%
	figure(6)
	%subplot(2,2,2)
	plot(w,spectrum,'g-',w,spectrumt+spectrumr,'c--', ...
             w,spectrumr,'r:',w,spectrumt,'b-.');
		title('sum of reflected and transmitted spectrum')
		xlabel('energy (w)')
		ylabel('sigma(w)')
		axis ([0 60 0 1.0])
		axis ('on')
	
% named- traa7.eps
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculating transmition and reflection coefficients
%     in three alternate ways:
%   
%  1) via the spectrum ratio
%  2) the analytical coefficient
%  3) the ratio of the k-representation WP for the reflected coefficient and
%     the ratio of the k-representation WP multiplied 
%                    by k/k` for the transmitted coefficient.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%	
%
% 2) calculate analytical Eckart barrier reflection (rho) 
%       and transmission coefficients (tau).
%
	cc=1./8.*a.*a;
	alpha1=0.5.*sqrt(w./cc);
	beta=0.5.*sqrt((w-v1)./cc);
	delta=0.5.*sqrt((v0-cc)./cc);
	rho1=cosh(2.*pi.*(alpha1-beta))+cosh(2.*pi.*delta);
	rho2=cosh(2.*pi.*(alpha1+beta))+cosh(2.*pi.*delta);
	rho=rho1./rho2;
	tau=1-rho;
%
% 3) creating vectors of k=sqrt(2*m*E/(hbar*hbar))  kk
%              and k`=sqrt(2*m*(E-v1)/(hbar*hbar))  kkt
%                   from even spaced energy vector
%
        
        energy=ak.^2/2; energy=energy(energy>v1 & energy<wmax);

        kk=sqrt(2*energy);
        mkk=-kk;
        kkt=sqrt(2*(energy-v1));

%		
% finding place of psi,psir,psit in k,k` vectors
%
	psift=fft(psinit);
	refft=fft(psir);
	transft=fft(psit);
		
		
	psikk=cspline(ak,psift,kk).*dx./sqrt(2*pi);
	refkkm=cspline(ak,refft,mkk).*dx./sqrt(2*pi);
	trakkt=cspline(ak,transft,kkt).*dx./sqrt(2*pi);

%
% Calculating reflection/transmission coefficients 
%	
        reflect=abs(refkkm).^2./abs(psikk).^2;
        transmit=(kk./kkt).*abs(trakkt).^2./ abs(psikk).^2; 
        
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%print -djpeg traa5.jpeg;   
        
        
%	
% plotting the coefficients calculated in the three different ways:
%
%  1) via the spectrum ratio
%  2) the analytical coefficient
%  3) the ratio of the k-representation WP for the reflected coefficient and
%     the ratio of the k-representation WP multiplied 
%                    by k/k` for the transmitted coefficient.
%
	figure(7)
	subplot(2,1,1)
	plot(w,spectrumr./spectrum,'g-',w,rho,'k:',energy,reflect,'c--')
		legend('spectrum coeff.','analytic coeff.','k-rep. coeff.',2)
		title('Reflection coefficient')
		xlabel('energy (w)')
		ylabel('reflection coeff.')
		axis ([10 40 0 1])
		axis ('on')
	subplot(2,1,2)
	plot(w,spectrumt./spectrum,'g-',w,tau,'k:',energy,transmit,'c--')
		legend('spectrum coeff.','analytic coeff.','k-rep. coeff.',2)
		title('Transmission coefficient')
		xlabel('energy (w)')
		ylabel('transmission coeff.')
		axis([10 40 0 1])
		axis ('on')
%		
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%print -djpeg traa6.jpeg;