
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%                                       %
	%     FILE NAME: control_3level_diss    %
	%            				            %
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%   
%  PROPAGATE
%   A pedagogical code for propagating the density matrix of a three level 
%   Lambda system governed by the shrodinger equation for the density matrix
%   with an additional dissipative operator Gamma:
%   (d/dt)rho=(-i/hbar)[H,rho]+Gamma(rho)
%   where rho is the density matrix H the hamiltonian and Gamma is the 
%   dissipative operator. 
%                                                             ----- |3>
%                                                              / \
%                                                             /   \
%                                                     |1>   ---   --- |2>
%
%   rho0:  Initial density matrix;
%   H0:    Energies of unperturbed system
%   mu:    Interaction cuoplings (the dipole strengths)
%   Einit: Initial guess of electric field.  
%   gam1,2:Decay rates to corresponding ground levels.
%
%  Comments:
%    a) A Lindblat form is used for the Decay operator.
%               [ gam1*rho_33       0               -0.5*gam*rho_13 ]
%         Gamma=[     0             gam2*rho_33     -0.5*gam*rho_23 ]
%               [ -0.5*gam*rho_31  -0.5*gam*rho_32  -gam*rho_33     ]
%    b) Propagation method: Brute force diagonalization:
%       Equation for density matrix is written in Liouville space as
%          rho=L *rho  where L is a N^2 x N^2 matrix.
%       The solution is then formally:
%          rho(t+dt)=exp(L dt)rho(t)
%       which is solved by diagonalizing L=U Z inv(U) 
%               where Z is diagonal.
%       We can now write:
%          rho(t+dt)=U exp(Zdt)inv(U)rho(t)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55

  hbar=1;
  E1=0; E2= 3; E3=0;

  rho0=[ 1/3 0 0;
	  0 1/3 0;
	  0 0 1/3];

  H0=diag([E1 E2 E3]);
  mu=[-1 -1];
  B=[1 0 0;
     0 0 0;
     0 0 0];

  Niter=100;
  Nt=2^11;
  dt=.02;
  lambda=1;
  b=zeros(9,Nt); rho=zeros(9,Nt);
  t=[0:Nt-1]*dt;


  % Einit=.000001*rand(Nt,1);
  % Einit=2*(exp(i*10*t)+exp(-i*10*t));
  % Einit=cos(10*t); 
  
  Einit=.2*exp(-i*3*t);  
  
  gam1=0.05; gam2=0.05; gamma=gam1+gam2;

  E=Einit;

% Initial objective:
  Goal=trace(B*rho0);

  rho(:,1)=rho0(:); 
 %%%% Forward propagation
  for idx=1:Nt-1;

   % Hamiltonian
    muE=mu*conj(E(idx));
    H1= [zeros(2,2) muE';muE 0];
    H=(H0+H1);

    L=L3lvl(H,gam1,gam2);


   % rho=L rho      
   % rho(t)=exp(Lt)rho=U^-1 exp(Zt) U rho(0)

     [U,Z]=eig(L);
     rho(:,idx+1)=U*(exp(diag(Z)*dt).*(inv(U)*rho(:,idx)));

  end;

figure;

  subplot(211);plot(t,E);
    title('Field profile');
    xlabel('Time')
    ylabel('Field intensity')

  subplot(212);plot(t,rho([1 5 9],:)');
    title('populations');
    xlabel('time');
    ylabel('Populations');
    legend('\rho_{11}','\rho_{22}','\rho_{33}')
 

 return;
