% Homodyne measurement  of Sz for  a  qubit
% wave function  formulation in the ideal case 
% October 2018 pierre.rouchon@mines-paristech.fr
clear all;
%% qubit  operators
i=sqrt(-1); ket_g=[1 0]'; ket_e=[0 1]';
Sz=ket_e*(ket_e')-ket_g*(ket_g'); Id=Sz*Sz;
Tmes=1; g=1/Tmes; % measurement strength
L=sqrt(g)*Sz; % Lindblad operator
%% simulation loop
Tsim=3*Tmes; dt=Tmes/100; t=[0:dt:Tsim];
psi=(ket_g+ket_e)/sqrt(2) ; % initial state 
z=0*t; Y=0*t;   % buffer for graphics
for ii=1:1:length(t)
    dW=sqrt(dt)*randn; dY= dt*(psi'*(L+L')*psi)+dW;  %measurement outcome
    M_dY=Id-L'*L*dt/2 + dY*L; 
    psi=M_dY*psi; psi=psi/norm(psi); 
    if (ii>1);Y(ii)=Y(ii-1)+dY;end % integration of the measurement outcome
    z(ii)= psi'*Sz*psi;
end
%% graphics
Y=real(Y);z=real(z);
figure(1);clf;
subplot(211)
plot(t,Y,'linewidth',4)
title('Homodyne measurement signal y(t)','fontsize',16)
subplot(212)
plot(t,z,'-k','linewidth',4)
axis([min(t),max(t),-1.05,+1.05]);
xlabel('time');
title('<\psi|\sigma_z|\psi> ','fontsize',18)




