from pylab import *

def planck(lam, T):
    
    c= 2.99792458e8
    k= 1.380662e-23
    h= 6.626180e-34
    return 2*pi*h*c*c/(lam**5 * (exp(h*c/(lam*k*T))-1.))

# irradiance from 300 to 500 nm
lw_spectrum_sur=loadtxt('thermal_sur.out')
lw_spectrum_toa=loadtxt('thermal_toa.out')

figure(1, figsize=(10,7))

plot(lw_spectrum_sur[:,0], lw_spectrum_sur[:,3], label='surface')
plot(lw_spectrum_toa[:,0], lw_spectrum_toa[:,3], label='TOA')
plot(lw_spectrum_toa[:,0], planck(lw_spectrum_toa[:,0]*1e-9, 288.2)*1e-9, label='Planck')
xlim(3500,60000)
title('LW spectrum')
xlabel('wavelength [nm]')
ylabel('irradiance [W/ (m^2 nm)]')
legend(loc=1)

savefig('olr.png')