from pylab import *

# Planck function to calculate radiance for given wavelength and temperature. 
# Note the difference (factor pi) to task 1 where irradiances were calculated.
def planck(lam, T):
 
    c= 2.99792458e8 # speed of light
    k= 1.380662e-23 # Boltzmann constant
    h= 6.626180e-34 # Planck constant
 
    return 2*h*c*c/(lam**5 * (exp(h*c/(lam*k*T))-1.))


# simulated radiance
bt=loadtxt('radiance.dat')
# atmospheric profile
atm=loadtxt('/local/emde/libRadtran/data/atmmod/afglus.dat')

plot(bt[:,2], bt[:,1], 'o', label='simulated radiance')
plot(planck(10e-6,atm[:,2])*1e-9, atm[:,0], label='B( T(z) )')
ylim(0,10)
xlim(0,0.015)
ylabel('cloud top height [km]')
xlabel('radiance [W/(m^2 nm sr)]')
legend()

savefig('cloudheight.png')