%  LKB Photon Box: open-loop/closed-loop  simulation of a single  trajectory
%  Physical description and feedback scheme:
%  I. Dotsenko, M. Mirrahimi,M. Brune, S. Haroche, J.-M. Raimond, P. Rouchon: 
%   Quantum feedback by discrete quantum non-demolition measurements: 
%   towards on-demand generation of photon-number states. 
%   Physical Review A,80: 013805-013813, 2009 (preliminary version http://arxiv.org/abs/0905.0114 ).
%  Convergence proof:   
%  H.~Amini, M.~Mirrahimi, and P.~Rouchon.
%    Stabilization of a delayed quantum system: the Photon Box case-study. 2010.
%     Submitted ( http://arxiv.org/abs/0905.0114 )
% 
%  IHP course 
% mazyar.mirrahimi@inria.fr
% pierre.rouchon@mines-paristech.fr
% October 2010
clear all; i=sqrt(-1); rand('twister',sum(100*clock));
%% nominal parameter
Tpulse=85.e-6; % intervalle between to atom
Tcav=0.13e20; % decay time of the cavity (0.15)
nThermal = 0.05e-10; % number of thermal photons (0.05)
PulseOccupancy =1;% probability to have an atom in a pulse (0.3)
Ndelay=0; % number of atom/pulse between the cavity and the detector (4)
DetectorEfficiency=1; %  detection probability (0.8)
DetectorFault= 0; % fault detection probability (0.1)
FlagFeedback=0 ; % 0 open-loop, 1 closed-loop
%  operator and matrices
NphotonMax=9; % max number of photons
Ph0=zeros(NphotonMax+1,1);Ph0(1)=1; % 0 photon state
A=diag(sqrt([1:1:NphotonMax]),1); % annihilation operator
NN=diag([0:1:NphotonMax]); % photon  number operator
Id=eye(NphotonMax+1);
%  Kraus operators
Ngoal=3; % mid fringe photon number
PhGoal=zeros(NphotonMax+1,1);PhGoal(1+Ngoal)=1; rhoGoal=PhGoal*PhGoal';
phi=0.2; % dephasing per photon
G=diag(cos(pi/4+phi*([0:1:NphotonMax]-Ngoal-1))); % Generalized  measure operator if g
E=diag(sin(pi/4+phi*([0:1:NphotonMax]-Ngoal-1))); % Generalized measure operator if e
%% Controller setting
alphaKick=1;
radius=([0:0.1:1].^2)*alphaKick; phase=exp(sqrt(-1)*[-1:0.1:1]*pi);
Aaa=radius'*phase; Aaa=reshape(Aaa,numel(Aaa),1);
PhAaa=zeros(length(PhGoal),length(Aaa));
for ii=1:length(Aaa)
    PhAaa(:,ii)=expm(-Aaa(ii)*A'+(Aaa(ii)')*A)*PhGoal;
end
%% initial state
alphaCav=sqrt(Ngoal); PhCav=expm(alphaCav*A'-alphaCav'*A)*Ph0; % initial  cavity  state (classical with Ngoal photon in average)
rhoCav=PhCav*PhCav'; % density operator
rhoEst=rhoCav;
%% loop init and  graphics buffers
Npulse=floor(100/PulseOccupancy); % number of step
Iatom=[1:Npulse];  % sampling atom indice
Yatom=zeros(Npulse,2); % detection -1=g and 1=e
Alpha=zeros(Npulse,1); % the control
Ycavp=[]; Icavp=[]; % photon jump +1
Ycavm=[]; Icavm=[]; % photon jump -1
Photon=zeros(NphotonMax+1,Npulse); % photons population
figure(3);
fig3=bar([0:1:NphotonMax],diag(rhoCav));set(gca,'fontsize',13);
set(gca,'YLimMode','manual','Ylim',[-0.01 1.01]);
set(gca,'XLimMode','manual','Xlim',[-0.5 NphotonMax+0.5]);
xlabel('Photon Number'),Ylabel('Population'); 
if (FlagFeedback <0.5 )
    title(['Open-Loop (step number ',num2str(0),')'],'fontsize',14);
else
    title(['Closed-Loop N_g_o_a_l=',num2str(Ngoal),'    (step number ',num2str(0),')'],'fontsize',14);
end
pause;
for ia=1:1:Npulse;
    Photon(:,ia)=(diag(rhoCav));
    % prediction over Ndelay step  (stochastic version of Smith predictor)
    rhoPG=rhoEst;rhoPE=rhoEst;
    for ip=max(1,ia-Ndelay):max(1,ia-1);
        UU=expm( Alpha(ip)*(A')-(Alpha(ip)')*A );
        rhoPG=UU*rhoPG*(UU');rhoPE=UU*rhoPE*(UU');
        if (ip==1);
            rhoPG=G*rhoPG*(G'); rhoPE=E*rhoPE*(E');
        else
            rhoPG=G*rhoPG*(G')+E*rhoPG*(E');
            rhoPE=G*rhoPE*(G')+E*rhoPE*(E');
        end
    end
    rhoP=rhoPG+rhoPE; rhoP=rhoP/trace(rhoP);
    alpha=trace((A*rhoGoal-rhoGoal*A)*rhoP)/(2*Ngoal+1);
    if ( real(trace(rhoGoal*rhoP))<0.1);
        buf=0; alpha=0;
        for iaaa=1:length(Aaa)
            sss= (PhAaa(:,iaaa)'*rhoPG*PhAaa(:,iaaa))*(PhAaa(:,iaaa)'*rhoPE*PhAaa(:,iaaa));
            % sss= (PhAaa(:,iaaa)'*rhoP*PhAaa(:,iaaa));
            if(sss>buf)
                alpha=Aaa(iaaa);
                buf=sss;
            end,
        end
    end
    %%
    Alpha(ia)=FlagFeedback*alpha;
    % injection of the  classical field (with delay)
    alpha=Alpha(max(1,ia-Ndelay));
    UU=expm(alpha*(A')-(alpha')*A);
    
    %% Cavity
    rhoCav=UU*rhoCav*UU';
    rhoCavG=G*rhoCav*G'; Pg=(trace(rhoCavG));rhoCavG=rhoCavG/Pg;
    rhoCavE=E*rhoCav*E'; Pe=(trace(rhoCavE));rhoCavE=rhoCavE/Pe;
    if (rand < PulseOccupancy )
        % pulse containing an atom
        if (rand < DetectorEfficiency )
            % atom detected
            if (rand < Pg)
                % collapse  in g
                rhoCav=rhoCavG;
                if (rand < DetectorFault )
                    Yatom(ia,2)=1;% wrong detection in e
                else
                    Yatom(ia,1)=-1;% correct detection in g
                end;
            else
                % collapse in e
                rhoCav=rhoCavE;
                if (rand < DetectorFault )
                    Yatom(ia,2)=-1;% wrong detection in g
                else
                    Yatom(ia,1)=1;% correct detection in e
                end;
            end
        else
            % atom not detected, collapse on e and g
            if (rand < Pg)
                % collapse in g without detection
                rhoCav=rhoCavG;
            else
                % collapse in e without detection
                rhoCav=rhoCavE;
            end
        end
    end
    % Cavity decay + thermal photons
    Pc1=1-trace((Id-0.5*Tpulse/Tcav*NN*(1+nThermal))*rhoCav*(Id-0.5*Tpulse/Tcav*NN*(1+nThermal)));
    Pc2=1-trace((Id-0.5*Tpulse/Tcav*A*A'*nThermal)*rhoCav*(Id-0.5*Tpulse/Tcav*A*A'*nThermal));
    ss=rand;
    if ( ss < Pc1 ); %  -1 photon
        rhoCav= A*rhoCav*A'; rhoCav=rhoCav/(trace(rhoCav));
        Ycavm=[Ycavm, 0]; Icavm=[Icavm, ia];
    elseif (ss < Pc1+Pc2 ) % +1 thermal photon
        rhoCav= A'*rhoCav*A; rhoCav=rhoCav/(trace(rhoCav));
        Ycavp=[Ycavp, 0]; Icavp=[Icavp, ia];
    else % no jump
        buf=(1+nThermal)*NN+nThermal*A*A';
        rhoCav= (Id-0.5*Tpulse/Tcav*buf)*rhoCav*(Id-0.5*Tpulse/Tcav*buf);
        rhoCav=rhoCav/(trace(rhoCav));
    end
    %% quantum filter
    rhoEst=UU*rhoEst*UU';
    rhoEstG=G*rhoEst*(G'); PgEst=real(trace(rhoEstG)); rhoEstG=rhoEstG/PgEst;
    rhoEstE=E*rhoEst*(E'); PeEst=real(trace(rhoEstE)); rhoEstE=rhoEstE/PeEst;
    if (Yatom(ia) < 0)
        % detection in g
        % probability to have a wrong detection knowing a detection in g.
        Pfalse=DetectorFault*PeEst/((1-DetectorFault)*PgEst+DetectorFault*PeEst);
        rhoEst=(1-Pfalse)*rhoEstG+Pfalse*rhoEstE;
    elseif (Yatom(ia) > 0 )
        % detection in e with probability 1-DetectorFault
        % probability to have a wrong detection knowing a detection in e.
        Pfalse=DetectorFault*PgEst/((1-DetectorFault)*PeEst+DetectorFault*PgEst);
        rhoEst=(1-Pfalse)*rhoEstE+Pfalse*rhoEstG;
    else
        % no detection:
        %  probability  to have an atom  knowing no detection
        Patom=PulseOccupancy*(1-DetectorEfficiency)/(1.e-10+1-PulseOccupancy*DetectorEfficiency);
        rhoEst=(1-Patom)*rhoEst + Patom*(G*rhoEst*(G')+E*rhoEst*(E'));
    end
    buf=(1+nThermal)*NN+nThermal*A*A';
    rhoEstNo=(Id-0.5*Tpulse/Tcav*buf)*rhoEst*(Id-0.5*Tpulse/Tcav*buf);
    rhoEstNo=rhoEstNo/(trace(rhoEstNo));
    PcEst=(trace(buf*rhoEst))*Tpulse/Tcav;
    rhoEst=(1-PcEst)*rhoEstNo+(Tpulse/Tcav)*((1+nThermal)*A*rhoEst*A'+nThermal*A'*rhoEst*A);
    
    %% graphics
    set(fig3,'Ydata',Photon(:,ia)');
    if (FlagFeedback <0.5 )
        title(['Open-Loop (step number ',num2str(ia),')'],'fontsize',14);
    else
        title(['Closed-Loop N_g_o_a_l=',num2str(Ngoal),'    (step number ',num2str(ia),')'],'fontsize',14);
    end
    pause(0.2);
end
%% visualisation
figure(1);clf;
subplot(311);set(gca,'fontsize',12);
hb=bar([1:Npulse],Yatom,'stack');
set(gca,'Xlim',[1 Npulse]);
title(['Detection (-1=g, 0=none,  +1=e). Total of clicks: ',num2str(sum(abs(Yatom(:,1))))...
    ,' (',num2str(sum(abs(Yatom(:,2)))),' wrong in red)']...
    ,'fontsize',12);

subplot(312);set(gca,'fontsize',12);
bar([1:Npulse],real(Alpha));
set(gca,'Xlim',[1 Npulse]);% set(gca,'Ylim',[-alphaLim alphaLim])
title('Control Re(\alpha)','fontsize',12);

subplot(313);set(gca,'fontsize',12);
Ph_inf=sum(Photon([1:Ngoal],:),1);
Ph_sup=sum(Photon([Ngoal+2:NphotonMax+1],:),1);
Ph_goal=Photon(Ngoal+1,:);
plot([1:Npulse],Ph_inf,'b-',...
    [1:Npulse],Ph_goal,'k-',...
    [1:Npulse],Ph_sup,'r-',...
    Icavp,Ycavp-0.05,'r*',Icavm,Ycavm-0.05,'b*','linewidth',2);
set(gca,'Xlim',[1 Npulse]);set(gca,'Ylim',[-0.05 1.05]);
ylabel('proba.');
title(['Photon: blue'...
    ,'(0-',num2str(Ngoal-1),'); black(' ...
    ,num2str(Ngoal),'); red(' ...
    ,num2str(Ngoal+1),'-',num2str(NphotonMax),'). Jump: *red(+1); *blue(-1).'],'fontsize',12);
xlabel('Step number');


figure(2);
bar3(Photon');
set(gca,'fontsize',13);
Xlab=[]; for n=0:NphotonMax; Xlab=[Xlab; num2str(n)];end;
set(gca,'XTickLabel',Xlab);
ylabel('Step Number');
xlabel('Photon Number'),zlabel('Probability');
    if (FlagFeedback <0.5 )
        title(['Open-Loop'],'fontsize',14);
    else
        title(['Closed-Loop N_g_o_a_l=',num2str(Ngoal)],'fontsize',14);
    end


return
