from math import pi,exp

kB=1.38e-23
h=6.626e-34
uam=1.66e-27

m=28*1e-3*uam
T=300
c=2.998e+8

d=1.13e-10
mC=12*uam
mO=16*uam
m=mC*mO/(mC+mO)
I=m*d**2
hbar=h/(2*pi)

for J in range(0,7):
    a=(2*J+1)*hbar**2
    b=2*I*kB*T
    c=exp(-((J*(J+1)*hbar**2)/b))
    f=(a/b)*c
    print("Nj/N pt J=", J, f)