function L=L3lvl(H,gam1,gam2)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   L3lvl
%   Utility function for optimize.m and propagate,m
%   function creates Liouville space operator for the 3 level lambda system
%   The Liouvillian L is a linear map of density matrices to density matrices
%   operating in the liouville space of density matrices
%   rho_{ij}=L_{ijkl}rho_{kl}
%
%   The Hamiltionian part:  L_{ijkl}=H_{ij}delta_{kl}-delta_{ij}H_{ji}
%   The Dissipative part is more complex
%
%   Input parameters:
%    H- Hamiltonian operator in Hilbert space.
%    gam1,gam2- decay rates to each of the ground states.
%   Output parameters
%    L- Liouvillian.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



     %     gam1=.05; gam2=.05;
     N=3;
    
     %Ld=sparse(N^2,N^2);
     Ld=zeros(N^2);

     % Superoperator for Ld(rho)
     for idx=1:N;
        for jdx=3:N;
           x=idx+N*(jdx-1); y=jdx+N*(idx-1);
           Ld(x,x)= Ld(x,x)-.5*(gam1+gam2);
           Ld(y,y)= Ld(y,y)-.5*(gam1+gam2);
        end;
     end;      

     x=1:N+1:N^2; 

     Ld(x(1),x(3))=gam1; Ld(x(2),x(3))=gam2;

     % L=sparse(N^2,N^2); 
     L=zeros(N^2);
     
     % Superoperator for [H,rho]
     for idx=0:N-1;   
         L(N*idx+[1:N],N*idx+[1:N])=L(N*idx+[1:N],N*idx+[1:N])+H;
         L(idx+[1:N:N^2],idx+[1:N:N^2])=L(idx+[1:N:N^2],idx+[1:N:N^2])-H.';
     end;

     L= (-i*L + Ld);
return;
