Cookbook / 物理 / 06
60 分钟参考常数与模型已执行的网页阅读版
下载并在 VibeIt 中运行所有教程

核对普朗克、维恩与斯忒藩—玻尔兹曼定律

从冻结 SI 常数计算光谱辐亮度，定位波长谱峰，并核对积分辐射出射度。

物理 · 第 06 课🔋 离线核心分析⏱ 60 分钟参考常数与模型📦 numpy · pandas · scipy · matplotlib · IPython

目标与运行方法

本课适合具备 Python 基础的本科生。先阅读每一步的方法与图注，再修改参数。AI 练习为可选部分，不需要 AI 账号即可完成核心课程。数据写在 notebook metadata 中，不需要另下载数据文件。请保持 notebook 已保存到当前工作文件夹。

学习任务是理解本课的研究问题、执行透明的计算、校验结果，并说明结论的边界。

环境与参数

In [1] · Python
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

读取内嵌快照

In [2] · Python
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']

方法函数

In [3] · Python
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

数据与模型边界

单摆、振子、抛体和温度轨迹是明确声明的数值模型，不是实验室测量。使用的 g0 为约定标准重力加速度。GW150914 使用真实冻结应变和已发表事件时间；滤波曲线不是独立的显著性或源参数分析。黑体曲线假定理想辐射体。

分步分析

1. 以 SI 单位计算波长谱辐亮度

使用冻结的 h、c、k_B。B_lambda 按每米波长计，乘 10⁻⁹ 才换为每 nm。波长形式与频率形式需要 Jacobian，谱峰坐标不同。

In [4] · Python
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()
Computed figure from this code cell

2. 定位连续波长极大值

有界标量优化器以米为单位搜索波长，容差属于数值设置。用波长形式的维恩常数核对 lambda_peak·T；仅选网格最高点会增加位置离散误差。

In [5] · Python
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()
Computed figure from this code cell

3. 将辐亮度积分为出射度

理想各向同性表面的出射度为波长积分辐亮度乘 pi。无量纲变量转换给出 x³/(exp(x)−1) 的积分，等于 pi⁴/15，可避免无限波长区间的尺度问题。

In [6] · Python
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)
Computed figure from this code cell
Derived sigma (W m⁻² K⁻⁴): 5.6703744191844314e-08 ; dimensionless integral: 6.493939402266821

4. 核对正值、谱峰和积分单位

核对独立积分恒等式及 T⁴ 指数。这些比较验证理想定律和数值实现，不会替真实材料验证发射率、温度或探测器校准。

In [7] · Python
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

自主练习与参考答案

  1. 为何绘图前要把每米辐亮度换为每 nm？
  2. 为何波长和频率表示的谱峰不同？

参考解释。 谱密度随坐标微分转换，波长/频率 Jacobian 改变密度及极大值位置。

先保留当前 notebook 副本，再修改参数。把新图与原图一起比较，并记录改变了哪些输入、哪些计算和哪些结论。

可选的 AI 编程练习

提示词 1

仅用 NumPy、SciPy、Matplotlib 加入理想 6500 K 光谱及谱峰，不改变冻结常数。

提示词 2

使用同样的包实现频率形式普朗克辐亮度，注明 W m^-2 sr^-1 Hz^-1 并解释谱峰差别。

人工验收。 光谱为正，单位换算明确，积分恒等式及 sigma T^4 校验通过。

检查修改后的源代码，再手动运行受影响单元及后续校验。要求助手解释修改的目的、使用的包和数据来源，并用实际输出核对解释。

可选联网扩展

In [8] · Python
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

来源、快照与下一步

  • 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.

下一步。 在保留本课的数据检查、参数记录和结果校验之后，将相同方法用于自己的公开或已获授权数据。记录研究设计与方法限制，不把示例结果当成新实验的验证。完整数据许可与生成说明见 cookbook 网页。