% November,  2021
% pierre.rouchon@mines-paristech.fr
% homodyne dispersive measurement of a qubit
%
clear all;
%% Pauli operators
Sm=[0 0; 1 0];Sp=Sm';
Sx=Sm+Sp;Sz=Sp*Sm-Sm*Sp;Sy=-1i*Sz*Sx;Id=Sx*Sx;
%% quadrature measurements
Tmes=1; g=1/Tmes;  % measurement strength
eta=0.5;  % measurement efficiency 
L=sqrt(g/4)*Sz; % Lindblad operator for homodyne measurement 
weg=3/Tmes;
Tb=10*Tmes; nth=0.1; 
Lp=sqrt(nth/Tb)*Sp; % gain of a thermal photon from ground state 
Lm=sqrt((1+nth)/Tb)*Sm; % lost of a photon from excited state  
%% initialization
Tsim=100*Tmes; dt=Tmes/20; t=[0:dt:Tsim];
U=expm(-1i*(dt/2)*weg*Sz);
%rho=Id/trace(Id); % initial state
 rho=(Id+(Sx-Sy)/sqrt(2))/2;
x=0*t;y=x;z=x; Y=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
    %
    M=Id-(L'*L+Lm'*Lm+Lp'*Lp)*dt/2+dY*L; 
    rho= U*rho*U';
    rho=M*rho*(M')+ dt*((1-eta)*L*rho*L'+Lm*rho*Lm'+Lp*rho*Lp');
    rho= U*rho*U';
    rho=rho/trace(rho);
end
x=real(x);y=real(y);z=real(z);Y=real(Y);
%%
figure(1);clf;
subplot(211)
plot(t,Y,'linewidth',2)
legend('Y'); title('Homodyne measurement signal Y_t')
subplot(212)
plot(t,x,'-r',t,y,'-b',t,z,'-k','linewidth',4)
legend('x','y','z');
axis([min(t),max(t),-1.05,+1.05]);
xlabel('time')
title('Bloch-sphere coordinates')



