% November,  2015
% 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; g=1/Tmes; % measurement strength
LI=sqrt(g/2)*Sz; % Lindblad operator for I
LQ=i*sqrt(g/2)*Sz;  % Lindblad operator for Q
%% initialization
Tsim=10*Tmes; dt=Tmes/20; t=[0:dt:Tsim];
rho=Id/trace(Id); % initial state
%rho=(Id-Sx)/2;
x=0*t;y=x;z=x; I=x; Q=x;  % buffer for graphics
for k=1:1:length(t)
    %
    dWI=randn*sqrt(dt);dWQ=randn*sqrt(dt);
    dI=trace((LI+LI')*rho)*dt + dWI;
    dQ=trace((LQ+LQ')*rho)*dt + dWQ;
    %
    x(k)=trace(rho*Sx);y(k)=trace(rho*Sy);z(k)=trace(rho*Sz);
    if (k>1);I(k)=I(k-1)+dI;Q(k)=Q(k-1)+dQ;end;
    %
    M=Id-(LI'*LI/2+LQ'*LQ/2)*dt+dI*LI+dQ*LQ;
    rho=M*rho*(M');
    rho=rho/trace(rho);
end
x=real(x);y=real(y);z=real(z);I=real(I);Q=real(Q);
%%
figure(1);clf;
subplot(211)
plot(t,I,t,Q,'linewidth',2)
legend('I','Q');
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')



