% LKB photon box: QND measurement of photon
%
% feedback scheme principle of Zhou,....,Haroche: Field locked to Fock state by quantum feedback with single photon corrections 
% 2012, Physical Review Letter, Vol. 108, p. 243602
%
% November  2023
% pierre.rouchon@minesparis.psl.eu
%
clear all;
%% operator for the quantum harmonic oscillator
i=sqrt(-1); Nmax=7; % 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
%%  measurement operators
varnu=pi/8; % dephasing per photon
ph0=pi/10; % phase of the second Ramsey pulse
th0=2*pi/sqrt(Nmax+1); % Rabi angle (trapping state condition)
nbar= 3; % goal photon number state
Mgm=diag(cos(th0*sqrt((0:1:Nmax))/2)); % u=-1
Mem=-diag(sin(th0*sqrt((1:1:Nmax))/2),1); %u=-1
Mg0=diag(cos((varnu*(0:1:Nmax)+ph0)/2)); % u=0
Me0=diag(sin((varnu*(0:1:Nmax)+ph0)/2)); % u= 0
Mgp=Mem';  % u=1
Mep=diag(cos(th0*sqrt((1:1:Nmax+1))/2)); % u=1
Lyap=diag(abs((0:1:Nmax)-nbar).^2);
%% initial state
Niter=300;
%rho=0*Id;rho(1,1)=1; ;% 0 photon; 
Vaccum=zeros(Nmax+1,1);Vaccum(1)=1; 
Psi= expm(sqrt(3)*(A-A'))*Vaccum; % initial state coherent state with 3 photons 
rho=Psi*(Psi');
rho=Id;rho=rho/trace(rho);
%% 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);
    % Lyapunov feedback law 
    Vp=trace(Lyap*(Mgp*rho*Mgp'+Mep*rho*Mep'));
    Vm=trace(Lyap*(Mgm*rho*Mgm'+Mem*rho*Mem'));
    Vmin=trace(Lyap*rho);Mg=Mg0;Me=Me0;u(ii)=0;
    if (Vp < Vmin);Vmin=Vp;Mg=Mgp;Me=Mep;u(ii)=1;end;
    if (Vm < Vmin);Vmin=Vm;Mg=Mgm;Me=Mem;u(ii)=-1;end;
    % 
    rhog=Mg*rho*Mg';rhoe=Me*rho*Me';
    if (rand < trace(rhog))
        y(ii)=-1;
        rho=rhog/trace(rhog);
    else
        y(ii)=1;
        rho=rhoe/trace(rhoe);
    end
end
%% graphics
figure(1),clf,
subplot(211);
bar([1:Niter],u);
ylabel('Control  u ');
axis([0 Niter+1 -1.02 1.02 ]);
ax = gca; ax.FontSize = 14;
subplot(212);
plot([1:Niter],Pop,'linewidth',4);
legend('0','1','2','3','4','5','6','7');
axis([1 Niter -0.02 1.02 ]);
ax = gca; ax.FontSize = 14;
ylabel('populations');
xlabel('step number');
title(['Target of ',num2str(nbar),' photon(s)']); 



