Modify the cells below so that they use lib/planck.planckwavenum to find the wavenumber nmax at which the planck function is maximum. Note: the peak wavenumber is not equal to 1/(peak wavelength) !!!


In [0]:
import numpy as np
from matplotlib import pyplot as plt
from __future__ import print_statement

In [0]:
%matplotlib inline

In [0]:
c=2.99792458e+08  #m/s -- speed of light in vacumn
h=6.62606876e-34  #J s  -- Planck's constant
kb=1.3806503e-23  # J/K  -- Boltzman's constant

def planckwavelen(wavel,Temp):
    """input wavelength in microns and Temp in K, output
    bbr in W/m^2/micron
    """
    wavel=wavel*1.e-6  #convert to meters
    c1=2.*h*c**2.
    c2=h*c/kb
    Elambda=1.e-6*np.pi*c1/(wavel**5.*(np.exp(c2/(wavel*Temp)) -1))
    return Elambda

In [0]:
the_wavelengths=np.linspace(0.1,60,5000) #microns
the_temp=300 #K
flux=planckwavelen(the_wavelengths,the_temp)
fig1=plt.figure(1)
axis1=fig1.add_subplot(1,1,1)
axis1.plot(the_wavelengths,flux)

In [0]:
dlambda=np.diff(the_wavelengths)
dflux=np.diff(flux)
the_deriv=dflux/dlambda
fig2,axis2=plt.subplots(1,1)
axis2.plot(the_wavelengths[:-1],the_deriv,label="approx derivative")
wiens_law=2897/the_temp
axis2.plot(wiens_law,0,'k+',markersize=20,label="Wien's law")
axis2.legend()
axis2.set_title("Wien's law demonstration")
axis2.set_ylabel("$dE_lambda/dlambda\ (W/m^2/micron/micron)$")
axis2.set_xlabel('wavelength lambda (microns)')
#savefig('wien_answer.png')

Borrow code from plotplanck.ipynb


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)
from planck import planckwavenum

In [0]:
wavenum_icm=np.arange(25,2500,20) #in inverse cm
wavenum_im=wavenum_icm*100.  #in inverse m
flux=planckwavenum(wavenum_im,the_temp)*np.pi*100 #convert to W/m^2/cm^-1

In [0]:
fig=plt.figure(2,figsize=(9,9))
fig.clf()
ax1=fig.add_subplot(111)
ax1.plot(wavenum_icm,flux)

Borrow code from the wien's law wavelength solution


In [0]:
dlambda=np.diff(wavenum_icm)
dflux=np.diff(flux)
the_deriv=dflux/dlambda
zero_crossings=np.where(np.diff(np.sign(the_deriv)))
num_max=wavenum_icm[zero_crossings[0]][0]
print("found a maximum at {} cm-1".format(num_max))

In [0]:
wiens_law_wavelength=2897./the_temp
wiens_law_inverse_microns=1./num_max*1.e4 #convert to wavelength in microns
print_tup=(wiens_law_wavelength,num_max,wiens_law_inverse_microns)
text="""note that for this temperature Wien's law for wavelength says that the
        maximum occurs at {:6.3f} microns but taking the inverse of the wavenumber maximum
        shows that this peak is located at 1/{} cm = {:6.3f} microns""".format(*print_tup)
print(text)

In [0]: