% March,  2018
% pierre.rouchon@mines-paristech.fr
% dispersive measurement of a qubit
% stabilization of |e>
% realitistic model with decoherence and  detection efficiency
%
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; 
%% Qubit decoherence
Te=20*Tmes; % life time of |e> (z=+1)
Le=sqrt(1/Te)*Sm; % 
%% initialization
Tsim=200*Tmes; dt=Tmes/20; t=[0:dt:Tsim];
%rho=Id/trace(Id); % initial state
rho=(Id+(Sx-Sy)/sqrt(2))/2;
%rho=(Id)/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;
    % Feedback stabilization of z=1
    alpha=2/Tmes;beta=0.9*sqrt(8*eta*alpha);
    u(k)=alpha*2*y(k)+beta*(1-z(k));
    v(k)=0;
    %
    H=(u(k)*Sx+v(k)*Sy+Weg*Sz)/2;
    M=Id-(i*H+Le'*Le/2+L'*L/2)*dt+sqrt(eta)*dY*L;
    rho=M*rho*(M')+(1-eta)*L*rho*(L')*dt+Le*rho*(Le')*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(311)
plot(t,Y,'linewidth',2)
legend('Y');
subplot(312)
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]);
subplot(313)
plot(t,u,'linewidth',2)
legend('u');
xlabel('time')


