% LKB photon box: QND measurement of photon
% measurement-based Feedback stabilization of a target  photon-number state
%  1- detection efficiency
% February 2016
% pierre.rouchon@mines-paristech.fr
% alain.sarlette@inria.fr
clear all;
%% operator for the quantum harmonic oscillator
i=sqrt(-1); Nmax=7; % truncation to Nmax
Ket0=zeros(Nmax+1,1);Ket0(1)=1; % vacuum state
Id=eye(Nmax+1);
N=diag((0:1:Nmax)); % photon  number operator
A=diag( sqrt((1:1:Nmax)),1); % annihilation operator
%% Lyapunov-function operator
ntarget= 3; % target  photon-number state
Lyap=diag(((0:1:Nmax)-ntarget).^2);
%%  measurement operators depending on the control u=0,-1,+1
% u=0 subscript _0 
ph0=pi/7; phR=0; 
Mg_0=diag(cos((ph0*(0:1:Nmax)+phR)/2)); 
Me_0=diag(sin((ph0*(0:1:Nmax)+phR)/2)); 
% u=-1 subscript _m
th0_m=2*pi/sqrt(ntarget-1);
Mg_m=diag(cos(th0_m*sqrt((0:1:Nmax))/2)); 
Me_m=-diag(sin(th0_m*sqrt((1:1:Nmax))/2),1); 
% u=+1 subscript _p
th0_p=2*pi/sqrt(ntarget+2);
Mg_p=( diag(sin(th0_p*sqrt((1:1:Nmax))/2),1) )'; 
Me_p=diag(cos(th0_p*sqrt((1:1:Nmax+1))/2)); 
%% measurement imperfections
eta=1/3; % detection efficiency
%% simulation parameter
Niter=300*1.5; 
Psi=expm(sqrt(3)*(A'-A))*Ket0; % initial state coherent state with an average of 3 photons
rho=Psi*(Psi'); % initial state (no information) 
%% simulation loop
Pop=zeros(Nmax+1,Niter);y=zeros(1,Niter);u=zeros(1,Niter); % buffer for graphics
for ii=1:Niter
    Pop(:,ii)=diag(rho);
    % quantum-state feedback
    V_0=trace(Lyap*(Mg_0*rho*Mg_0'+Me_0*rho*Me_0'));
    V_m=trace(Lyap*(Mg_m*rho*Mg_m'+Me_m*rho*Me_m'));
    V_p=trace(Lyap*(Mg_p*rho*Mg_p'+Me_p*rho*Me_p'));
    Vmin=V_0;Mg=Mg_0;Me=Me_0;u(ii)=0;
    if (V_p < Vmin);Vmin=V_p;Mg=Mg_p;Me=Me_p;u(ii)=1;end;
    if (V_m < Vmin);Vmin=V_m;Mg=Mg_m;Me=Me_m;u(ii)=-1;end;
    % measurement with Mg and Me
    rhog=Mg*rho*Mg';rhoe=Me*rho*Me';
    if(rand < eta)
        if (rand < trace(rhog))
            y(ii)=-1;
            rho=rhog/trace(rhog);
        else
            y(ii)=1;
            rho=rhoe/trace(rhoe);
        end
    else
        y(ii)=0;
        rho=rhog+rhoe;
    end
    rho=rho/trace(rho);
end
%% graphics
figure(1),clf,
subplot(211);
bar([1:Niter],u);
ylabel('Control  u ');
axis([0 Niter+1 -1.02 1.02 ]);
subplot(212);
plot([1:Niter],Pop,'linewidth',2);
legend('0','1','2','3','4','5','6','7');
axis([1 Niter -0.02 1.02 ]);
ylabel('populations');
xlabel('step number');
title(['Target of ',num2str(ntarget),' photon(s)']); 



