Reconcile Planck, Wien and Stefan–Boltzmann laws
Compute spectral radiance from frozen SI constants, locate wavelength peaks and verify integrated exitance.
Goal and how to run
This lesson assumes basic Python knowledge. Read each method and caption before changing parameters. AI exercises are optional; no AI account is needed for the core workflow. Data lives in notebook metadata, so there is no companion data download. Keep the notebook saved in the active working folder.
Your goal is to understand the research question, execute transparent calculations, check the results, and state the limits of the conclusion.
Setup and parameters
from pathlib import Path
import base64, gzip, hashlib, json, platform
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from IPython.display import display, HTML
from scipy.optimize import minimize_scalar
from scipy.integrate import quad
from IPython.display import display as _display
def display(value):
if isinstance(value,pd.DataFrame):
table=value.to_html(max_rows=12,escape=True)
style="<style>.vibeit-readable-table{max-width:100%;overflow-x:auto}.vibeit-readable-table table{width:max-content;max-width:none;border-collapse:collapse}.vibeit-readable-table td,.vibeit-readable-table th{white-space:nowrap!important;word-break:normal!important;overflow-wrap:normal!important}</style>"
return _display(HTML(style+'<div class="vibeit-readable-table">'+table+'</div>'))
return _display(value)
RECIPE_ID='ph06-blackbody'
OUTPUT_DIR=Path.cwd()/"results"/RECIPE_ID
plt.rcParams.update({"figure.dpi":120,"font.size":10,"axes.spines.top":False,"axes.spines.right":False,
"axes.prop_cycle":plt.cycler(color=["#18785f","#416aa6","#bc6541","#9974a5"])})
print("Python:",platform.python_version(),"; NumPy:",np.__version__,"; Pandas:",pd.__version__)
print("Working folder:",Path.cwd())
print("Core data mode: embedded snapshot; no network requests")
Python: 3.13.13 ; NumPy: 2.5.3 ; Pandas: 3.0.6 Working folder: /Users/zhiluo/Developer/vibeit-cookbook/cookbook/downloads Core data mode: embedded snapshot; no network requests
Load the embedded snapshot
def load_snapshot(recipe_id, expected_sha256, snapshot_key):
"""Read embedded data from this notebook, including a renamed copy."""
for path in sorted(Path.cwd().glob("*.ipynb")):
try:
notebook = json.loads(path.read_text(encoding="utf-8"))
metadata = notebook.get("metadata", {}).get("vibeit_cookbook", {})
blob = metadata.get("snapshots", {}).get(snapshot_key)
if metadata.get("recipe_id") != recipe_id or not blob:
continue
if blob.get("sha256") != expected_sha256:
continue
raw = gzip.decompress(base64.b64decode(blob["gzip_base64"], validate=True))
if hashlib.sha256(raw).hexdigest() != expected_sha256:
raise ValueError("Embedded snapshot checksum mismatch: " + snapshot_key)
return json.loads(raw)
except (OSError, json.JSONDecodeError):
continue
raise FileNotFoundError(
"Save/open the complete .ipynb in the active working folder. "
"The embedded snapshot could not be located for " + recipe_id)
SNAPSHOT_HASHES={'constants': '2533802a1a546c612d06ba2e931e4f32a418be854c75a5e9c327ff57d5b56126'}
def snapshot(key):
return load_snapshot(RECIPE_ID,SNAPSHOT_HASHES[key],key)
print("Embedded snapshots:",list(SNAPSHOT_HASHES))
Embedded snapshots: ['constants']
Read the method functions
def planck_radiance(wavelength_m,temperature_K,h,c,k_B):
"""Spectral radiance per metre of wavelength: W m^-2 sr^-1 m^-1."""
wavelength=np.asarray(wavelength_m,dtype=float)
if not np.isfinite(wavelength).all() or np.any(wavelength<=0) or temperature_K<=0 or min(h,c,k_B)<=0:
raise ValueError('Planck law requires positive finite wavelengths, temperature and constants')
exponent=h*c/(wavelength*k_B*temperature_K);values=np.zeros_like(wavelength)
safe=exponent<700
values[safe]=2*h*c*c/(wavelength[safe]**5*np.expm1(exponent[safe]))
return values
print("Method ready: planck_radiance")
Method ready: planck_radiance
Data and model boundary
Pendulum, oscillator, projectile and heat trajectories are declared numerical models, not laboratory measurements. The gravitational acceleration is the conventional standard g0 when used. GW150914 uses actual frozen detector strain and published event timing; a filtered curve is not an independent significance or source-parameter analysis. Blackbody curves assume an ideal emitter.
Steps
1. Compute wavelength spectral radiance in SI units
Use frozen h, c and k_B. B_lambda is per metre of wavelength; multiplying by 10⁻⁹ expresses it per nm. The wavelength and frequency versions require a Jacobian and have different peak coordinates.
constants=snapshot("constants")["constants"]
h=constants["h"]["value"];c=constants["c"]["value"];k_B=constants["k_B"]["value"]
temperatures=[3000,5000,8000];wavelength=np.linspace(100e-9,3000e-9,1001)
fig,ax=plt.subplots(figsize=(9,4))
for temperature in temperatures:ax.plot(1e9*wavelength,1e-9*planck_radiance(wavelength,temperature,h,c,k_B),label=f"{temperature} K")
ax.set(xlabel="Wavelength (nm)",ylabel="Radiance (W m⁻² sr⁻¹ nm⁻¹)",title="Ideal blackbody wavelength spectra");ax.legend();plt.tight_layout();plt.show()
2. Locate continuous wavelength maxima
A bounded scalar optimizer searches wavelengths in metres. Its tolerance is a numerical setting. Confirm lambda_peak·T against the wavelength Wien constant; simply taking the highest sampled grid point adds grid-location error.
peaks=[]
for temperature in temperatures:
optimum=minimize_scalar(lambda wave:-planck_radiance(np.array([wave]),temperature,h,c,k_B)[0],bounds=(50e-9,3000e-9),method="bounded",options={"xatol":1e-15})
assert optimum.success
peaks.append(optimum.x)
peaks=np.array(peaks)
fig,ax=plt.subplots(figsize=(7,4));ax.plot(temperatures,1e9*peaks,"o-",label="Continuous spectral maximum")
ax.plot(temperatures,1e9*2.897771955e-3/np.asarray(temperatures),"--",label="Wavelength Wien reference")
ax.set(xlabel="Model temperature (K)",ylabel="Peak wavelength (nm)",title="Peak-coordinate cross-check");ax.legend();plt.tight_layout();plt.show()
3. Integrate radiance into exitance
For an ideal isotropic surface, exitance is pi times wavelength-integrated radiance. A dimensionless change of variable yields the integral of x³/(exp(x)−1), equal to pi⁴/15. This avoids a poorly scaled wavelength quadrature over an infinite interval.
integral,estimated_error=quad(lambda x:x**3/np.expm1(x),0,100,epsabs=1e-10)
sigma=2*np.pi**5*k_B**4/(15*h**3*c**2)
grid_temperature=np.linspace(1000,10000,101)
numeric_exitance=2*np.pi*k_B**4*grid_temperature**4/(h**3*c**2)*integral
reference_exitance=sigma*grid_temperature**4
fig,ax=plt.subplots(figsize=(7,4));ax.loglog(grid_temperature,numeric_exitance,label="Dimensionless quadrature")
ax.loglog(grid_temperature,reference_exitance,"--",label="Derived sigma T⁴")
ax.set(xlabel="Model temperature (K)",ylabel="Exitance (W/m²)",title="Integrated spectral law and T⁴ scaling");ax.legend();plt.tight_layout();plt.show()
print("Derived sigma (W m⁻² K⁻⁴):",sigma,"; dimensionless integral:",integral)
Derived sigma (W m⁻² K⁻⁴): 5.6703744191844314e-08 ; dimensionless integral: 6.493939402266821
4. Reconcile positivity, peaks and integrated units
Verify the independent integral identity and T⁴ exponent. These comparisons concern ideal laws and numerical implementation; they do not validate emissivity, temperature or detector calibration for a real material.
assert np.isclose(integral,np.pi**4/15,rtol=1e-10)
assert np.allclose(peaks*np.asarray(temperatures),2.897771955e-3,rtol=1e-5)
assert np.allclose(numeric_exitance,reference_exitance,rtol=1e-10)
assert np.isclose(np.polyfit(np.log(grid_temperature),np.log(numeric_exitance),1)[0],4)
OUTPUT_DIR.mkdir(parents=True,exist_ok=True)
pd.DataFrame({"wavelength_nm":wavelength*1e9,**{f"radiance_per_nm_{t}K":1e-9*planck_radiance(wavelength,t,h,c,k_B) for t in temperatures}}).to_csv(OUTPUT_DIR/"blackbody-spectra.csv",index=False)
pd.DataFrame({"temperature_K":temperatures,"peak_wavelength_nm":peaks*1e9}).to_csv(OUTPUT_DIR/"spectral-peaks.csv",index=False)
print("Integral, wavelength-peak and scaling checks passed; exports:",OUTPUT_DIR)
Integral, wavelength-peak and scaling checks passed; exports: /Users/zhiluo/Developer/vibeit-cookbook/cookbook/downloads/results/ph06-blackbody
Exercises and reference explanation
- Why convert per-metre radiance to per-nm before plotting?
- Why do wavelength and frequency representations peak differently?
Reference explanation. Spectral density transforms with the coordinate differential. The wavelength/frequency Jacobian changes density and maximum coordinates.
Keep a copy before changing parameters. Compare old and new figures and record which inputs, calculations and conclusions changed.
Optional AI coding exercises
Prompt 1
Using NumPy/SciPy/Matplotlib, add an ideal 6500 K spectrum and peak without changing the frozen constants.
Prompt 2
With the same packages, implement frequency-domain Planck radiance with explicit W m^-2 sr^-1 Hz^-1 units and explain the peak difference.
Human acceptance. Spectra remain positive, unit conversions are stated, the integral identity and sigma T^4 checks pass.
Review the edited source, then manually run affected cells and subsequent checks. Ask the assistant to explain the purpose, packages and data provenance of its changes, and reconcile that explanation with actual output.
Optional online extension
RUN_ONLINE_UPDATE=False
if RUN_ONLINE_UPDATE:
import requests
try:
response=requests.get('https://physics.nist.gov/cuu/Constants/',timeout=20)
response.raise_for_status()
print("Source:",response.url)
print(response.text[:700])
except requests.RequestException as error:
print("Online extension unavailable:",type(error).__name__,str(error)[:160])
else:
print("Online extension skipped; embedded core data unchanged")
Online extension skipped; embedded core data unchanged
Sources, snapshots and next steps
- NIST CODATA physical constants — NIST public reference data; credit NIST/CODATA; retrieved 2026-10-09. SHA-256
77fb90e66c40db3e6eb16630bc9c88e4c7c8beddbe5e71be406f2f26e3f67e67. Preserve official ASCII table; extract seven explicitly named SI constants and their units.
Next steps. Transfer these methods to your own public or authorized data while retaining input checks, recorded parameters and result validation. Document the new experimental design and method limits; this example does not validate a new experiment. Full data licenses and generation instructions accompany the cookbook.