% November  2023
% pierre.rouchon@minesparis.psl.eu
% dispersive measurement of a qubit
% homodyne detection
% feedback stabilization of |e>  based on 
% Koji Tsumura. Global stabilization of n-dimensional quantum spin systems via continuous feedback. 
% In American Control Conference, ACC 2007, New York, NY, USA, 9-13 July, 2007, pages 2129–2134. IEEE, 2007.
%
% realistic  model
%
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
T1=10*Tmes; % life time of |e> (z=+1)
L1=sqrt(1/T1)*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;  dI=x; u=x; % buffer for graphics
for k=1:1:length(t)
    %
    dWI=randn*sqrt(dt);
    dI(k)=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(k);end
    % Feedback stabilization of z=1
    alpha=1; beta=2*sqrt(alpha); u(k)=alpha*y(k)+beta*(1-z(k)); 
    %
    H=(u(k)*Sx)/2;
    M=Id-(i*H+L1'*L1/2+LI'*LI/2)*dt+sqrt(eta)*dI(k)*LI;
    rho=M*rho*(M')+(1-eta)*LI*rho*(LI')*dt+L1*rho*(L1')*dt;
    rho=rho/trace(rho);
end
x=real(x);y=real(y);z=real(z);I=real(I);dI=real(I);u=real(u);
%%
figure(1);clf;
subplot(311)
plot(t,I,'linewidth',2)
legend('I');
ax = gca; ax.FontSize = 14;
subplot(312)
plot(t,u,'linewidth',2)
legend('u');
ax = gca; ax.FontSize = 14;
subplot(313)
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')
ax = gca; ax.FontSize = 14;



