% March,  2018
% 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
L=sqrt(1/Tmes)/2*Sz; % Lindblad operator for Y
eta=0.6; % efficiency 
Weg=2*pi/Tmes;
%% initialization
Tsim=10*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; Y=x; u=x;v=x; % buffer for graphics
%%
for k=1:1:length(t)
    %
    dW=randn*sqrt(dt);
    dY=sqrt(eta)*trace((L+L')*rho)*dt + dW;
    %
    x(k)=trace(rho*Sx);y(k)=trace(rho*Sy);z(k)=trace(rho*Sz);
    if (k>1);Y(k)=Y(k-1)+dY;end;
    %
    H=Weg*Sz/2;
    M=Id-(i*H+L'*L/2)*dt+sqrt(eta)*dY*L;
    rho=M*rho*(M')+(1-eta)*L*rho*(L')*dt;
    rho=rho/trace(rho);
end
%%
x=real(x);y=real(y);z=real(z);Y=real(Y);u=real(u);v=real(v);
%
figure(1);clf;
subplot(211)
plot(t,Y,'linewidth',2)
legend('Y');
subplot(212)
plot(t,x,'-r',t,y,'-b',t,z,'-k','linewidth',2)
legend('x','y','z');
axis([min(t),max(t),-1.05,+1.05]);
xlabel('time')


