python，绘制通过数组循环的普朗克曲线

Question

我试图让自己熟悉 python 编程，但刚刚开始并努力解决以下问题。也许有人可以给我一个提示，或者我可以在哪里寻找一个好的解决方案。

我想通过一个循环在 6 个不同温度下绘制 132 个波长的普朗克曲线。函数 planckwavel 接收两个参数，波长和温度，我将它们分开在两个循环中。

到目前为止，我设法使用了列表，这很有效，但可能没有以优雅的方式解决：

plancks = []
temp = [280, 300, 320, 340, 360, 380]
temp_len = len(temp)

### via fun planckwavel
for i in range(temp_len):
    t_list = [] # list nach jeder j schleife wieder leeren
    for j in range(wl_centers_ar.shape[0]):
        t = planckwavel(wl_centers_ar[j],temp[i])
        t_list.append(t)      
    plancks.append(t_list)


### PLOT Planck curves

plancks = np.array(plancks).T # convert list to array and transpose

view_7 = plt.figure(figsize=(8.5, 4.5))
plt.plot(wl_centers_ar,plancks)
plt.xticks(rotation='vertical')

但我想使用包含列表的数组，因为我想在之后继续使用更大尺寸的图像。所以我对数组进行了同样的尝试，但不幸的是，这段代码失败了：

plancks_ar = zeros([132,6], dtype=float ) # create array and fill with zeros
temp_ar = array([273, 300, 310, 320, 350, 373])

for i in range(temp_ar.shape[0]):
    t_ar = np.zeros(plancks_ar.shape[0])
    for j in range(plancks_ar.shape[0]):
        t = planck(wl_centers_ar[j]*1e-6,temp[1])/10**6
        np.append(t_ar,t)        
    np.append(plancks_ar, t_ar)

plt.plot(wl_centers_ar,plancks)

如果有人能给我一些建议，我将非常感激。

谢谢， 最好的问候，

彼得

Answer 1

我想你是在问如何使用 NumPy 的 broadcasting and vectorization。这是一种删除显式 Python 循环的方法：

import numpy as np

# Some physical constants we'll need
h, kB, c = 6.626e-34, 1.381e-23, 2.998e8

def planck(lam, T):
    # The Planck function, using NumPy vectorization
    return 2*h*c**2/lam**5 / (np.exp(h*c/lam/kB/T) - 1)

# wavelength array, 3 - 300 um
lam = np.linspace(3, 75, 132)
# temperature array
T = np.array([280, 300, 320, 340, 360, 380])
# Remember to convert wavelength from um to m
pfuncs = planck(lam * 1.e-6, T[:,None])

import pylab
for pfunc in pfuncs:
    pylab.plot(lam, pfunc)
pylab.show()

我们想要为每个波长和每个 T 计算 planck，因此我们需要在两个阵列上广播计算。按照上面链接的文档中列出的规则，我们可以通过向温度数组添加一个新轴来做到这一点（使用T[:, None]）：

lam:       132
T      6 x   1
--------------
       6 x 132

T[:, None] 的最终维度是 1，因此 lam 的 132 个值可以通过它广播以生成一个 6 x 132 数组：6 行（每个 T 一个），包含 132 个值（波长）。

Answer 2

我尝试使用倒数（亮度温度）再次检查上述普朗克方程。以您的代码为导向，我定义了以下函数，并期望获得 300 开尔文（@ 10 微米，对于 10 W/m2/str/microns）：

def planckInv(lam, rad): 
     rad=rad*1.e6 #convert to W/m^2/m/sr 
     lam=lam*1.e-6 #convert wavelength to m 
     return (h*c/kB*lam)*( 1/ np.log( (2*h*c**2/lam**5) / rad +1 )) 

但收到了一个奇怪的结果

planckInv(10,10) - - > 3.0039933569668916e-08

任何建议我的亮度温度功能有什么问题？

谢谢， 彼得