How to calculate and plot Planck's law with python ?

Planck's law describes the spectral energy distribution of electromagnetic radiation emitted by a black body in thermal equilibrium. Calculating and plotting Planck's law can be done with python using various methods and libraries.

Get fundamental physical constants

To calculate the radiance emitted from temperature using Planck's principle, we need to incorporate fundamental physical constants such as the speed of light, Boltzmann Constant and Plank Constant. To easily access physical constants in your Python code, a great solution is to use Scipy and its import constants feature.

from scipy import constants

c = constants.speed_of_light
h = constants.h
k = constants.k

Create a Planck function

Now, let's craft a python Planck function which will take in both wavelength and temperature as inputs and return the radiance:

k1 = 2.0*h*c**2
k2 = h*c/k

def planck_function(wl, T, k1, k2):
    radiance = k1 / ( (wl**5) * (np.exp(k2/(wl*T)) - 1.0) )
    return radiance

k1 and k2 are constants that remain unchanged, so there is no need to compute them each time the function runs.

Plot Planck function

The last task is to utilize Matplotlib for visually representing Planck's function across a range of temperatures and wavelengths.

import matplotlib.pyplot as plt
import numpy as np

wls = np.linspace(1e-9, 3e-6, 1000)

temps = [3000,4000,5000,6000]

for temp in temps:
    plt.plot(wls, planck_function(wls, temp, k1, k2))

plt.grid('..')

plt.title("How to calculate and plot Planck's law with python ?")

plt.savefig('plank_law_01.png', dpi=100, facecolor='white')

plt.show()

Ouput

Aggregated code

from scipy import constants

import matplotlib.pyplot as plt
import numpy as np

c = constants.speed_of_light
h = constants.h
k = constants.k

k1 = 2.0*h*c**2
k2 = h*c/k

def planck_function(wl, T, k1, k2):
    radiance = k1 / ( (wl**5) * (np.exp(k2/(wl*T)) - 1.0) )
    return radiance

wls = np.linspace(1e-9, 3e-6, 1000)

temps = [3000,4000,5000,6000]

for temp in temps:
    plt.plot(wls, planck_function(wls, temp, k1, k2))

plt.grid('..')

plt.title("How to calculate and plot Planck's law with python ?")

plt.savefig('plank_law_01.png', dpi=100, facecolor='white')

plt.show()

References