% LKB photon box: 
% quantum trajectories associated to 
%     QND measurement of photon
% ideal  model  without cavity decoherence
% initial state: Id
%
% March 2015
% 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
Mg=diag(cos((ph0*(0:1:Nmax)+phR)/2));
Me=diag(sin((ph0*(0:1:Nmax)+phR)/2));
%% Lyapunov feedback  
nbar=3; % goal state
F0=diag([0.6  0.4 0.2 0 0.25 0.47 0.65 0.8]); 
F1=(A'-A)*F0-F0*(A'-A);
F2=(A'-A)*F1-F1*(A'-A);
umax=0.05;
%% simulation parameter
Niter=300; 
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);
    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
    % feedback
    a1=trace(rho*F1);a2=trace(rho*F2);
    if (a2>0)
        u(ii)=min(umax,max(-umax,a1/a2));
    else
        if (a1 < 0);
            u(ii)=-umax;
        else
            u(ii)=umax;
        end
    end
    % displacement 
    Du=expm(u(ii)*(A'-A));
    rho=Du*rho*(Du');
    rho=rho/trace(rho);
end
%% graphics
figure(1),clf,
subplot(311);
bar([1:Niter],y);
ylabel('measurement outcome');
axis([0 Niter+1 -1.02 1.02 ]);
subplot(312);
plot([1:Niter],Pop,'linewidth',2);
legend('0','1','2','3','4','5','6','7');
axis([0 Niter+1 -0.02 1.02 ]);
ylabel('populations');
subplot(313)
plot([1:Niter],u,[1 Niter],-[umax umax],'k--',[1 Niter],[umax umax],'k--','linewidth',2);
axis([0 Niter+1 -umax*1.1 umax*1.1 ])
ylabel ('control u')
xlabel('step number');



