% LKB photon box: QND measurement of photon
% November  2023
% pierre.rouchon@minesparis.psl.eu
% density operator formulation, ideal case
%
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; % tuned by the phase of the second pi/2 Ramsey pulse
Mg=diag(cos((varnu*(0:1:Nmax)+ph0)/2));
Me=diag(sin((varnu*(0:1:Nmax)+ph0)/2));
%% simulation parameter
Niter=300; 
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/(Nmax+1);
%% simulation loop
Pop=zeros(Nmax+1,Niter);y=zeros(1,Niter); % buffer for graphics
for ii=1:Niter
    Pop(:,ii)=diag(rho);
    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],y);
ylabel('Measurement outcome');
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');



