% Photon number distribution (Sec. IV.A)
% P_n = |<n|psi>|^2. Poissonian for a coherent state, monotonic for a
% thermal state. Two bar plots are produced.

clear;          % Clear memory
clc;            % Clear the command window/screen
d = 15;         % dimension of the field
I = eye(d);

%%%%%%%% Coherent state
alpha = 2;      % amplitude of the coherent state
Coh = 0;        % initialization
for x = 0:d-1
  Coh = Coh + exp(-norm(alpha)^2/2)*alpha^x/sqrt(prod(1:x))*I(:,x+1);
end
for n = 0:d-1
  PnC(n+1) = norm(I(:,n+1)'*Coh)^2;   % Probability of |n> in coherent state
end

%%%%%%%% Thermal state
nth = 0.5;      % average number of photons in the thermal state
RhoTh = 0;
for x = 0:d-1
  RhoTh = RhoTh + nth^(x)/(1+nth)^(x+1)*I(:,x+1)*I(:,x+1)';
end
for n = 0:d-1
  PnTh(n+1) = I(:,n+1)'*RhoTh*I(:,n+1);   % Probability of |n> in thermal state
end

%%%%%%%%%% Bar plots
n = 0:d-1;
figure(1)
bar(n, PnC)
title('Coherent state \alpha = 2');  xlabel('n'); ylabel('P_n')

figure(2)
bar(n, PnTh)
title('Thermal state n_{th} = 0.5'); xlabel('n'); ylabel('P_n')