from matplotlib import pyplot as plt
import numpy as np
k=0.695 #In wavenumbers/K. We can change it to other units if we want to use $\varepsilon$ in other units such as 1.38x10^-23 J/K
Na=6.022045*10**23
h=6.262e-34
c=2.9979e10
def q(eps,T):
sum=0
for i in range(0,len(eps)):
sum+=np.exp(-eps[i]/k/T)
return sum
T=np.arange(1,1500,1)
#@title Default title text
upperlevelenergy = 500#@param {type:"integer"}
eps1=[0,upperlevelenergy] #energy levels at 0 and 50 cm^-1
plt.plot(T,q(eps1,T))
plt.ylim(1,2)