12. Plotting the Planck function - solution#

import numpy as np
from matplotlib import pyplot as plt

12.1. Write a function to compute Stull 2.13#

import numpy as np
#
# get Stull's c_1 and c_2 from fundamental constants
#
# c=2.99792458e+08  #m/s -- speed of light in vacuum
# h=6.62606876e-34  #J s  -- Planck's constant
# k=1.3806503e-23  # J/K  -- Boltzman's constant

def Flambda(wavel, Temp):
    """
    Calculate the blackbody radiant exitence (Stull 2.13)

    Parameters
    ----------

      wavel: float or array
           wavelength (meters)

      Temp: float
           temperature (K)

    Returns
    -------

    Flambda:  float or arr
           monochromatic radiant exitence (W/m^2/m)
    """
    c, h, k = 299_792_458.0, 6.626_070_04e-34, 1.380_648_52e-23
    c1 = 2.0 * h * c ** 2.0
    c2 = h * c / k
    Flambda_val = c1 * np.pi / (wavel ** 5.0 * (np.exp(c2 / (wavel * Temp)) - 1))
    return Flambda_val

12.2. run the function for a single temperature#

npoints = 10000
Temp = 255  # K
wavelengths = np.linspace(0.1, 500.0, npoints) * 1.0e-6  # meters
Fstar = Flambda(wavelengths, Temp)
fig, ax = plt.subplots(1, 1, figsize=(7, 7))
#
#  change wavelength to microns and flux to W/m^2/micron
#
ax.plot(wavelengths * 1.0e6, Fstar * 1.0e-6)
ax.set(xlim=[0, 50])
ax.grid(True)
ax.set(
    xlabel="wavelength (m)",
    ylabel=r"$F_\lambda^*\ (W\,m^{-2}\,\mu m^{-1}$)",
    title=f"Monochromatic blackbody flux at Temp={Temp} K",
);

12.3. Convert flux to radiance#

This uses the reading Finding the flux given the radiance

# convert isotropic flux to radiance
Blambda = Fstar / np.pi
fig, ax = plt.subplots(1, 1, figsize=(8, 8))
#convert wavelength from meters to microns
ax.plot(wavelengths * 1.0e6, Blambda * 1.0e-6)
ax.set(xlim=[0, 50])
ax.grid(True)
ax.set(
    xlabel="wavelength ($\mu m$)",
    ylabel="$B_\lambda\ (W\,m^{-2}\,sr^{-1}\,\mu m^{-1}$)",
    title=f"Monochromatic blackbody radiance at Temp={Temp} K",
);

12.4. Reproduce W&H figure 4.6#

In the cell below add lines and change the axis limits to reproduce the high temperature emission spectrum in the W&H figure

temps =[5000, 6000, 7000] #K
# convert isotropic flux to radiance
npoints=1000
wavelengths = np.linspace(0.01,2,npoints)*1.e-6  #meters
fig, ax = plt.subplots(1, 1, figsize=(8, 8))
for the_temp in temps:
    Fstar = Flambda(wavelengths, the_temp)
    Blambda = Fstar / np.pi
    #convert wavelength from meters to microns
    #convert radiance to megawatts per meter^2 per micron bin width per sr
    ax.plot(wavelengths * 1.0e6, Blambda*1.e-12,label=f"{the_temp} K")
    ax.set(xlim=[0., 2])
ax.grid(True)
ax.legend()
ax.set(
    xlabel="wavelength ($\mu m$)",
    ylabel="$B_\lambda\ (MW\,m^{-2}\,sr^{-1}\,\mu m^{-1}$)",
    title=f"Monochromatic blackbody radiance vs temperature"
);