% July ,  2019
% pierre.rouchon@mines-paristech.fr
% dispersive measurement of a qubit
%
clear all;
%% Pauli operators
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;
%% quadrature measurements
Tmes=1; % measurement time constant
LI=sqrt(1/(Tmes))*Sz; % Lindblad operator for I
eta=0.4; % efficiency 
%% Qubit decoherence
Te=10*Tmes; % life time of |e> (z=+1)
Le=sqrt(1/Te)*Sm; % 
%% initialization
Tsim=30*Tmes; dt=Tmes/20; t=[0:dt:Tsim];
%rho=Id/trace(Id); % initial state
rho=(Id+(Sx+Sy)/sqrt(2))/2;
x=0*t;y=x;z=x; I=x; % buffer for graphics
for k=1:1:length(t)
    %
    dWI=randn*sqrt(dt);dWQ=randn*sqrt(dt);
    dI=sqrt(eta)*trace((LI+LI')*rho)*dt + dWI;
     %
    x(k)=trace(rho*Sx);y(k)=trace(rho*Sy);z(k)=trace(rho*Sz);
    if (k>1);I(k)=I(k-1)+dI;end;
   
    M=Id-(Le'*Le/2+LI'*LI/2)*dt+sqrt(eta)*dI*LI;
    rho=M*rho*(M')+(1-eta)*LI*rho*(LI')*dt+Le*rho*(Le')*dt;
    rho=rho/trace(rho);
end
x=real(x);y=real(y);z=real(z);I=real(I);
%%
figure(1);clf;
subplot(211)
plot(t,I, 'linewidth',2)
legend('measurement I');
subplot(212)
plot(t,x,'-r',t,y,'-b',t,z,'-k','linewidth',4)
legend('x','y','z');
axis([min(t),max(t),-1.05,+1.05]);
xlabel('time')


