#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Somme d'état (ensemble canonique) de rotation (rotateur biatomique)

Les impressions sont à récrire avec l'instruction format() de python 3
"""

from math import exp    # on a juste besoin de l'exponentielle
import matplotlib.pyplot as plt  # directive d'importation standard de Matplotlib

T = 100. # (température réduite = T / Theta)
Zrot = 0.  # somme d'état
Jmax = 30  # valeur maximale de J
print('---------------------------------------------------')
print('|   J | dégénerescence | J*(J+1) |   terme de Z   |') 
print('---------------------------------------------------')
Js=[]
terms=[]
for J in range(Jmax+1):
    g = 2*J+1  # dégénerescence
    term= (2*J+1)*exp(-J*(J+1)/T)
    print('| %3d |      %4d      |  %4d   | %.8e |' % (J, g, J*(J+1), term))
    Zrot += term
    Js.append(J)
    terms.append(term)
print('---------------------------------------------------')

print( "Somme d'état à T=%3d jusqu'au niveau %2d : %.8e " % (T,Jmax,Zrot))

fig = plt.figure()
plt.title(u"Rotateur biatomique et somme d'état à T = "+str(T))
plt.plot(Js, terms, 'bo', label=r'$(2J+1)\exp{(-J(J+1)\theta_{rot}/T)}$')
plt.legend(loc='upper right')
plt.vlines(Js, [0], terms)
plt.xlabel(u"Niveaux J")  
plt.ylabel(u"Termes de la somme d'état")
eq = r'$Z_I^{Rot} = \sum_{J} (2J+1) \exp{(-J(J+1)\theta_{rot}/T)} = '+str(round(Zrot,3))+"$"
size = 14
x,y = 12,7
plt.text(x, y, eq, fontsize=size, clip_on=True)

plt.show()  # vue interactive de la figure