from pylab import *

def planck(lam, T):

    c= 2.99792458e8 # speed of light
    k= 1.380662e-23 # Boltzmann constant
    h= 6.626180e-34 # Planck constant
   
    return 2*pi*h*c*c/(lam**5 * (exp(h*c/(lam*k*T))-1.))

lam=arange(0.05, 1000, 0.01)*1e-6

d= 1.49e11      # distance between sun and earth
r= 1.3914e9/2.  # radius of sun

figure(1, figsize=(8,5))
loglog(lam*1e6, planck(lam, 6000)*1e-6, label='6000 K')
loglog(lam*1e6, planck(lam, 6000)*(r/d)**2*1e-6, ':' )
loglog(lam*1e6, planck(lam, 300)*1e-6, label='300 K')
ylim(1e-1,1e9)
legend()
xlabel('wavelength [$\mu$m]')
ylabel('irradiance [W/(m$^2$ $\mu$m]')

savefig('planck_functions.png')