% Thermal state (Sec. II.D)
%   rho_th = 1/(1+nth) * sum_n (nth/(1+nth))^n |n><n|
% A mixed state: a square density matrix, diagonal in the number basis.

clear;          % Clear memory
clc;            % Clear the command window/screen
d = 10;
I = eye(d);
nth = 0.5;      % average number of photons in the thermal state
RhoTh = 0;      % initialization
for x = 0:d-1
  RhoTh = RhoTh + nth^(x)/(1+nth)^(x+1)*I(:,x+1)*I(:,x+1)';
end
RhoTh               % display the thermal state
N_th = trace(RhoTh) % check trace (should be 1)