% LKB photon box: QND measurement of photon
% October 2014
% pierre.rouchon@mines-paristech.fr
%
clear all;
%% operator for the quantum harmonic oscillator
i=sqrt(-1); Nmax=7; % truncation to Nmax
Id=eye(Nmax+1);
N=diag((0:1:Nmax)); % photon  number operator
A=diag( sqrt((1:1:Nmax)),1); %  operator
%%  measurement operators
ph0=pi/8; % dephasing per photon
phR=pi/10; % phase of the second Ramsey pulse
Mg=diag(cos((ph0*(0:1:Nmax)+phR)/2));
Me=diag(sin((ph0*(0:1:Nmax)+phR)/2));
%% measurement imperfections
Eff=0.4; % efficiency
Err=0.1; % error rate
%% cavity decoherence  
nth=0.05; % thermal photon
Tcav=140e-3; % photon life time (s)
Tsampling=80e-6; % sampling period (s)
Lm=sqrt((1+nth)/Tcav)*A; % Lindblad operator associated to photon lost 
Lp=sqrt(nth/Tcav)*(A'); % Lindblad operator associated to photon gain
M0=Id-Tsampling*(Lm'*Lm+Lp'*Lp)/2; % Kraus operator for  zero photon annihilation
Mm=sqrt(Tsampling)*Lm; % Kraus operator to one photon annihilation
Mp=sqrt(Tsampling)*Lp; % Kraus operator to one photon creation
%% simulation parameter
Niter=5000;
rho=0*Id;rho(Nmax,Nmax)=1;    % 6 photons
%buf=expm(sqrt(3)*(A-A'))*Vacuum; rho=buf*buf';
%rho=Id/trace(Id); % initial state (no information) 

%% simulation loop
Pop=zeros(Nmax+1,Niter);y=zeros(1,Niter); % buffer for graphics
for ii=1:Niter
    Pop(:,ii)=diag(rho);
    rhog0=Mg*rho*(Mg');rhoe0=Me*rho*(Me');
    rhog=(1-Err)*rhog0+Err*rhoe0;
    rhoe=Err*rhog0+(1-Err)*rhoe0;
    if(rand < Eff)
        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=M0*rho*(M0')+Mm*rho*(Mm')+ Mp*rho*(Mp');
    rho=rho/trace(rho);
end
%% graphics
figure(1),clf,
subplot(211);
bar((1:Niter),y);
ylabel('Measurement outcome');
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');

figure(2),clf;
[buf,Nphoton]=max(Pop);Nphoton=Nphoton-1; 
 plot([1:Niter],Nphoton,'linewidth',2);
%plot([1:Niter],sum(N*Pop,1),'linewidth',2);
axis([0 Niter+1 -0.1 7.1 ]);
xlabel('step number');
title('Most likely photon-number')






