% Coherent population trapping
% stabilisation du detuning Raman à partir du signal de fluorescence 
% pierre.rouchon@mines-paristech.fr
% le 17 septembre 2008
% simulation présentées à la conf. 
%   Quantum/Classical Control in Quantum Information: \\
%   Theory and experiments. September 13-20, 2008. Otranto, Italy
%
% d rho/dt = - i [H,rho] 
%                   + G1 ( L1 rho L1'  - L1' L1 rho/2 - rho L1' L1/2)
%                   + G2 ( L2 rho L2'  - L2' L2 rho/2 - rho L2' L2/2)
%  y= G1 trace(L1' L1 rho) +  G2 trace(L2' L2 rho)
% avec rho matrice densite d'un trois niveaux
%  1 et 2 niveau stable
%   3=e niveau instable: 
%        L1= |1><e| durée vie 1/G1
%        L2= |2><e| durée vie 1/G2
% 
% H = d |e><e| + D |2><2| + W1 |1><e| + cc + W2 |2><e| + hc
% excitation Raman via deux lasers, detuning Raman: D 
%     1 <-> e d'amplitude W1  detuning  transition: De 
%      2<->e d'ampltiude  W2  
%    
% Hypothèse lent/rapide: G1,G2 >> |W1|, |W2|, D, De 
% 
% modèle de contrôle, le modèle lent suivant
% 
% d rhos/dt = - i [Hs,rhos] 
%                   + G1s (L1s rhos L1s'  - L1s' L1s rhos/2 - rhos L1s' L1s/2)
%                   + G2s (L2s rhos L2s'  - L2s' L2s rhos/2 - rhos L2s' L2s/2)
%  y= G1s trace(L1s' L1s rhos) +  G2s trace(L2s' L2s rhos)
%
% Hs=D |2><2|, L1s= |1><a|, L2s=|2><a| 
% |a> = (W1 |1> + W2|2>/(|W1|^2+[W2|^2)^(1/2)
% G1s= 4G1*(|W1|^2+[W2|^2)/(G1+G2)^2, 
% G2s= 4G2*(|W1|^2+[W2|^2)/(G1+G2)^2, 
% 
% Objectif faire converger D vers 0 
% avec le contrôle d D/dt = u
% on peut utiliser W1 et W2 comme on veut 
% simulation d'une trajectoire pour un atome unique 
%  quantum Monte-Carlo simulation de Dalibard
% mise en place d'un terme integral pour ramener le detuning à zero à partir des
% temps entre deux détections de clic
% 
clear all; 
% parametre nominal   systeme
De=2; W1=1; W2=1;  
G1=10; G2=10;
% temps simu et pas discretisation 
G=G1+G2; Tfast=1/(G1+G2); 
g=4*(W1^2+W2^2)/(G1+G2); Tslow=1/g; 
ne=50; dt=Tfast; 
Tsimu=ne*800*Tslow;  nT=floor(Tsimu/dt);
% gain 
eps=sqrt(g/G);  
omega=g/eps;
K=eps*g/25;
sig2=eps*g*K/4;
Dmax=g/2;
% bruit mesure fausse detection et retard
eta_noise= 0; % eps*eps/2;
eta_detec=1; % 0.5; 
delay=0; % (pi/4)/omega;
% %  vecteurs et operateurs,
i=sqrt(-1); 
e=[0 0 1]'; g1= [1 0 0]'; g2=[0 1 0]'; 
P=e*e'; P12=g1*g1'+g2*g2'; L1=g1*e';L2=g2*e'; 
%  initial condition 
t=0;
%psi=(g1+g2)/sqrt(2);  % etat brillant
psi=(g1-g2)/sqrt(2); % etat noir 
y=0;
D=g/2; 
% buffer stockage 
t_b=zeros(1,floor(nT/ne)+1); y_b=t_b; D_b=t_b;
iie=0; xi=0; Tjump=0; 
for ii=1:nT
    ss=eps*cos(omega*t); 
    W1t=W1-i*ss*W2;
    W2t=i*ss*W1+W2;
    % 
    pL1=dt*G1*psi'*L1'*L1*psi; 
    pL2=dt*G2*psi'*L2'*L2*psi; 
    pnoise=dt*eta_noise/Tslow; 
    p0=1-pnoise-pL1-pL2;
    pn=p0+pnoise; 
    p1=pn+pL1;
    r=rand; 
    if (r < p0)
        % no jump no detector click 
        H=W1t*g1*e'+W2t*g2*e';H=H+H'; H=H-(De+D/2)*e*e'+D*g2*g2';
        He= H/2 - i* (G1*L1'*L1+G2*L2'*L2)/2;  
        psi = (eye(3) - i*dt *He)*psi; psi=psi/sqrt(psi'*psi);
        dy=0; 
    elseif (r<pn)
        % no jump but a wrong detector click 
        H=W1t*g1*e'+W2t*g2*e';H=H+H'; H=H-(De+D/2)*e*e'+D*g2*g2';
        He= H/2 - i* (G1*L1'*L1+G2*L2'*L2)/2;  
        psi = (eye(3) - i*dt *He)*psi; psi=psi/sqrt(psi'*psi);
        dy=1; 
    elseif (r<p1)
        % jump  L1
         psi=L1*psi;psi=psi/sqrt(psi'*psi);
        if (rand <=eta_detec) % efficacité  du détecteur 
             dy=1;
        else
            dy=0;
        end
    else
        % jump L2
        psi=L2*psi;psi=psi/sqrt(psi'*psi);
        if (rand <=eta_detec) % efficacité  du détecteur 
            dy=1;
        else
            dy=0;
        end
    end
    y=y+dy;
    D=D-K*cos(omega*(t-delay))*2*W1*W2/(W1^2+W2^2)*dy;
    t=t+dt;    
    % save
    iie=iie+1;
    if ( ne*floor(iie/ne) == iie)
        iie=0;
        iii=floor(ii/ne)+1;
        t_b(iii)=t;
        y_b(iii)=y;
        D_b(iii)=D;
     end;
end;        
Rms=ones(1,length(t_b))*sqrt(sig2);
figure(1),clf,
  subplot(211)
    plot(omega*t_b/(2*pi),y_b,'linewidth',2); grid on; set(gca,'fontsize',14);
    legend('Click number'); 
    xlabel('\omega t / 2 \pi'); 
  subplot(212)
    plot(omega*t_b/(2*pi),D_b/sqrt(sig2),'linewidth',2); grid on; set(gca,'fontsize',14)
    legend('Normalized de-tuning \Delta/\sigma_K'); 
    xlabel('\omega t / 2 \pi'); set(gca,'Ylim',[-10,40]);
    
    
return

    
    
    