Show the luminosity by Planck law, numerically integrate to compute the S-B constant for different temperatures

Planck_and_Stefan-Boltzmann_laws

#!/usr/bin/env python
#-*- coding: utf-8 -*-

## Import common moduli
import matplotlib, sys, os, time
import matplotlib.pyplot as plt
import numpy as np
from   scipy.constants import c, h, hbar, pi, k, e
import scipy.integrate

f = np.linspace(c/100000e-9, c/100e-9, 100000) ## define the spectral range for integration
for T in [300., 5780.]:         ## define the temperatures you wish to inspect
    #T = 9e2
    L=2*h*f**3/c**2 / ((np.exp(h*f/k/T))-1)
    totalL = scipy.integrate.trapz(L, f) * pi    # integrate by frequency
    print "Radiosity at T=%g K (integrated from f=%g to %g): %g W/m^2" % (T, np.min(f), np.max(f), totalL)
    print "\t\tCompared to the radiosity of black body over entire spectrum ratio (Stefan-Boltzmann law): %.4g" % (totalL / 5.670367e-8 /T**4)
    plt.plot(f,L, label="T = %g K" % T)

plt.xlabel(u"Radiation frequency [Hz]"); 
plt.ylabel(u"Spectral radiosity [W/m$^2/Hz$]" ); 
plt.grid()
plt.legend(prop={'size':10}, loc='upper right')
plt.savefig("output.png", bbox_inches='tight')