% Average number of photons (Sec. IV.B)
%   <a'a> = <psi|a'a|psi>  (pure states),   Tr(a'a rho)  (mixed states)
% Number state |4> -> 4, coherent (|alpha|^2 = 3) -> 3, thermal (nth) -> nth.

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

%%%% Number operator
A = diag(sqrt(1:d-1), 1);   % Annihilation operator
Ad = A';                    % Creation operator
AdA = Ad*A;                 % Number operator

%% Number state |4>
Ket4 = I(:,5);                  % four-photon state |4>
AdANumber = Ket4'*AdA*Ket4      % average photon number in the number state

%% Coherent state
alpha = sqrt(3);    % 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
AdACoherent = Coh'*AdA*Coh      % average photon number in coherent state (= |alpha|^2)

%%% Thermal state
nth = 0.85;     % assumed 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
AdAThermal = trace(AdA*RhoTh)   % average photon number in thermal state (= nth)