% this program calculates numerically the solution of the generalized
% Rabi problem, i.e. two potential energy curves coupled by a strong
% field of a single frequency coupled to the transition dipole moment.
  % For simplicity, the program uses harmonic potentials.
% setting the parameter xd=0 makes the displacement between the excited
  % and ground states equal to 0; as a result the only transition possible
% is 0 <-> 0, the problem reduces to that of a two-level system, and the
% Rabi solution is recovered.
 
% parameters of the calculation

clear

nop=128;                    % no of points /2
dx=sqrt(pi/nop);
xdim=nop*dx;
dpl=4;	        	    % pulse duration
dt=.002;                   % time step in each iteration
nstep=10;                   % no of time steps per iteration
niter=dpl/(dt*nstep);       % no of iteration
alpha_g = 2.0 ;
%alpha_e = 4.0 ;
alpha_e = 2.0 ;
xd=-3.0;                    % the shift between V_g and V_e
xd=0;
%x0=5.0;                     % initial shift of the wp, so it can oscillate
x0=0.0;
e0=4.0;                     % field amplitude
v_g = 10 ;
v_e = 50 ;
% omega=40.3001;              % field frequency
omega=v_e-v_g;          % zero detuning
mue=1.0;
t1=0.0;
c=-i;                       % constant copmlex -i
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%             
        x=[-1*nop:nop-1]*dx;  

        n=length(x);

% calculate k vector and vectors of kinetic and potential energy

        ak1=[0:nop-1]*pi/xdim;
        ak2=[-1*nop:-1]*pi/xdim;
        ak=[ak1 ak2];

% m=1, hbar=1, w=sqrt(2.0)

        kine=1.*ak.*ak/2;
        ekine=exp(c*dt.*kine);

%   ground state potential
        potg= v_g + 0.5 * alpha_g *x.*x;

%   excited state potential
        pote= v_e + 0.5 * alpha_e * (x-xd).*(x-xd);

        dsum=pote+potg;
        difeg=pote-potg;
        difeg2=difeg.*difeg;
        edif=exp(0.5*c*dt.*dsum);
        
        

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   initial wavefunctions: for coherent state alpha=mw/(2*hbar)
%    on ground state
        psinitg1=exp(-sqrt(2.0)/2.0 * (x+x0).*(x+x0));
%    on excited state
        psinite=0.*psinitg1;
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
%  normalize psint;
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

psisqr=psinitg1.*conj(psinitg1);
normg=(psisqr*ones(n,1));
psinitg=psinitg1/sqrt(normg*dx);

	psig=psinitg;
	psie=psinite;

        pnorg(1:n)=psinitg.*conj(psig);
	corr_g(1)=(pnorg(1:n)*ones(n,1))*dx  % =1
 
        pnore(1:n)=psinite.*conj(psie);
	corr_e(1)=(pnore(1:n)*ones(n,1))*dx  % =0
    
        corr_tot(1)=corr_g(1) + corr_e(1);

        laser(1)=-mue*e0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%	

        figure(1);clf;
        plot(x,potg,'k-',x,pote,'k')
        title('Ground and excited  potentials')
        axis([-15 15 0 300])
        hold on
        drawnow
      
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%	

	psig1=psig;
	psie1=psie;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% calculate the propagation for ntime time steps
%
       h2=figure(2);  %------------------------------------
       set(h2,'position',[50 50 300 250]); %----------------
       set(h2,'color','w');  %------------------------------

       for iter=2:niter
         for istep=1:nstep

%          drabi=e0*exp(-c*omega*t1)*(sin(pi*t1/dpl))^2;
%          drabi1=drabi.*dipm;
	   drabi1=-mue*e0*cos(omega*t1);
	   dis=sqrt(difeg2+4.*drabi1.*conj(drabi1));
	   disdt=.5*dt.*dis;

           sin_overD=0.5*dt*sinc(disdt);
           V11=(cos(disdt)-c*sin_overD.*difeg);
           V12=2*c*sin_overD.*drabi1;
           V21=V12;
           V22=(cos(disdt)+c*sin_overD.*difeg);

           psig=V11.*psig1 + V12.*psie1;
           psig=edif.*psig;
           
           psie=V21.*psig1 + V22.*psie1;
           psie=edif.*psie; 

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

	   psig1=fft(psig);
       psig1=psig1.*ekine;
       psig1=ifft(psig1);
       psie1=fft(psie);
	   psie1=psie1.*ekine;
	   psie1=ifft(psie1);

	   t1=t1+dt;

         end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

       
         plot(x,abs(psig1),'b-',x,abs(psie1),'r-')
         %title('psi_g (solid), psi_e (dotted)')
         legend('\psi_g','\psi_e');  %-----------------
         axis([-7 7 0 1])
         drawnow
         
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%	   
 
%        calculate  N(t)  

         pnorg(1:n)=psig1.*conj(psig1);
	     corr_g(iter)=(pnorg(1:n)*ones(n,1))*dx;
 
         pnore(1:n)=psie1.*conj(psie1);
	     corr_e(iter)=(pnore(1:n)*ones(n,1))*dx; 

         corr_tot(iter)=corr_g(iter) + corr_e(iter);
         laser(iter)=drabi1;
         
      
       end
             
    

       t=[0.0:nstep*dt:dpl-nstep*dt];
       
       
          figure(3);    %-------------------------------------    
          plot(t,corr_g,t,corr_e,'k:',t,corr_tot,'k--')
          title('population: g(solid), e(dotted), sum(dashed)')
          axis([0 dpl 0 1.2])
        
          figure(4);   %-------------------------------------
          plot(t,laser)
          title('the laser pulse')
          axis([0 1 -4.5 4.5]) 
     
      
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      