积分黑体光谱以获取太阳的全光度

Question

我尝试把黑体光谱（叫做BBS的函数）整合进来，以计算太阳的总光度（在主函数中叫Lbol），这个值应该大约是3.85乘以10的26次方瓦特。但我得到的结果只有这个值的三分之一。

import numpy as np
from scipy.integrate import quad

global h, c, k # ISU
h = 6.62607e-34
c = 2.998e8
k = 1.38065e-23

global mu_min, mu_max
mu_min, mu_max = 3e10, 3e18
# hertz, corresponds to 1 ångström to 1e8 ångström
# while the sun's spectrum peak at 5000 ångström

global Rsun
Rsun = 6.955e8 # meter

def BBS(mu, tempe):
    i = 2.*h/(c**2.) * (mu**3.) / (np.exp(h/k*mu/tempe)-1.)
    return i

def Teff2Lbol(Teff):
    I = quad(BBS, mu_min, mu_max, args=(Teff,))[0]
    return I

def main():
    T = 5800 # Kelvin
    Lbol = Teff2Lbol(T) * (4*np.pi*Rsun**2.)

Answer 1

你的代码是对的，但物理知识显然不太准确。光谱辐射强度的单位是瓦特每平方米每赫兹每球面度，简写为 W m^-2 Hz^-1 sr^-1。这里的 sr^-1 是因为辐射强度是按单位立体角来计算的，你需要对覆盖发射点的整个半球进行积分。在计算这个积分时，要记住黑体是兰伯特型的，也就是说它们的辐射是遵循余弦定律的：I(theta) = I0*cos(theta)，其中 theta 是表面法线和辐射方向之间的角度。

要得到每单位表面积的总辐射强度，你需要将频率上的积分乘以上半球的 cos(theta) dOmega 的积分（在球坐标系中）。这个积分可以很容易地用解析方法计算，其值恰好是 pi。因此，你需要重新定义 Teff2Lbol 为：

def Teff2Lbol(Teff):
    I = quad(BBS, mu_min, mu_max, args=(Teff,))[0]
    return np.pi*I

另外，请注意，当你在频率范围内积分时，实际上是跨越了8个数量级，而98%的辐射能量实际上是在10¹³ Hz到10¹⁵ Hz之间。幸运的是，QUADPACK 是一个非常好的积分工具，能够处理这种情况。