
% A program to calculate the evolution of a wave packet 

%      hitting  the Eckart barrier using the split operator propagator

% Parameters of the calculation


xdim=250;                   % range of x dimension 
nop=2^10;                   % no. of points  
xrange=nop/xdim;
x0=70;                      % the starting point (50) 
dt=0.002;                   % time step in each iteration 
time=6;                     %12 the time for the whole process (45)
nstep=time/dt;              % no. of time steps per iteration 
tmax=time;                  % maximum time of propagation
dw=pi/(2*tmax);             % see NUMERICAL RECIPES about fft
wmax=pi/dt;
c=-i;                       % constant complex -i
t=0;
k0a=7;
k0b=sqrt(29);               % initial momentum
a=1;	  	            % Eckart barrier parameter
v0=75.0;		    % barrier height (initially I chose 100)
v1=10;		            % potential on right side of barrier 

                            % (initially 0.3 for assymmetric,0.0 for symmetric)
q0=70;		            % barrier coordinate
dx=xdim/nop;
alpha=0.2;         	    % initial Gaussian width (0.2)  
enstep=0.10;                % Numerical energy step for creating energy vector 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



	x=[1-nop/2:nop/2]/xrange;      % x range vector
	x=x+70*ones(1,nop);

    
% calculate k vector and vectors of kinetic and Eckart potential energy
	dk=2*pi/xdim;
        ak=[-nop/2:nop/2-1]*dk;
	kine=fftshift(ak.^2/2);
	ekine=exp(c*dt.*kine);
	ekinei=exp(-1*c*dt.*kine);
	xi=-exp((x-q0)/a);
	pote_1=zeros(1,nop);
	pote_2=10*ones(1,nop);
	pote=(-v1.*xi./(1-xi)-v0.*xi./((1-xi).*(1-xi)));
	epote_1=exp(c*dt.*pote_1);
	epote_1i=exp(-1*c*dt.*pote_1);
	epote_2=exp(c*dt.*pote_2);
	epote_2i=exp(-1*c*dt.*pote_2);
	epote=exp(c*dt.*pote);
	epotei=exp(-1*c*dt.*pote);
	epote2=exp(.5*c*dt.*pote);
	epotei2=exp(-.5*c*dt.*pote);




%	Initial Gaussian

	psinita=exp(-alpha.*(x-x0).*(x-x0)-c.*k0a.*(x-x0));
	psinitb=exp(-alpha.*(x-x0).*(x-x0)-c.*k0b.*(x-x0));

%

% Normalize psinit

%


	psiasqr=conj(psinita).*psinita;
	psibsqr=conj(psinitb).*psinitb;
	norma=(psiasqr*ones(nop,1))*dx;
	normb=(psibsqr*ones(nop,1))*dx;
	psinita=psinita/sqrt(norma);	
	psinitb=psinitb/sqrt(normb);
	phi_a=psinita;
	phi_am=conj(phi_a);
	phi_b=psinitb;

%
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   	
% creating phi_alpha- & phi_beta+ at the region of the transition state x=70,

% by operating moller operators.

% first operate exp(+-iHot),wavepackets outgoing of transition state region.
    
   in = 1;   %--------------------------------------------
	for istep=1:nstep
		phi_a=fft(phi_a);
    	phi_b=fft(phi_b);
		phi_am=fft(phi_am);
		phi_a=phi_a.*ekinei;
		phi_am=phi_am.*ekine;
		phi_b=phi_b.*ekine;
		phi_a=ifft(phi_a);
		phi_am=ifft(phi_am);
		phi_b=ifft(phi_b);
		phi_a=phi_a.*epote_1i;
		phi_am=phi_am.*epote_1;
		phi_b=phi_b.*epote_2;

        h1=figure(1);%------------------handle 4 the movie
        set(h1,'position',[0, 50,449,400])%----------------------------
		if mod(istep,50)==0
			drawnow
			xtick=[-12:3:12];
			ytick=[0:1]*.2;
         
			subplot(3,1,1);
			plot(x,pote./25,'r',x,abs(phi_a),'g-',x,abs(phi_b),'b:')
			title('Outwards propagation of \psi under e^{iH_0t}')
			xlabel('x')
			ylabel('|\phi_a & \phi_b|')
			axis([-70 220 0 1.0])
            
          
			subplot(3,1,2);
			plot(ak,abs(fftshift(fft(phi_a)).*dx./sqrt(2*pi)))
			%title('k propagation of \phi_a')
			xlabel('k')
			ylabel('|\phi_a|')
            
			subplot(3,1,3);
			plot(ak,abs(fftshift(fft(phi_b)).*dx./sqrt(2*pi)))
			%title('k propagation of \phi_b')
			xlabel('k')
			ylabel('|\phi_b|')

            
            M1(in) = getframe(h1);
            in = in + 1;

		end

	end
    
   % movie2avi(f1,'trab1.avi','compression','cinepak'); %--------------------
    

% Operate exp(-+iHt), wave packets incoming to transition region.

     in = 1;%----------------------------------------------
     h2=figure(2)  %------------------handle 4 the movie
     set(h2,'position',[50, 50,449,400])%----------------------------  
    
     for istep=1:nstep
         
	phi_a=fft(phi_a);
	phi_am=fft(phi_am);
	phi_b=fft(phi_b);
	phi_a=phi_a.*ekine;
	phi_am=phi_am.*ekinei;
	phi_b=phi_b.*ekinei;
	phi_a=ifft(phi_a);
	phi_am=ifft(phi_am);
	phi_b=ifft(phi_b);
	phi_a=phi_a.*epote;
	phi_am=phi_am.*epotei;
	phi_b=phi_b.*epotei;


     	if mod(istep,50)==0
		drawnow
		xtick=[-12:3:12];
		ytick=[0:1]*.2;
        
		subplot(3,1,1);
        plot(x,pote./25,'r',x,abs(phi_a),'g-',x,abs(phi_b),'b:')
		title('Inwards propagation of \psi under e^{\pm iHt}')
		%xlabel('x')-----------------------------------------------
		ylabel('|\phi_a & \phi_b|')
		axis([-70 220 0 1.0])
        
        
        subplot(3,1,2);
		plot(ak,abs(fftshift(fft(phi_a)).*dx./sqrt(2*pi)))
		%title('k propagation of phi_a')
		xlabel('k')
		ylabel('|\phi_a|')
        
        
		subplot(3,1,3);
		plot(ak,abs(fftshift(fft(phi_b)).*dx./sqrt(2*pi)))
		%title('k propagation of phi_b')
		xlabel('k')
		ylabel('|\phi_b|')


        M2(in) = getframe(h2);%--------------------------------
        in = in+1;%------------------------------------------
     	end

     end
     
   
     %movie2avi(M2,'trab2.avi','compression','cinepak');
     
     %movie2avi(f2,'trab2.avi','compression','cinepak');%%----------------------------------------------- saving movie in file
     %movie2avi(f21,'trab21.avi','compression','cinepak');%%%--------------------------------------------
     %movie2avi(f22,'trab22.avi','compression','cinepak');%%%--------------------------------------------
 
      phi_ap=phi_a;
      phi_api=phi_a;
      phi_bm=phi_b;
      phi_bmi=phi_bm;
%     phi_am=conj(psinita);
      phi_bp=conj(phi_b);
      phi_bpi=conj(phi_b);




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Calculating the correlation functions:

%  => C_ab(t)=<phi_a-|exp(-iHt)|phi_b+>

%  => C_ba(t)=<phi_b-|exp(-iHt)|phi_a+>

%  => C_aa(t)=<phi_a-|exp(-iHt)|phi_a+>

%  => C_bb(t)=<phi_b-|exp(-iHt)|phi_b+>

%  -inf <  t  < +inf  => we have to propagate promote the w.p. 

%          both forwards and backwards.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



	tzero=nstep+1;
	j=tzero;
	c_ab(j)=conj(phi_am).*phi_bp*ones(nop,1)*dx;
	c_ba(j)=conj(phi_bm).*phi_ap*ones(nop,1)*dx;
	c_aa(j)=conj(phi_am).*phi_ap*ones(nop,1)*dx;
	c_bb(j)=conj(phi_bm).*phi_bp*ones(nop,1)*dx;
	j=1;

	

	for istep=1:nstep

		phi_ap=fft(phi_ap);
		phi_api=fft(phi_api);
		phi_bp=fft(phi_bp);
		phi_bpi=fft(phi_bpi);
		phi_ap=phi_ap.*ekine;
		phi_api=phi_api.*ekinei;
		phi_bp=phi_bp.*ekine;
		phi_bpi=phi_bpi.*ekinei;
		phi_ap=ifft(phi_ap);
		phi_api=ifft(phi_api);
		phi_bp=ifft(phi_bp);
		phi_bpi=ifft(phi_bpi);
		phi_ap=phi_ap.*epote;
		phi_api=phi_api.*epotei;
		phi_bp=phi_bp.*epote;
		phi_bpi=phi_bpi.*epotei;
		c_ab(tzero+istep)=trapz(conj(phi_am).*phi_bp)*dx;
		c_ab(tzero-istep)=trapz(conj(phi_am).*phi_bpi)*dx;
		c_ba(tzero+istep)=trapz(conj(phi_bm).*phi_ap)*dx;
		c_ba(tzero-istep)=trapz(conj(phi_bm).*phi_api)*dx;
		c_aa(tzero+istep)=trapz(conj(phi_am).*phi_ap)*dx;
		c_aa(tzero-istep)=trapz(conj(phi_am).*phi_api)*dx;
		c_bb(tzero+istep)=trapz(conj(phi_bm).*phi_bp)*dx;
		c_bb(tzero-istep)=trapz(conj(phi_bm).*phi_bpi)*dx;

	end



     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
	figure(3)
	subplot(4,1,1);
	plot([-nstep:nstep]*dt,real(c_aa),'b');
    hold on;
	plot([-nstep:nstep]*dt,imag(c_aa),'r');
	xlabel('t');
	ylabel('Caa (Re,Im)');
    title( 'Correlation functions');



	subplot(4,1,2);
	plot([-nstep:nstep]*dt,real(c_bb),'b');
    hold on;
	plot([-nstep:nstep]*dt,imag(c_bb),'r');
	xlabel('t');
	ylabel('Cbb (Re, Im)');



	subplot(4,1,3);
	plot([-nstep:nstep]*dt,real(c_ab),'b');
    hold on;
	plot([-nstep:nstep]*dt,imag(c_ab),'r');
	xlabel('t');
	ylabel('Cab (Re, Im)');



	subplot(4,1,4);
	plot([-nstep:nstep]*dt,real(c_ba),'b');
	hold on;
	plot([-nstep:nstep]*dt,imag(c_ba),'r');
	xlabel('t');
	ylabel('Cba(Re, Im)');



	
  %print -djpeg -r60 'trab3.jpeg';  %-----------------------------------------------
   
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%         	 	

% Calculate the S matrix elements S_ba(E)=F.T.(C_ba)/(conj(eta_b)*eta_a)

%

% Calculate eta_a & eta_b

% creating vectors of k=sqrt(2*m*E/(hbar*hbar)) 

% and k`=sqrt(2*m*(E-v1)/(hbar*hbar))


% from even spaced energy vector

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% Calculating Eta_a & Eta_b by using F.F.T.

		figure(4)
        eta=fft(psinita).*dx/sqrt(2*pi);  eta=fftshift(eta);
        etb=fft(psinitb).*dx/sqrt(2*pi);  etb=fftshift(etb);
	 	energy=[v1+enstep:enstep:wmax/8];



%energy is the total energy i.e. internal energy + external energy.

		kay1=sqrt(2.*energy);
		kay2=sqrt(2.*(energy-v1));
		eta_a=cspline(ak, eta, kay1);
		eta_b=cspline(ak, etb, kay2);



		subplot(2,1,1);
		plot(energy,abs(eta_a));
		xlabel('E');
		ylabel('|\eta_a(E)|');



		subplot(2,1,2);
		plot(energy,abs(eta_b));
		xlabel('E');
		ylabel('|\eta_b(E)|');


		eta_a=eta_a./sqrt(abs(2*pi*kay1));
		eta_b=eta_b./sqrt(abs(2*pi*kay2));

       %print -djpeg -r60 'trab4.jpeg';  %-----------------------------------------------
    
       
        figure(6)
		subplot(2,1,1);      
		plot(energy,abs(eta_a));
		title('\eta_a & \eta_b');
		ylabel('\eta_a');
		xlabel('E');
        axis([0 50 0 0.4]); %----------------------
        

		subplot(2,1,2);     
		plot(energy,abs(eta_b));
		ylabel('\eta_b');
		xlabel('E');
        axis([0 50 0 0.4]); %----------------------
		eta_aa=conj(eta_a).*eta_a;
		eta_bb=conj(eta_b).*eta_b;
		eta_ab=conj(eta_a).*eta_b;
		eta_ba=conj(eta_ab);



% Calculating the S.F.T. of the correlation functions.

		sqpi=1/sqrt(2*pi);
		t=[-tmax:dt:tmax];
        lent=length(t);
		w=[0.0:dw/2:pi/dt];
		f=2*pi*[0:(length(t)-1)]/(dt * length(t));
		s_aa=ifft(c_aa)*lent*dt/sqrt(2*pi);
		s_bb=ifft(c_bb)*lent*dt/sqrt(2*pi); 
		s_ab=ifft(c_ab)*lent*dt/sqrt(2*pi); 
		s_ba=ifft(c_ba)*lent*dt/sqrt(2*pi); 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        

                saa=cspline(f,s_aa,energy);
                sbb=cspline(f,s_bb,energy);
                sab=cspline(f,s_ab,energy);
                sba=cspline(f,s_ba,energy);





	%	subplot(4,1,3);
 %       
	%	plot(energy,abs(saa));
%
%		ylabel('FT(Caa)');
%
%		xlabel('E');

	%	subplot(4,1,4);
        
		%plot(energy,abs(sbb));

		%ylabel('FT(Cbb)');

		%xlabel('E');
       
        %print -djpeg -r60 'trab6.jpeg';  %-----------------------------------------------
        
     
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% select region where the eta functions are > than epsilon=1e-7;

%   thus eliminating  numerical divergances.

%


           select=eta_aa.*eta_bb.*abs(eta_ab).^2>1e-14;
           energy=energy(select);
           eta_aa=eta_aa(select);
           eta_bb=eta_bb(select);
           eta_ba=eta_ba(select);
           eta_ab=eta_ab(select);
           saa=saa(select);
           sbb=sbb(select);
           sab=sab(select);
           sba=sba(select);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



	 figure(7)
		subplot(4,1,1);
		%plot(energy,abs(saa),'b-',energy,abs(eta_aa),'r:');
        plot(energy,abs(saa),'b-');    
		%title('F.T(C(t)) & \eta');
        title('F.T(C(t))'); 
		%ylabel('S_{aa} & \eta_{aa}');
        ylabel('C_{aa}');

		subplot(4,1,2);
		%plot(energy,abs(sbb),'b-',energy,abs(eta_bb),'r:');
        plot(energy,abs(sbb),'b-');
	    %ylabel('S_{bb} & \eta_{bb}');
        ylabel('C_{bb}');

		subplot(4,1,3);
		%plot(energy,abs(sab),'b-',energy,abs(eta_ab),'r:');
        plot(energy,abs(sab),'b-'); 
	    %ylabel('S_{ab} & \eta_{ab}');
        ylabel('C_{ab}');

		subplot(4,1,4);
		%plot(energy,abs(sba),'b-',energy,abs(eta_ba),'r:');
        plot(energy,abs(sba),'b-');
	    %ylabel('S_{ba} & \eta_{ba}');
        ylabel('C_{ba}');
        s_aa=saa;s_ab=sab;
		saa=saa./eta_aa/(2*pi);
		sbb=sbb./eta_bb/(2*pi);
		sab=sab./eta_ab/(2*pi);
		sba=sba./eta_ba/(2*pi);


      %print -djpeg -r60 'trab7.jpeg';  %-----------------------------------------------
       
        

	 figure(8)
		subplot(4,1,1);
		plot(energy,saa.*conj(saa)/(2*pi));
		title('Reflection & Transmision coefficients');
		ylabel('|Saa|^2');

		subplot(4,1,2);
		plot(energy,sbb.*conj(sbb)/(2*pi));
		ylabel('|Sbb|^2');
		subplot(4,1,3);
		plot(energy,sab.*conj(sab)/(2*pi));
		ylabel('|Sab|');
	
        subplot(4,1,4);
        plot(energy,sba.*conj(sba)/(2*pi));
		ylabel('|Sba|');
		xlabel('E');


       %print -djpeg -r60 'trab8.jpeg';  %-----------------------------------------------
       
        

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	 figure(9)

             subplot(321); plot(t,real(c_aa));
	           xlabel('t');
	           ylabel('Re(Caa)');
               
             subplot(322); plot(t,real(c_ab));
                   xlabel('t');
	           ylabel('Re(Cab)');



             subplot(323); plot(energy,abs(s_aa));
                   ylabel('FT(Caa)');
		   xlabel('E');

             subplot(324); plot(energy,abs(s_ab));
		   ylabel('FT(Cab)');
		   xlabel('E');

     subplot(325); plot(energy,saa.*conj(saa)/(2*pi));
		   ylabel('|Saa|^2');
     subplot(326);plot(energy,sab.*conj(sab)/(2*pi));
		   ylabel('|Sab|^2');



    %print -djpeg -r60 'trab9.jpeg';  %-----------------------------------------------