% dispersive measurement of a qubit
% pierre.rouchon@mines-paristech.fr, October 2018
clear all;
%% operators and parameters
i=sqrt(-1);
Sm=[0 0; 1 0];Sp=Sm';
Sx=Sm+Sp;Sz=Sp*Sm-Sm*Sp;Sy=-i*Sz*Sx;Id=Sx*Sx;
Tmes=1; g=1/Tmes; % measurement strength
L=sqrt(g)*Sz; % Lindblad operator for direct quadature  measurement
eta=0.7;% measurement efficiency
%% initialization
Tsim=3*Tmes; dt=Tmes/200; t=[0:dt:Tsim];
rho=Id/2 ; % initial state fully depolarized (no information)
z=0*t; Y=0*t;   % buffer for graphics
for k=1:1:length(t)
    % measurement outcome: 
    dW=randn*sqrt(dt); dY=sqrt(eta)*trace((L+L')*rho)*dt + dW;
    % measurement backaction
    M=Id-(L'*L/2)*dt+dY*L;
    rho=M*rho*(M') + (1-eta)*L*rho*(L')*dt;   rho=rho/trace(rho);
    z(k)=trace(rho*Sz); if (k>1);Y(k)=Y(k-1)+dY;end; % graphics
end
z=real(z);Y=real(Y); % avoid infinitesimal  imaginary  part
%% graphics
figure(1);clf;
subplot(211)
plot(t,Y,'linewidth',3)
title('Measurement signal y(t)','fontsize',16)
subplot(212)
plot(t,z,'-k','linewidth',3)
axis([min(t),max(t),-1.05,+1.05]);
xlabel('time'); title('tr(\rho \sigma_z)','fontsize',16)




