# QuTiP lecture: Single-Atom-Lasing

I have tried to run one of the examples in QuTip lectures: single atom lasing. Although I have used the prepared code represented in reference, I have received error and I have not been able to achieve the correct result. I have focused on the "Steady state average photon occupation in cavity as a function of pump rate". The code that I have used is as follows

from qutip import *
import matplotlib.pyplot as plt
import numpy as np
from numpy import sqrt
w0 = 1.0  * 2 * pi  # cavity frequency
wa = 1.0  * 2 * pi  # atom frequency
g  = 0.05 * 2 * pi  # coupling strength
kappa = 0.04        # cavity dissipation rate
gamma = 0.00        # atom dissipation rate
Gamma = 0.35        # atom pump rate
N = 50              # number of cavity fock states
n_th_a = 0.0        # avg number of thermal bath excitation

tlist = np.linspace(0, 150, 101)

# intial state

# operators
a  = tensor(destroy(N), qeye(2))
sm = tensor(qeye(N), destroy(2))
sx = tensor(qeye(N), sigmax())

# Hamiltonian
H = w0 * a.dag() * a + wa * sm.dag() * sm + g * (a.dag() + a) * sx

def calulcate_avg_photons(N, Gamma):

# collapse operators
c_ops = []

rate = kappa * (1 + n_th_a)
if rate > 0.0:
c_ops.append(sqrt(rate) * a)

rate = kappa * n_th_a
if rate > 0.0:
c_ops.append(sqrt(rate) * a.dag())

rate = gamma
if rate > 0.0:
c_ops.append(sqrt(rate) * sm)

rate = Gamma
if rate > 0.0:
c_ops.append(sqrt(rate) * sm.dag())

# Ground state and steady state for the Hamiltonian: H = H0 + g * H1

# cavity photon number
n_cavity = expect(a.dag() * a, rho_ss)

# cavity second order coherence function
g2_cavity = expect(a.dag() * a.dag() * a * a, rho_ss) / (n_cavity ** 2)

return n_cavity, g2_cavity

Gamma_max = 2 * (4*g**2) / kappa
Gamma_vec = np.linspace(0.1, Gamma_max, 50)

n_avg_vec = []
g2_vec = []

for Gamma in Gamma_vec:
n_avg, g2 = calulcate_avg_photons(N, Gamma)
n_avg_vec.append(n_avg)
g2_vec.append(g2)

fig, axes = plt.subplots(1, 1, figsize=(12,6))

axes.plot(Gamma_vec * kappa / (4*g**2), n_avg_vec, color="blue", alpha=0.6, label="numerical")

axes.set_xlabel(r'$$\Gamma\kappa/(4g^2)$$', fontsize=18)
axes.set_ylabel(r'Occupation probability $$\langle n \rangle$$', fontsize=18)
axes.set_xlim(0, 2);

fig, axes = plt.subplots(1, 1, figsize=(12,6))

axes.plot(Gamma_vec * kappa / (4*g**2), g2_vec, color="blue", alpha=0.6, label="numerical")

axes.set_xlabel(r'$$\Gamma\kappa/(4g^2)$$', fontsize=18)
axes.set_ylabel(r'$$g^{(2)}(0)$$', fontsize=18)
axes.set_xlim(0, 2)
axes.text(0.1, 1.1, "Lasing regime", fontsize=16)
axes.text(1.5, 1.8, "Thermal regime", fontsize=16);



The error that I received is

TypeError
Traceback (most recent call last)
<ipython-input-3-a1345bdca7b1> in <module>()
50 g2_vec = []
51 for Gamma in Gamma_vec: ---> 52
n_avg, g2 = calulcate_avg_photons(N, Gamma) 53
n_avg_vec.append(n_avg) 54
g2_vec.append(g2)
TypeError: 'NoneType' object is not iterable


Can you please guide me how to fix the problem and achieve the correct result?