% 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
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
ph0=pi/7; % dephasing per photon
phR=0; % 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((ph0*(0:1:Nmax)+phR)/2)); % u=0
Me0=diag(sin((ph0*(0:1:Nmax)+phR)/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);
%% measurement imperfections
Eff=4/10; % efficiency
Err=1/10; % 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
%% initial state
rho=0*Id;rho(1,1)=1; Niter=800;% 0 photon; 
%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);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;
    % 
    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],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 ]); grid on;
ylabel('populations');
xlabel('step number');
title(['Target of ',num2str(nbar),' photon(s)']); 
%%
figure(2),clf;
subplot(211);
bar([1:Niter],u);
ylabel('Control  u ');
axis([0 Niter+1 -1.02 1.02 ]);
subplot(212);
moy = sum(N*Pop,1); 
dev = sqrt( sum(N^2*Pop)-moy.^2);
plot([1:Niter],moy,'k-',...
     [1:Niter],moy+dev,'r:',...
     [1:Niter],moy-dev,'r:','linewidth',2);
axis([0 Niter+1 -0.1 7.1 ]); grid on;
xlabel('step number');
title(['average number of photon(s) (solid line) with +/- standard deviation (dot lines) (target ',num2str(nbar),' photon(s))'])


