This notebook demonstrates three different variations of Wallace and Hobbs equation 4.10 All three are contained in the planck.py module, which is in the lib folder in the classcode repository.
The following cell adds the python library directory, lib, to your sys.path so that python will be able to import it's modules. It assumes that you have kept the same folder structure as my repository, i.e. that you are in a folder called "notebooks" and that there is a folder at the same level called "lib"
In [0]:
from __future__ import print_function
import os,site
currdir=os.getcwd()
head,tail=os.path.split(currdir)
libdir=os.path.join(head,'lib')
site.addsitedir(libdir)
import sys
print('python will search the following directories in order for modules: {}'.format(sys.path))
the cell below imports planck and then reloads it, in case I've made changes to the library and we need to rebuild library.pyc
In [0]:
import planck
reload(planck)
from planck import planckwavelen, planckwavenum, planckfreq
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
sigma=5.67e-8 #W/m^2/K^4
c=3.e8 #m/s
the next cell tells ipython to put the graphs right below the cell, instead of popping a new window
In [0]:
%matplotlib inline
In [0]:
fig=plt.figure(1,figsize=(9,9))
plt.clf()
ax1=fig.add_subplot(1,1,1)
#check to see how good this approx. is:
wavel=np.arange(1,50.,0.1)*1.e-6 #meters
TempList=[260,270,280,290,300]
TempNames=["{:d} K".format(item) for item in TempList]
the_temps=zip(TempList,TempNames)
for Temp,name in the_temps:
bbr=planckwavelen(wavel,Temp)*1.e-6*np.pi #convert to W/m^2/micron
ax1.plot(wavel*1.e6,bbr,label=name)
ax1.grid(True)
ax1.set_ylabel('$E^*_\lambda\ {(W\,m^{-2}\,\mu m^{-1})}$')
ax1.legend()
ax1.set_title('Planck function $E^*_\lambda\ {(W\,m^{-2}\,\mu m^{-1})}$ for 5 blackbody temperatures')
ax1.set_xlabel('wavelength $\lambda\ {\mu m}$')
#fig.savefig('/home/phil/Dropbox/lecture/planckII.png')
Temp=300 #K
bbr=planckwavelen(wavel,Temp)
dn=np.diff(wavel)
integ=np.sum(bbr*dn[0]*np.pi)
stefan=sigma*Temp**4.
the_string="""integrated bbr as a function of wavelength (microns)
at 300 K: {:8.3f} W/m^2/micron; Stefan Boltzman: {:8.3f} W/m^2/micron"""
print(the_string.format(integ,stefan))
In [0]:
fig=plt.figure(2,figsize=(9,9))
fig.clf()
ax1=fig.add_subplot(111)
TempList=[260.,270,280,290,300,]
wavenum_icm=np.arange(25,2500,20) #in inverse cm
wavenum_im=wavenum_icm*100. #in inverse m
bList=[]
for Temp in TempList:
bbr_wavenum=planckwavenum(wavenum_im,Temp)*np.pi*100 #convert to W/m^2/cm^-1
ax1.plot(wavenum_icm,bbr_wavenum)
ax1.set_ylabel('$E^*_n\ {(W\,(m^{-2}\,cm^{-1})}$')
ax1.grid(b='off',linewidth=1,linestyle='-',which='both')
ax1.grid(b=False)
ax1.figure.canvas.draw()
xminorLocator = matplotlib.ticker.MaxNLocator(nbins=25)
ax1.xaxis.set_minor_locator(xminorLocator)
yminorLocator = matplotlib.ticker.MaxNLocator(nbins=25)
ax1.yaxis.set_minor_locator(yminorLocator)
ax1.legend(('260 K','270 K','280 K','290 K','300 K'))
ax1.set_title('Planck function $E^*_n\ {(W\,/(m^{2}\,cm^{-1})}$ for 5 blackbody temperatures')
ax1.set_xlabel('wavenumber $n\ {cm^{-1}}$')
#fig.savefig('/home/phil/Dropbox/lecture/q3_retrievalIIb.png',dpi=150)
Temp=300.
bbr=planckwavenum(wavenum_im,Temp)
dn=np.diff(wavenum_im)
integ=np.sum(bbr*dn[0]*np.pi)
stefan=sigma*Temp**4.
print("integrated bbr as a function of wavenumber (m^{{-1}}) at 300 K: {:8.3f} Stefan Boltzman: {:8.3f}".format(integ,stefan))
freq=wavenum_im*c
bbr=planckfreq(freq,Temp)
df=np.diff(freq)
df=df[0]
integ=np.sum(bbr*df*np.pi)
print("integrated bbr as a function of frequency (Hz) at 300 K: {:8.3f} Stefan Boltzman: {:8.3f}".format(integ,stefan))
In [0]:
fig=plt.figure(4,figsize=(10,10))
fig.clf()
bbr_wavenum=planckwavenum(wavenum_im,300.)*100.*np.pi #convert to W/m^2/cm^{-1}
ax1=fig.add_subplot(1,1,1)
fig.subplots_adjust(top=0.8)
first_fill=ax1.fill_between(wavenum_icm,0,bbr_wavenum) #miliWatts/m^2/cm^-1
xlab=ax1.set_xlabel('$wavenumber\ (cm^{-1})$')
xlab.set_fontsize(16)
theTextTitle=plt.figtext(0.5,0.9,'300 K blackbody irradiance $(W/\,(m^{2}\,cm^{-1})$', ha='center')
theTextTitle.set_fontsize(25)
ylab=ax1.set_ylabel('$E_n^*\ (W/(m^{2}\,cm^{-1})$')
ylab.set_fontsize(16)
theText=ax1.text(1500,0.2,'$\sigma T^4 = 460\ W\,m^{-2}$')
theText.set_fontsize(20)
fig.canvas.draw()
#savefig('fig41.png',dpi=300)
sb=5.67e-8*300**4.
level=sb/2500.
ax1.plot([0,2500],[level,level],'r-')
theFill=ax1.fill_between(wavenum_icm,0,level,j
facecolor='red')
theFill.set_alpha(0.2)
#fig.savefig('fig42.png',dpi=300)
ax2 = fig.add_axes(ax1.get_position(), frameon=False)
ax2.set_ylim(ax1.get_ylim())
ax2.set_xlim(ax1.get_xlim())
topticks=np.asarray([1,2,3,4,5,7,10.,12,16])
tickLabels=["%d" % i for i in topticks]
topticks=topticks*1.e-4
ticloc=1./topticks
ax1.xaxis.set_ticks_position('bottom')
ax2.xaxis.set_ticks(ticloc)
ax2.xaxis.set_ticklabels(tickLabels)
ax2.xaxis.tick_top()
ax2.yaxis.tick_right()
ax2.set_xlim(ax1.get_xlim())
ax2.xaxis.set_label_position('top')
xlab2=ax2.set_xlabel('$wavelength\ (\mu m)$')
xlab2.set_fontsize(16)
fig.canvas.draw()
#fig.savefig('fig43.png',dpi=300)
ax1.fill_between(wavenum_icm,0,bbr_wavenum, where=np.logical_and(wavenum_icm > 600,wavenum_icm < 800),
facecolor='green')
#savefig('fig44.png',dpi=300)
In [0]: