黑体光谱的普朗克公式答案

【问题标题】：Plancks Formula for Blackbody spectrum黑体光谱的普朗克公式
【发布时间】：2014-04-20 11:32:56
【问题描述】：

我正在尝试为给定温度 T=200K 的强度与波长图编写一个简单的 Python 代码。到目前为止，我有这个......

import scipy as sp
import math
import matplotlib.pyplot as plt
import numpy as np
pi = np.pi
h = 6.626e-34
c = 3.0e+8
k = 1.38e-23

def planck(wav, T):
    a = 2.0*h*pi*c**2
    b = h*c/(wav*k*T)
    intensity = a/ ( (wav**5)*(math.e**b - 1.0) )
    return intensity

我不知道如何定义波长（wav），从而产生普朗克公式图。任何帮助将不胜感激。

【问题讨论】：

大概你想要一个波长范围。你知道什么范围吗？查看range。
通常波长只是偶数间隔矢量，可以使用np.arange 或np.linspace。记住普朗克函数中的缩进。

标签： python numpy matplotlib scipy

【解决方案1】：

您可能还想使用RADIS 库，它允许您根据波长或频率/波数绘制普朗克函数（如果需要）！

from radis import sPlanck

sPlanck(wavelength_min=135, wavelength_max=3000, T=4000).plot()
sPlanck(wavelength_min=135, wavelength_max=3000, T=5000).plot(nfig='same')
sPlanck(wavelength_min=135, wavelength_max=3000, T=6000).plot(nfig='same')
sPlanck(wavelength_min=135, wavelength_max=3000, T=7000).plot(nfig='same')

【讨论】：

【解决方案2】：

这是一个基本情节。要使用plt.plot(x, y, fmt) 绘图，您需要两个大小相同的数组 x 和 y，其中 x 是要绘制的每个点的 x 坐标，y 是 y 坐标，fmt 是描述如何绘制数字的字符串。

因此，您需要做的就是创建一个均匀分布的波长阵列（我将其命名为 np.arraywavelengths）。这可以通过arange(start, end, spacing) 来完成，它将创建一个从头到尾（不包括在内）以spacing 分隔的数组。

然后使用您的函数计算数组中每个点的强度（将存储在另一个np.array 中），然后调用plt.plot 绘制它们。注意 numpy 让您在 vectorized form 中快速对数组进行数学运算，这将是计算高效的。

import matplotlib.pyplot as plt
import numpy as np

h = 6.626e-34
c = 3.0e+8
k = 1.38e-23

def planck(wav, T):
    a = 2.0*h*c**2
    b = h*c/(wav*k*T)
    intensity = a/ ( (wav**5) * (np.exp(b) - 1.0) )
    return intensity

# generate x-axis in increments from 1nm to 3 micrometer in 1 nm increments
# starting at 1 nm to avoid wav = 0, which would result in division by zero.
wavelengths = np.arange(1e-9, 3e-6, 1e-9) 

# intensity at 4000K, 5000K, 6000K, 7000K
intensity4000 = planck(wavelengths, 4000.)
intensity5000 = planck(wavelengths, 5000.)
intensity6000 = planck(wavelengths, 6000.)
intensity7000 = planck(wavelengths, 7000.)


plt.plot(wavelengths*1e9, intensity4000, 'r-') 
# plot intensity4000 versus wavelength in nm as a red line
plt.plot(wavelengths*1e9, intensity5000, 'g-') # 5000K green line
plt.plot(wavelengths*1e9, intensity6000, 'b-') # 6000K blue line
plt.plot(wavelengths*1e9, intensity7000, 'k-') # 7000K black line

# show the plot
plt.show()

你会看到：

您可能需要清理坐标轴标签、添加图例、在同一图上绘制多个温度下的强度，等等。请咨询relevant matplotlib documentation。

【讨论】：

为什么是math.e**b 而不仅仅是np.exp(b)
因为我从user3421850的问题中获取了函数，而他们使用了math.e**b。
@Zhenya - 现在我看了一下，他们的公式不正确，因为他们在a 的公式中有一个额外的 pi。（注意 h 是普朗克常数，而不是 h-bar（普朗克约化常数 ħ = h/2pi）。）
Re 单位：总是使用无量纲单位而不是（任意）SI 单位要好。将公式写在一张纸上，使用替换（这里，x = \hbar * \omega / k T），用无量纲变量的无量纲函数表示 rhs 和 lhs，绘制它们（并使用轴标签来说明确切的变量）。
@DevonDahon - Matplotlib 2.1 已于 2017 年 10 月弃用。此评论早于该评论三年。虽然会编辑。我相信只需要删除对hold() 的调用，因为这是默认行为。

【解决方案3】：

只是想指出，OP 想要在 astropy 中做的事情似乎与此等价：

https://docs.astropy.org/en/stable/api/astropy.modeling.physical_models.BlackBody.html

不幸的是，我还不太清楚如何获得基于波长与频率的表达式。

【讨论】：