% November  2023
% pierre.rouchon@minesparis.psl.eu
% dispersive measurement of a qubit
% homodyne  detection
% ideal model
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
LI=sqrt(g)*Sz; % Lindblad operator for  detection quadrature I
%% 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;dI=x; I=x;   % buffer for graphics
for k=1:1:length(t)
    %
    dWI=randn*sqrt(dt);
    dI(k)=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 % 
    %
    M=Id-(LI'*LI/2)*dt+dI(k)*LI;
    rho=M*rho*(M');
    rho=rho/trace(rho);
end
x=real(x);y=real(y);z=real(z);I=real(I);dI=real(dI);
%%
figure(1);clf;
subplot(311)
plot(t,dI,'linewidth',2)
legend('dI','dQ');
ax = gca; ax.FontSize = 14;
subplot(312)
plot(t,I,'linewidth',2)
legend('I','Q');
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;





