% Reservoir  stabilisation  of  a cat qubit 
% via the decoherence  dynamics 
%    drho/dt =  D_L0(rho) +   k1 D_L1(\rho) 
% with Lindblad operators 
% L0=a^2 - alpha^2,  alpha real large,  
% L1= a   with  k1 \ll \alpha^2  
% 
% December 2021
% pierre.rouchon@mines-paristech.fr
%
clear all; 
%% parameters
alpha=sqrt(12.5); 
k1= 1/10; % error rate for L1
Nmax=floor(alpha^2+15*alpha); % truncation to Nmax
%% operator and matrices
Id=eye(Nmax+1); Vacuum=zeros(Nmax+1,1);Vacuum(1)=1;
N=diag((0:1:Nmax)); % photon  number operator
A=diag( sqrt((1:1:Nmax)),1);% annihilation operator
Q=(A+A')/sqrt(2); P=(A-A')/sqrt(2)/1i; % [Q,P]= i 
%% 
L0= A*A- alpha*alpha*Id; %  two-photon pumping
L1=sqrt(k1)*A; %  error operator 
%% code space and associated operators
[V,T]=schur(P); D=diag(T); Dalpha= V*diag(exp(-1i*sqrt(2)*alpha*D))*V'; %shift of alpha along q
Ket0 = Dalpha*Vacuum; % |alpha>
Ket1 = Dalpha'*Vacuum; % |-alpha>
CatEven=Ket0+Ket1; CatEven=CatEven/norm(CatEven); % 
CatOdd=Ket0-Ket1; CatOdd=CatOdd/norm(CatOdd); % 
Pcode=CatEven*CatEven' + CatOdd*CatOdd';   % Projection on code space 
%% initial condition 
Psi=Ket0; % Vacuum CatEven  Ket0  initial pure state 
rho0=Psi*Psi'; rho0=rho0/trace(rho0); WignerPlot(rho0,5+alpha*sqrt(2),1); 
%%
figure(2);clf
points = (-5*alpha:0.01:5*alpha);
q_quad = q_quadrature(points,Psi); p_quad = p_quadrature(points,Psi);
theta=0.5*angle(-sum((imag(q_quad).^2-real(q_quad).^2)/2+1i*real(q_quad).*imag(q_quad)));
q_quad=exp(1i*theta)*q_quad; p_quad=exp(1i*theta)*p_quad; 
if(sum(real(q_quad))<0); q_quad=-q_quad; p_quad=-p_quad; end
plot(points,[real(q_quad),imag(q_quad)],'linewidth',2);
xlabel('$q$','Interpreter','latex');
legend(["$Re(\psi(q))$","$Im(\psi(q))$"], 'Interpreter','latex')
title('Initial pure state')
grid on
drawnow;  
%% Simulation from  initial state  rho0
Tsim= 10/(alpha^2); % convergence time to code space around 1/alpha^2
dt=0.001/(alpha^2); nt=floor(Tsim/dt);
M0=Id-(L0'*L0+L1'*L1)*dt/2;
fid=zeros(nt,1);pcode=zeros(nt,1);
rho=rho0; 
disp('Simulation  from 0 to Tsim'); Int=floor(1.01*nt/100);
for it=1:nt
     fid(it)=trace(rho*rho0); pcode(it)=trace(Pcode*rho); 
     rho=M0*rho*M0'+dt*L0*rho*L0'+dt*L1*rho*L1'; % Kraus map numerical scheme 
     rho=rho/trace(rho);
     if (it/Int==floor(it/Int));disp([' progress ',num2str(it/Int),'%']);end    
end
pcode=real(pcode);fid=real(fid);
%%
[V,T]=schur(rho); D=abs(diag(T));[pop,I]=sort(D,'descend'); Vrho=V(:,I); 
% save('rho_catqubit.mat', 'pop', 'Vrho');  
WignerPlot(rho,5+alpha*sqrt(2),10);
%%
figure(11),clf; 
plot(dt*(0:nt-1)*alpha^2,fid,'linewidth',2);
xlabel('$\alpha^2 t$','Interpreter','latex');
title('Fidelity between $\rho_0$ and $\rho_t$','Interpreter','latex');
%%
figure(12);clf
subplot(411)
Psi=Vrho(:,1);
points = (-5*alpha:0.01:5*alpha);
q_quad = q_quadrature(points,Psi); p_quad = p_quadrature(points,Psi);
theta=0.5*angle(-sum((imag(q_quad).^2-real(q_quad).^2)/2+1i*real(q_quad).*imag(q_quad)));
q_quad=exp(1i*theta)*q_quad; p_quad=exp(1i*theta)*p_quad; 
if(sum(real(q_quad))<0); q_quad=-q_quad; p_quad=-p_quad; end
plot(points,[real(q_quad),imag(q_quad)],'linewidth',2);
xlabel('$q$','Interpreter','latex');
legend(["$Re(\psi_1(q))$","$Im(\psi_1(q))$"], 'Interpreter','latex')
title(['First dominant pure state $p_1=$',num2str(pop(1))],'Interpreter','latex')
grid on
subplot(412)
Psi=Vrho(:,2);
q_quad = q_quadrature(points,Psi); p_quad = p_quadrature(points,Psi);
theta=0.5*angle(-sum((imag(q_quad).^2-real(q_quad).^2)/2+1i*real(q_quad).*imag(q_quad)));
q_quad=exp(1i*theta)*q_quad; p_quad=exp(1i*theta)*p_quad; 
if(sum(real(q_quad))<0); q_quad=-q_quad; p_quad=-p_quad; end
plot(points,[real(q_quad),imag(q_quad)],'linewidth',2);
xlabel('$q$','Interpreter','latex');
legend(["$Re(\psi_2(q))$","$Im(\psi_2(q))$"], 'Interpreter','latex')
title(['Second dominant pure state $p_2=$',num2str(pop(2))],'Interpreter','latex')
grid on
subplot(413)
Psi=Vrho(:,3);
q_quad = q_quadrature(points,Psi); p_quad = p_quadrature(points,Psi);
theta=0.5*angle(-sum((imag(q_quad).^2-real(q_quad).^2)/2+1i*real(q_quad).*imag(q_quad)));
q_quad=exp(1i*theta)*q_quad; p_quad=exp(1i*theta)*p_quad; 
if(sum(real(q_quad))<0); q_quad=-q_quad; p_quad=-p_quad; end
plot(points,[real(q_quad),imag(q_quad)],'linewidth',2);
xlabel('$q$','Interpreter','latex');
legend(["$Re(\psi_3(q))$","$Im(\psi_3(q))$"], 'Interpreter','latex')
title(['Third dominant pure state $p_3=$',num2str(pop(3))],'Interpreter','latex')
grid on
subplot(414)
Psi=Vrho(:,4);
q_quad = q_quadrature(points,Psi); p_quad = p_quadrature(points,Psi);
theta=0.5*angle(-sum((imag(q_quad).^2-real(q_quad).^2)/2+1i*real(q_quad).*imag(q_quad)));
q_quad=exp(1i*theta)*q_quad; p_quad=exp(1i*theta)*p_quad; 
if(sum(real(q_quad))<0); q_quad=-q_quad; p_quad=-p_quad; end
plot(points,[real(q_quad),imag(q_quad)],'linewidth',2);
xlabel('$q$','Interpreter','latex');
legend(["$Re(\psi_4(q))$","$Im(\psi_4(q))$"], 'Interpreter','latex')
title(['Fourth dominant pure state $p_4=$',num2str(pop(4))],'Interpreter','latex')
grid on
return

%% Computing q-quadrature given Fock states populations
% Lev-Arcady Sellem <lev-arcady.sellem@inria.fr>
%
function [q_quad] = q_quadrature(q_points, psi)
q_points = reshape(q_points,[numel(q_points),1]); % We need a column vector here

q_quad = zeros(numel(q_points),1);
F_0 = zeros(numel(q_points),1); % Pre-allocating
F_1 = zeros(numel(q_points),1);
%F_n = zeros(numel(x_points),1); % No point in pre-allocating here because
%the first computation of F_n will allocate the same space
for n = 0:numel(psi)-1
    if n == 0
        F_0 = pi^(-0.25) * exp(-q_points.^2/2) ;
        q_quad = q_quad + psi(n+1)*F_0;
    elseif n == 1
        F_1 = sqrt(2)*q_points.*F_0 ;
        q_quad = q_quad + psi(n+1)*F_1;
    else
        F_n = sqrt(2/n)*q_points.*F_1 - sqrt((n-1)/n)*F_0;
        q_quad = q_quad + psi(n+1)*F_n;
        F_0 = F_1; F_1 = F_n; % Updating so that F_0,F_1 retain the last two polynomials
    end
end
end
%%
function [p_quad] = p_quadrature(p_points, psi)
p_points = reshape(p_points,[numel(p_points),1]); % We need a column vector here

p_quad = zeros(numel(p_points),1);
F_0 = zeros(numel(p_points),1); % Pre-allocating
F_1 = zeros(numel(p_points),1);
%F_n = zeros(numel(x_points),1); % No point in pre-allocating here because
%the first computation of F_n will allocate the same space
for n = 0:numel(psi)-1
    if n == 0
        F_0 = pi^(-0.25) * exp(-p_points.^2/2);
        p_quad = p_quad + psi(n+1)*F_0;
    elseif n == 1
        F_1 = -1i*sqrt(2)*p_points.*F_0;
        p_quad = p_quad + psi(n+1)*F_1;
    else
        F_n = -1i*sqrt(2/n)*p_points.*F_1 + sqrt((n-1)/n)*F_0;
        p_quad = p_quad + psi(n+1)*F_n;
        F_0 = F_1; F_1 = F_n; % Updating so that F_0,F_1 retain the last two polynomials
    end
end
end
%%
function  WignerPlot(rho,qpfig,nofig)
% plot on figure(nofig)
%  the  wigner function of the 30  dominant eigen-states for   density operator rho
%  with Nmax >30;
% pierre.rouchon@mines-paristech.fr
% october 2021
%
epsilon=1.e-20; % for a bug in wvd when rho=vaccum exactly 
[V,T]=schur(rho); D=abs(diag(T));[pop,I]=sort(D,'descend'); Vrho=V(:,I);
Nmax=length(pop)-1;
Npop=20; if (Nmax < Npop); disp('  Error in WignerPlot: Nmax < 20');  return; end
%Q=diag( sqrt((1:1:Nmax)),1); Q=(Q+Q')/sqrt(2); Q2av=real(trace(Q*Q*rho));
nW=2^10-1; dqp=sqrt(pi/(nW)); qpmax=sqrt(pi*nW/4); q=(-qpmax:dqp:+qpmax)'; 
Psi=Vrho(:,1)+epsilon*(rand+1i*rand);
Psi_q = q_quadrature(q,Psi); 
[W1,f,t] = wvd(Psi_q); W1=pop(1)*W1; % positive value of p
[W2,f,t] = wvd(Psi_q'); W2=pop(1)*W2;  % negative value of p
for k=2:1:Npop
    Psi_q = q_quadrature(q,Vrho(:,k)+epsilon*(rand+1i*rand));
    [Wk,f,t] = wvd(Psi_q); W1=W1+pop(k)*Wk;
    [Wk,f,t] = wvd(Psi_q'); W2=W2+pop(k)*Wk;
end
[n,m]=size(W1);n2=floor(n/2);
W=[flipud(W1((2:n2),:));W2((1:n2),:)];
f=[flipud(f((2:n2),:));-f(1:n2)];
q=qpmax*(2*t/(t(end)-t(1))-1); q=q+(q(2)-q(1))*1.5; p=2*qpmax*f;
W=W/(sum(sum(W))*abs((p(2)-p(1))*(q(2)-q(1))));  
figure(nofig),clf; 
imagesc(q,p,W)
xlabel('$q$','Interpreter','latex');
ylabel('$p$','Interpreter','latex');
colorbar
grid on
axis([-qpfig qpfig -qpfig qpfig]);axis square
caxis([-0.301,0.301])
colormap('jet');
end