% 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; % measurement time constant
LI=sqrt(1/(2*Tmes))*Sz; % Lindblad operator for I
LQ=i*sqrt(1/(2*Tmes))*Sz;  % Lindblad operator for Q
eta=0.4; % efficiency 
%% Qubit decoherence
Te=10*Tmes; % life time of |e> (z=+1)
Le=sqrt(1/Te)*Sm; % 
%% 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; I=x; Q=x; u=x;v=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;
    dQ=sqrt(eta)*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;
    % Feedback stabilization of z=1
    gain=(1-z(k))/Tmes;
    u(k)=gain; if(y(k)<0);u(k)=-gain;end;
    v(k)=-gain; if(x(k)<0);v(k)=gain;end;
    %
    H=(u(k)*Sx+v(k)*Sy)/2;
    M=Id-(i*H+Le'*Le/2+LI'*LI/2+LQ'*LQ/2)*dt+sqrt(eta)*dI*LI+sqrt(eta)*dQ*LQ;
    rho=M*rho*(M')+(1-eta)*LI*rho*(LI')*dt+(1-eta)*LQ*rho*(LQ')*dt+Le*rho*(Le')*dt;
    rho=rho/trace(rho);
end
x=real(x);y=real(y);z=real(z);I=real(I);Q=real(Q);u=real(u);v=real(v);
%%
figure(1);clf;
subplot(211)
plot(t,u,t,v,'linewidth',2)
legend('u','v');
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')


