% LKB photon box: 
% Optimal QND measurement of photon according to 
% Haroche-Brune-Raimond: Measuring photon numbers in a cavity 
% by atomic interferometry: optimizing the convergence procedure.
% J. Phys. II France, 1992
% https://doi.org/10.1051/jp2:1992157
%
% pierre.rouchon@minesparis.psl.eu
% density operator formulation
% Januray 2025
%
clear all;
%% operator for the quantum harmonic oscillator
i=sqrt(-1); m=4; Nmax=2^m-1; % truncation to Nmax
Vacuum=zeros(Nmax+1,1);Vacuum(1)=1;
Id=eye(Nmax+1);
N=diag((0:1:Nmax)); % photon  number operator
A=diag( sqrt((1:1:Nmax)),1); % annihilation operator
%% simulation parameter
Niter=m;  
%Psi= expm(sqrt(3)*(A-A'))*Vaccum; % initial state coherent state with 3 photons 
%rho=Psi*(Psi');
rho=Id/(Nmax+1);
%% simulation loop
Pop=zeros(Nmax+1,Niter+1);y=zeros(1,Niter);% buffer for graphics
u=zeros(Niter);v=u; f=y;Bm=y;
u(1)=0;v(1)=pi/2; uk=u(1);vk=v(1);
Pop(:,1)=diag(rho);
for el=0:m-1
    Bm(el+1)=2^el;
end
for k=1:Niter
    u(k)=uk;v(k)=vk;
    Mg=diag(cos((vk*(0:1:Nmax)+uk)));
    Me=diag(sin((vk*(0:1:Nmax)+uk)));
    rhog=Mg*rho*Mg';rhoe=Me*rho*Me';
    if (rand < trace(rhog))
        y(k)=-1; f(k)=0;
        rho=rhog;
    else
        y(k)=1; f(k)=1; 
        rho=rhoe;
    end
    rho=rho/trace(rho);
    Pop(:,k+1)=diag(rho);
    uk=-pi/2^(k+1)*sum(f.*Bm);
    vk=pi/2^(k+1);
end
%% graphics
figure(1),clf,
subplot(211);
bar([1:Niter],y);
ylabel('Measurement outcome');
axis([0 Niter+1 -1.02 1.02 ]);
ax = gca; ax.FontSize = 14;
subplot(212);
plot((0:Niter),Pop,'linewidth',4);
axis([0 Niter -0.02 1.02 ]);
ax = gca; ax.FontSize = 14;
ylabel('populations');
xlabel('step number');
ax = gca; ax.FontSize = 14;
title(['Final Fock state n=',num2str(sum(f.*Bm))],'Interpreter','latex')

%%
figure(2);clf;
%y = 0:pi/16:pi;
%z = [sin(y')/4 sin(y')/2 sin(y')];
%bar3(y,z,1,'r')
bar3((0:1:Nmax),Pop,1);
