使用Python基于Planck方程绘制强度-频率双对数图遇错求助

错误原因拆解 #

• np.logspace参数错误：
  • logspace(start, stop, num)的start和stop是10的指数值，不是实际频率值。你要生成10^9到1019 Hz的频率，应该传入start=9、stop=19，而非10e+9和10e+19。
  • num=1e10是浮点数，而该参数要求必须是整数；同时10亿个计算点会直接耗尽内存，完全没必要，改成1000个点足够绘图。
• Y轴范围设置不合理：你设置的plt.ylim([1,1e10])和Planck公式计算出的实际强度值差距极大，会导致图像显示异常。

修复后的完整代码 #

import matplotlib.pyplot as plt 
import numpy as np

# 物理常数
h = 6.626e-34  # 普朗克常数
c = 3.0e+8     # 光速
k = 1.38e-23   # 玻尔兹曼常数

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

# 生成频率范围：10^9 ~ 10^19 Hz，1000个对数间隔点
frequency = np.logspace(9, 19, 1000, endpoint=False) 

# 计算不同温度下的强度
intensities = {
    3000: planck(frequency, 3000.),
    5000: planck(frequency, 5000.),
    10000: planck(frequency, 10000.),
    20000: planck(frequency, 20000.),
    30000: planck(frequency, 30000.)
}

# 绘制双对数曲线
plt.loglog(frequency, intensities[3000], 'k-', label='3000K')
plt.loglog(frequency, intensities[5000], 'y-', label='5000K')
plt.loglog(frequency, intensities[10000], 'r-', label='10000K')
plt.loglog(frequency, intensities[20000], 'g-', label='20000K')
plt.loglog(frequency, intensities[30000], 'b-', label='30000K')

# 设置坐标轴范围和标签
plt.xlim([1e9, 1e19])
plt.ylim([1e-20, 1e-5])
plt.xlabel('频率 (Hz)')
plt.ylabel('辐射强度 (Joule/sec/m²/steradian/Hertz)')
plt.legend()

# 显示图像
plt.show()

额外说明 #

• 用字典存储不同温度的强度值，比单独定义变量更简洁，也方便后续扩展。
• 添加了坐标轴标签和图例，让图像可读性更强。
• 若遇到高温高频下np.exp(b)溢出的情况，可添加判断：当b > 700时，直接用np.exp(b)代替np.exp(b)-1（此时两者数值几乎相等），避免计算错误。

内容的提问来源于stack exchange，提问作者VISMAYA P V