% Hamiltonian for the electromagnetic field (Sec. III.B)
%   H = hbar*omega*(a'a + 1/2 I)
% We set hbar = omega = 1 to express energies in units of hbar*omega.

clear;          % Clear memory
clc;            % Clear the command window/screen
d = 6;          % dimension of the operator
I = eye(d);     % identity matrix
hbar = 1;
omega = 1;      % frequency of the cavity field
A = diag(sqrt(1:d-1), 1);   % Annihilation operator
Ad = A';                    % Creation operator
H = hbar*omega*(Ad*A + (1/2)*I)   % Hamiltonian (diagonal: 0.5, 1.5, 2.5, ...)