% Dispersive measurement of the energy of a qubit
% stabilization of the exicted state via quantum-state feedback 
% based on Lyapunov techniques 
% model with finite life time of the exicted state
% perfect measurement
% March 2016
% pierre.rouchon@mines-paristech.fr
% alain.sarlette@inria.fr
%
clear all;
%% Pauli operators
i=sqrt(-1);
Sm=[0 0; 1 0];Sp=Sm'; % sigma_minus and sigma_plus
Sx=Sm+Sp; Sz=Sp*Sm-Sm*Sp; Sy=-i*Sz*Sx; Id=Sx*Sx;
%% quadrature measurements I and Q
g=1; Tg=1/g; % gamma
LI=sqrt(g/2)*Sz; % Lindblad operator for measurement signal I
LQ=i*sqrt(g/2)*Sz;  % Lindblad operator for  measurement signal 
T1=20*Tg;Le=sqrt(1/T1)*Sm;
%% initialization
Tsim=100*Tg; dt=Tg/20; t=(0:dt:Tsim)'; Niter=length(t);
rho=(Id+Sx)/2; %  x=1 corresponds to the pure state (|g>+|e>)/sqrt(2)
x=zeros(Niter,1);y=zeros(Niter,1);z=zeros(Niter,1); % buffer for graphics
I=zeros(Niter,1); Q=zeros(Niter,1);  % buffer for  measured outputs
u=zeros(Niter,1); v=zeros(Niter,1);  % buffer for  control inputs
for ii=1:Niter
    dWI=randn*sqrt(dt);dWQ=randn*sqrt(dt); % Wiener increment
    % measured outputs
    dI=trace((LI+LI')*rho)*dt + dWI; %  increment of measured signal I
    dQ=trace((LQ+LQ')*rho)*dt + dWQ; %  increment of measured signal Q
    % control inputs 
    uii=2*g*(1-trace(rho*Sz)); if (trace(rho*Sy)<0); uii=-uii; end;
    vii=-2*g*(1-trace(rho*Sz)); if (trace(rho*Sx)<0); vii=-vii; end;
    % buffer for graphics
    x(ii)=trace(rho*Sx);y(ii)=trace(rho*Sy);z(ii)=trace(rho*Sz);
    u(ii)=uii;v(ii)=vii;
    if (ii>1);I(ii)=I(ii-1)+dI;Q(ii)=Q(ii-1)+dQ;end;
    % positivity and trace  perserving numerical integration scheme
    H=(uii*Sx+vii*Sy)/2;
    M_dy = Id -(i*H+LI'*LI/2+LQ'*LQ/2+Le'*Le/2)*dt + dI*LI + dQ*LQ;
    rho=M_dy*rho*(M_dy')+ dt*Le*rho*(Le');
    rho=rho/trace(rho);
end
x=real(x);y=real(y);z=real(z);I=real(I);Q=real(Q);u=real(u);v=real(v);
%%
figure(1);clf;
subplot(311)
plot(g*t,u/g,g*t,v/g,'linewidth',2)
legend('u/\gamma','v/\gamma');
subplot(312)
plot(g*t,x,'-r',g*t,y,'-b',g*t,z,'-k','linewidth',2)
legend('x','y','z');
axis([min(g*t),max(g*t),-1.05,+1.05]);
subplot(313)
plot(g*t,I,g*t,Q,'linewidth',2)
legend('I','Q');
xlabel('\gamma t')



