Python 图形学教程（四）

原文：Python Graphics

协议：CC BY-NC-SA 4.0

十、示例

在这一章中，你将应用在前面章节中开发的一些技术来制作一些有趣的图像。这些图片应该让您对 Python 图形可以完成的事情有所了解。

10.1 土星

土星以它的光环而闻名。虽然木星、土星、天王星和海王星也有光环，但土星是太阳系中最大、最亮、最著名的。它们由小到灰尘，大到巨石大小的物体组成。这些物体主要由冰组成，被认为是彗星或大型小行星与土星的一颗卫星相撞时产生的，两者都碎成了碎片。土星自古以来就为人所知，但在 1610 年，伽利略首次用望远镜观察到了它。这颗行星是以罗马农业之神土星命名的，就像我们的第六天，星期六一样。

清单 10-1 建立在更早的程序上，清单 7-2 来自第七章。除了引入构建土星环和土星在环上的阴影的算法之外，这个程序在这里基本保持不变。

图 10-1 到 10-5 显示了列表 10-1 产生的图像。它们处于不同的方位角度，这些都列在标题中。还列出了入射光线的单位矢量分量。比如 lx=+.707，ly=+.707，lz=0 表示左上象限的一个源；lx=-1，ly=0，lz=0 表示光源来自右侧。在图像中，请注意行星投在光环上的阴影，特别是图 10-5 ，它显示了行星的曲率。

作为比较，土星的照片可以在 www.jpl.nasa.gov/spaceimages/?search=saturn&category=#submit 找到。

图 10-5

Saturn with rings and shadow 5: Rx=20, Ry=0, Rz=30, lx=-1, ly=0, lz=0 (produced by Listing 10-1)

图 10-4

Saturn with rings and shadow 4: Rx=-10, ry=0, Rz=25, lx=-.707, ly=-.707, lz=0 (produced by Listing 10-1)

图 10-3

Saturn with rings and shadow 3, Rx=-20, Ry=0, Rz=25, lx=.707, ly=.707, lz=0 (produced by Listing 10-1)

图 10-2

Saturn with rings and shadow 2: Rx=-8, Ry=0, Rz=30, lx=.707, ly=.707, lz=0 (produced by Listing 10-1)

图 10-1

Saturn with rings and shadow 1: Rx=-20, Ry=0, Rz=-10, lx=1, ly=0, lz=0 (produced by Listing 10-1)

图 10-6 显示了用于构建环的模型。在第七章中，你通过首先创建一个直立球体开发了阴影球体算法。也就是说，经度是垂直的，纬度是水平的(即平行于 x，z 平面)。从这个起始方向开始，围绕 x、y 和 z 轴旋转球体。你在这里为戒指做同样的事情。创建平行于 x，z 平面的水平环，然后将它们与球形行星体一起旋转相同的角度。环位于通过球体中心的平面上，因此球体和环具有相同的旋转中心。

图 10-6

Rings model: top view of planet and rings looking down on the x,z plane with Rx=0, Ry=0, Rz=0

该波段被绘制为一系列相邻的同心圆，每个同心圆由短线段组成。参考图 10-6 和清单 10-1 ，程序行 42 和 43 设置环的内半径和外半径。第 44 行设置圆之间的距离。环被分成七个环形带(图 10-6 中未显示)以适应不同的颜色；它们的宽度是第 45 行中的δr。

每条线段都是单独旋转和绘制的。线 48 开始了从 r1 到 r2 的径向循环，绘制了圆段。线 49 开始沿圆周方向绘制的循环。线 50-61 进行旋转，在线 62 和 63 中产生全局标绘坐标 xpg 和 ypg。旋转功能与先前程序中的功能相同。

接下来，设置线段的颜色。这些光环被排列成不同颜色的条带，这是美国宇航局图像中看到的它们的物理组成的结果。这是在第 66-75 行完成的。从 r=r1 到 R1+δr 的第一个波段的颜色为 clr=(.63，. 54，… 18)，其余波段以此类推。你省略了第五个波段，它是空的；背景色清晰可见。第六波段是其他波段的两倍宽。这为七个波段提供了颜色。

对于给定的光线方向，在大多数方向上，行星的身体会在光环上投下阴影。参考图 10-7 ，你的目标是确定 p 点是在行星的阴影区之内还是之外。球形行星投下圆形阴影。阴影的直径将等于行星的大小，或者更准确地说，是球体的“大圆”这是用穿过球体中心的平面切割球体所能得到的最大圆。这就像把一个橘子切成两半；你看到的是橘子的大圆。在图 10-7 中，阴影可能是由这种大小的圆盘造成的，也可能是由球形行星造成的；无论哪种情况，阴影的大小都是一样的。土星大圆的侧视图显示为穿过飞机中心的粗线。从图 10-7 中的几何图形可以看出，如果 p 位于这样一个位置，使得|B| > rs，其中 rs 是土星的半径，它在阴影区之外；如果|B| < rs，p 在阴影区内。一旦你确定了 p 的位置，如果它在阴影区域内，当你画圆环的时候，你将把那个点涂成灰色。如果它在外面，你将在第 66-75 行中给它一个带颜色。

图 10-7

Shadow model

你现在的工作是对于 p 的给定位置得到|B|，你从图 10-7 看到

)

(10-1)

你知道

)

(10-2)其中 uˇ=-^ˇl .结合以上等式用| uˇ| = 1、

)

(10-3)

)

(10-4)

在清单 10-1 中，第 78 行确定了入射光矢量的长度，^ˇl。它应该等于 1，但如果第 23-25 行输入的分量不等于 1(即))，它可能不等于 1。如果需要，第 79-81 行重新建立组件。第 82-84 行建立了矢量 v 的分量。第 85-87 行计算 B 的分量。第 88 行给出了它的大小 magB=|B|。第 89 行确定 p 是否在阴影区内。如果是，则执行第 90 行。这是 V 与^ˇl 的点积，它决定了 p 是否位于行星朝向光源的一侧，在这种情况下，它与行星的黑暗面相对，不在阴影区。这是必要的，因为第 78-89 行的阴影算法没有进行这种区分。如果 p 确实位于阴影区内的暗侧，则颜色在第 91 行被设置为中灰色。

你会注意到在上面的图片中，光环内有一条黑色的带。这是因为土星环在那一带有一个空洞:那里没有粒子来反射光线；你看到的是背景色，“午夜蓝”，透过。这就产生了一个问题，因为阴影颜色会在空白区域覆盖背景颜色。第 93 行和第 94 行将其重建为“午夜蓝”。

现在，波段颜色已经确定，您可以绘制环。这是通过绘制短线段来完成的。第 97-100 行计算第一段的起始位置。参考图 10-6 ，线 100-101 确定该线段是否在行星前面，在这种情况下，绘制该线段。第 103-108 行确定它是否在行星后面，在这种情况下，它不被标绘。这是通过计算该点的全球坐标与行星中心的距离 c 来实现的。第 107 行说，如果 c 大于球的半径乘以 1.075，那么画出线段。包括因子 1.075 是为了防止线段蚕食球体的边缘。需要经过这个逻辑；否则，位于球体半径内的前方可见线段将不会被绘制。

关于清单 10-1 生成的上述图像，可以注意到两点。首先是颜色。美国宇航局的照片显示了灰色的色调，几乎没有颜色。但是许多土星的观察者描述它有金色的色调，因此我选择了它的颜色。任何摄影师都知道，在摄影图像中捕捉物体的真实颜色是很困难的；这很大程度上取决于入射光的颜色和图像捕捉介质。也许最好依靠观星者的观察。如果你不同意清单 10-1 生成的图像中的颜色，你可以通过改变程序中的 clr 定义来修改它们。第二件要注意的事情是在图 10-5 中跟随行星曲率的阴影曲率。它表明着色算法按预期工作。

关于程序的使用，你可以在第 24-26 行改变入射光的方向，在第 32-34 行改变旋转角度。清单 10-1 需要一些时间来运行，所以请耐心等待。

1   """
2   SATURN
3   """
4  
5   import numpy as np
6   import matplotlib.pyplot as plt
7   from math import sin, cos, radians, sqrt
8  
9   plt.axis([0,150,100,0])
10  plt.axis('off')
11  plt.grid(False)
12 
13  print('running')
14  #—————————————————parameters
15  g=[0]*3
16
17  xc=80 #———sphere center
18  yc=50
19  zc=0
20
21  rs=25 #———sphere radius
22
23  lx=-1 #———light ray unit vector components
24  ly=0
25  lz=0
26
27  IA=0
28  IB=.8
29  +n=2
30
31  Rx=radians(-20)
32  Ry=radians(0)
33  Rz=radians(30)
34
35  #————————same as SHADESPHERE—————–
36
37  #———————————————————rings
38  alpha1=radians(-10)
39  alpha2=radians(370)
40  dalpha=radians(.5)

41
42  r1=rs*1.5
43  r2=rs*2.2
44  dr=rs*.02
45  deltar=(r2-r1)/7 #———ring band width
46
47  #—————————————rotate ring point p which is at r, alpha
48  for r in np.arange(r1,r2,dr):
49      for alpha in np.arange(alpha1,alpha2,dalpha):
50          xp=r*cos(alpha)
51          yp=0
52          zp=-r*sin(alpha)
53          rotx(xc,yc,zc,xp,yp,zp,Rx)
54          xp=g[0]-xc
55          yp=g[1]-yc
56          zp=g[2]-zc
57          roty(xc,yc,zc,xp,yp,zp,Ry)
58          xp=g[0]-xc
59          yp=g[1]-yc
60          zp=g[2]-zc
61          rotz(xc,yc,zc,xp,yp,zp,Rz)
62          xpg=g[0]
63          ypg=g[1]
64
65  #—————————————————select ring band color
66      if r1 <= r < r1+1*deltar:
67          clr=(.63,.54,.18)
68      if r1+1*deltar <= r <= r1+2*deltar:
69          clr=(.78,.7,.1)
70      if r1+2*deltar <= r <= r1+3*deltar:
71          clr=(.95,.85,.1)
72      if r1+3*deltar <= r <= r1+4*deltar:
73          clr=(.87,.8,.1)
74      if r1+5*deltar <= r <= r1+7*deltar:
75          clr=(.7,.6,.2)
76
77  #———————————————————————shadow
78      magu=sqrt(lx*lx+ly*ly+lz*lz)
79      ux=-lx/magu
80      uy=-ly/magu
81      uz=-lz/magu
82      vx=xc-xpg
83      vy=yc-ypg
84      vz=zc-zpg
85      Bx=uy*vz-uz*vy
86      By=uz*vx-ux*vz
87      Bz=ux*vy-uy*vx
88      magB=sqrt(Bx*Bx+By*By+Bz*Bz)
89      if magB < rs: #—————————if in the shadow region
90          if vx*lx+vy*ly+vz*lz <= 0: #———if v points toward light source
91              clr=(.5,.5,.2) #———shadow color
92
93      if r1+4*deltar <= r <= r1+5*deltar: #———overplot empty band
94          clr='midnightblue' #———with background color
95
96  #——————————————————–plot line segment
97      if alpha == alpha1:
98          xstart=xpg
99          ystart=ypg
100     if zpg <= zc: #–front (z axis points into the screen)
101         plt.plot([xstart,xpg],[ystart,ypg],linewidth=2,color=clr)
102
103     if zpg >= zc: #–back
104         a=xpg-xc
105         b=ypg-yc
106         c=sqrt(a*a+b*b)
107         if c > rs*1.075: #——plot only the visible portion of rings
108             plt.plot([xstart,xpg],[ystart,ypg],linewidth=2,color=clr)
109         xstart=xpg
110         ystart=ypg
111
112 plt.show()

Listing 10-1
Program SATURN

10.2 太阳辐射

本节展示了一个典型的场景:使用 Python 来绘制和标记代表数学函数的曲线。这里你画出了马克斯·普朗克的辐射光谱。你用它来表示太阳发出的能谱，计算太阳的总功率输出，以及到达地球的量，这个量叫做太阳常数。这一部分的科学方面非常有趣，就像它们的发展历史一样。Python 编程的主要好处是可以看到程序如何执行数值积分、建立图表以及显示数值数据。

光子和太阳

像所有辐射体一样，太阳以光子的形式发射电磁能。我们知道光子以不同的频率或波长发射，波长与频率成反比，如

)

(10-5)其中λ是波长，s _m 是光在介质中的速度，ν是波在该介质中传播的频率。由于我们主要关心的是光在真空中的传播(即从太阳到地球)，s _m =c 其中 c 是光在真空中的速度，因此等式 10-5 变成了

)

(10-6)

像太阳能这样的函数，当用一系列频率或波长来表示时，叫做光谱。在电磁辐射的情况下，我们最关心的是不同频率的光的功率，或者说是波长。这被称为功率谱。例如图 10-8 所示的曲线，其中绘制了功率谱密度(通常简称为功率谱，S(λ))与波长λ的关系。这条曲线源于方程 10-7 。

图 10-8

Max Planck’s Solar Spectrum

太阳发出的大部分频率，从高频、短波长的紫外线到低频、长波长的红外线，我们的人眼是看不见的。我们只能看到光谱的一小部分，幸运的是，它位于太阳发射的功率谱的峰值附近。这一定让我们狩猎采集的祖先很高兴，因为这使他们能够在一天的早些时候和晚些时候狩猎采集。虽然我们可以为此感谢太阳功率谱的形状，但这也是我们眼睛生物学的一个特征，如果我们相信查尔斯·达尔文的话，这可能是在太阳最大功率输出附近的频率下优化进化的。

10.2.2 普朗克黑体辐射

如前所述，光是光子。但是什么是光子呢？我们知道光子是电磁能量的量子化形式。但是在 19 世纪晚期，这仍然是一个谜。人们曾多次试图解释受热物质发出的光谱。例如，当我们加热一根铁拨火棒时，起初我们看不到颜色有任何变化，但在达到一定温度后，我们可以看到它发出一系列颜色:暗红色、亮红色、橙色、黄色、白色、蓝色，然后是紫色。这些颜色对应于物体发出的电磁辐射的不同频率。早期试图解释这一现象是基于当时称为麦克斯韦方程组的经典理论。这些方程描述了一个电磁场，其中电磁能量被假设为一个平滑的连续体。尽管做了许多尝试，这种方法还是无法解释所观察到的现象。

当时许多科学家都在努力解决这个问题。然后在 1900 年，德国物理学家马克斯·普朗克给一位同事寄了一张明信片。在背面，他写了一个精确描述光谱的方程式。与当时流行的理论相反，普朗克的突破是假设电磁场不是能量的连续体。相反，他猜测电磁能以离散的包存在，而不是像麦克斯韦方程所假设的那样是一个连续的场。这导致了他的突破性公式，如方程式 10-7 所示。1900 年 12 月 14 日，他向德国物理学会提交了他的想法，这一天被称为量子力学的诞生。他的方程被称为黑体辐射公式。稍后你会看到更多。

1901 年，普朗克在《物理学年鉴》的一篇文章中发表了他的结果，在这篇文章中，他假设电磁能量只能由一个源，如太阳，以离散的能量包的形式释放，而不是连续的波。因为人们知道光表现出波动特性，1905 年阿尔伯特·爱因斯坦扩展了这一观点，提出普朗克的离散“包”只能以离散“波包”的形式存在。他称这种包为 Lichtquant 或“光量子”后来，在 1928 年，阿瑟·康普顿使用了“光子”一词，这个词来源于希腊语 Phos，意为光。

普朗克还假设波包的来源是热激发的电荷，每个电荷以特定的频率发射一包电磁能量，即光子。以相同频率发射光子的电荷越多，以该频率发射的光的功率就越大。他进一步提出理论，认为波包的能量只能出现在特定的固定能级或状态。

回到 1900 年，普朗克写在明信片背面的方程，它预言了一个黑体发出的光 S(λ)的功率谱是

)

(10-7)其中 c 是光速(m/s)，h 是普朗克常数(J s)，λ是波长(m)，K 是玻尔兹曼常数(J/K)，T 是温度(K)，ε是辐射体表面的发射率。ϵ本质上是对表面辐射能力有效性的一种度量。它的范围可以从 0 到 1。正如你可能想象的那样，这个等式的发展背后有很多想法，我在这里就不赘述了。

在过去的 100 多年里，这种关系经受住了时间的考验，并给出了非常准确的结果。如图 10-8 所示，通常被称为普朗克黑体辐射公式。它同样适用于所有辐射体和太阳。尽管太阳看起来肯定不像一个黑体，但就其辐射特性而言，它的行为就像一个黑体——一个非常热的黑体。作为一个类比，你可能会认为太阳是一个非常热(大约 5800 K)的黑色火炉，发出非常明亮的光。

图 10-8 显示了由方程 10-7 预测的太阳输出光谱(红色曲线)。这被称为功率谱密度，或简称为功率谱。曲线上的每个点给出了相应波长λ下的功率密度 S(λ)。稍后将解释图中所示的绿色带。

10.2.3 太阳总功率输出

图 10-8 中显示的量 S(λ)是功率密度。什么是功率密度，它与简单功率有何不同？注意方程 10-7 中 S(λ)的单位是每立方体积的功率。这些是密度的单位。你可能会认为这个“密度”类似于质量密度，单位是每立方体积的质量。在 S(λ)的情况下，你处理的是功率密度。

等式 10-7 中使功率谱类似于太阳输出而不是任何其他黑体输出的特征是温度 T。对于太阳，T 约为 5800 K。要获得太阳在λ ₁ 到λ ₂ 的带宽上发射的功率 P(λ)，需要对该频带上的 S(λ)求和。在微积分中，这相当于求 S 对λ的积分，相当于求 S(λ)曲线下这些极限之间的面积。

)

(10-8)

用方程式 10-7 这就变成了

)

(10-9)

等式 10-9 给出了太阳在λ ₁ 到λ ₂ 的带宽上发射的功率。它等于 S(λ)乘以无穷小带宽 dλ的积分。换句话说，如果你沿着 S(λ)曲线选择一个点，如图 10-9 所示，将它乘以 dλ，然后将从λ ₁ 到λ ₂ 的所有值相加，你将得到在λ ₁ 到λ ₂ 波段的波长所发射的总电磁功率。这是从λ ₁ 到λ ₂ 的 S(λ)曲线下的区域。

为了得到整个太阳光谱产生的能量，你从波长λ=0 延伸到λ=∞，对方程 10-9 进行积分。对于那些喜欢对方程 10-9 进行数学积分的人，我在附录 b 中展示了如何进行积分。积分就是找到 S(λ)曲线下的面积。你可以用数字来避免计算。为此，将无限小的波段 dλ替换为有限宽度∏λ的小波段，并将积分替换为求和，如在

)

(10-10)中，其中 I 是指以λ _i 为中心的第 i ^个波段，N 是λ _i 和λ _N 之间宽度∏λ的波段数。图 10-9 显示了宽度为∏λ的典型带。为了说明的目的，所示的带的宽度被夸大了。实际上，它应该显得更窄。

图 10-9

Numerical integration of power S(λ)dλ emitted by spectrum band_i across a .01 μm bandwidth at λ_i=1.5 (produced by Listing 10-2)

方程 10-10 是方程 10-9 的近似，因为它假设 S(λ)的值在每个波段∈λ的宽度上是常数。然而，如果选择的∏λ足够小，曲线 S(λ_I-∏λ/2)→S(λ_I+∏λ/2)可以用带宽∏λ上的常数值 S(λ _i )来近似，在这种情况下，结果可能相当精确。利用这个简单的积分方案，频带中的功率等于频带的矩形面积。虽然您可以使用更复杂的集成方案，但是这个方案很简单，易于编程，并且足以满足您的需求。

让我们计算小波段∏λ上的波长发射的功率 P(λ)。图 10-9 是图 10-8 中以λ= 1.5 μm 为中心的谱带的放大图。这可能被认为是等式 10-10 中的典型带 _i 。清单 10-2 评估了该带宽内波长产生的功率。曲线 S(λ)已根据方程式 10-7 生成。根据这个简化的积分方案，这个频带产生的功率，也就是它的矩形面积，由

)

(10-11)给出

在绘制图 10-9 的清单 10-2 中，条带的面积是根据方程 10-11 计算的。∏λ的大小是任意的。在程序中，它是参数 dla，设置为 0.01 X10⁶米或 0.01 μm。无论∏λ是大还是小，它发出的功率都是该带宽内波长辐射的功率。带宽越宽，产生的功率越大，带宽越窄，产生的功率越小。稍后，当您对整个 S(λ)曲线下的面积进行数值积分以获得太阳在其整个光谱上辐射的总功率时，选择较小的∏λ值将会导致更精确的结果。

在图 10-9 中，波段显示为λ=1.5 um。根据程序计算，对应的 S 值为 1.164x10 ⁷ MW/m ³ 。带宽为 0.01 微米，相当于 1.0x 10⁸米，这个频段产生的功率为(1.164 X10⁷)x(1 X10⁸)= . 1164 MW/m²，大约相当于一个小型发电厂的发电量。

注意，S(λ)的单位将与输入参数的单位一致:光速、普朗克常数、玻尔兹曼常数和波长λ。这些参数的单位应该彼此一致。为了避免混淆，在这项工作中，当计算方程 10-7 时，你将保持所有这些量在米的空间维度中。S(λ)将有单位(J/s)/m ³ ，与 W/m ³ 相同。如果需要另一种功率尺寸的输出，如 kW 或 MW，可以在计算 S(λ)后进行转换，即将 S(λ)乘以 10³得到千瓦，或将 S(λ)乘以 10⁶得到兆瓦。计算某个波段上发射的功率时，该波段的宽度∈λ也应该以米为单位。例如，1.5μm 应指定为 1.5x 10^-6m。将λ米乘以 10 ⁺⁶ 可以将米转换回微米微米，以便以后显示或用于其他目的。这在列表 10-2 中显示，在直线 lag=la*10**6 中绘制 s 曲线。

图中显示 S _i 的值为 1.164x10 ⁷ MW/m ³ 。S(λ)轴表示 11.64 MW/m³X10⁶的值，这表示 11.64 的值已经乘以 10⁶用于显示。这将使其实际值为 11.64x10 ⁺⁶ ，等于程序计算的值。这在图上显示为 1.164 X10⁺⁷MW/m³。

在创建了图 10-9 的列表 10-2 中，截面图 S 曲线求解方程 10-7 得到波长 la 的值，该值以增量 dla 从 la=lamin 到 lamax。代码中的注释描述了 s 单位的演变。如公式 10-7 所示，当参数如“建立参数”一节所示时，s 的单位以焦耳/秒/立方米开始(记住 s(λ)是密度)。因为 1 焦耳每秒定义了瓦特，所以单位是瓦特每立方米。这些被转换成兆瓦每立方米，然后换算成以单位(MW/m³)X10⁶表示的垂直轴，作为变量 sg。10⁶因子表示实际值已经乘以该数值。接下来，绘制绿色带，并显示温度和发射率的值。

使用公式 10-7 计算λ=1.5 时的 S(λ)值，转换为 MW/m ³ ，然后乘以带宽 dl = . 01x 10⁶，得到 pl MW/m ² ，即该带宽内的功率。程序的其余部分显示数据并清理地块。

"""
BANDINTEGRAL
"""

import numpy as np
import matplotlib.pyplot as plt

#---------------------------------------------------------- set up axes
ymax=20
plt.axis([1.,2.,0,ymax])
plt.xlabel('Wavelength $\lambda$ ($\mu$m)')
plt.ylabel('S($\lambda$) (MW/m$^{3}$) x 10$^{-6}$')
plt.grid(True)
plt.title('Max Planck's Solar Spectrum - Band Integral')

#------------------------------------------------ establish parameters
c=2.9979*(10.**8)           # speed of light in a vacuum m/s
h=6.63*(10.**-34)           # Planck's Constant J.s
kb=1.38*(10**-23)           # Boltzmann's Constant J/K

t=5800\.                     # temperature K
e=1.0                       # emissivity

lamin=.01*10**-6            # starting wavelength m
lamax=2.*10**-6             # ending wavelength m
dla=.01*10**-6              # incremental wavelength m

#———————————————————————————————————————————— plot s curve
for la in np.arange(lamin,lamax,dla):
    a1=2.*np.pi*c*c*h/(la**5.)
    a2=h*c/(la*kb*t)
    sl=e*a1/(np.exp(a2)-1.)                 # J/s/m³ = W/m³
    sl=sl*10**-6                            # MW/m³
    slg=sl*10**-6                           # scale plot at 10^-6 scale
    lag=la*10**6                            # scale to plot at 10⁶ scale
    plt.scatter(lag,slg,s=1,color='r')
#——————————————————————————————————————————————— plot band
plt.plot([1.495,1.495],[0.,11.64],color='g')
plt.plot([1.4975,1.4975],[0.,11.64],color='g')
plt.plot([1.5,1.5],[0,11.64],color='g')
plt.plot([1.5025,1.5025],[0.,11.64],color='g')
plt.plot([1.5005,1.505],[0.,11.64],color='g')

#——————————————————————————————— plot temperature and emissivity
d=str(t)
plt.text(1.6,15,'T=')
plt.text(1.65,15,d)
plt.text(1.6,14,'e=')
d=str(e)
plt.text(1.65,14,d)

#———————————————————— calculate s and band power pl at lambda=1.5
la=1.5*10**-6
a1=2.*np.pi*c*c*h/(la**5.)
a2=h*c/(la*kb*t)
sl=e*a1/(np.exp(a2)-1.)  # J/s/m³ = W/m³
sl=sl*10**-6             # MW/m³
dl=.01*10**-6            # bandwidth m
pl=sl*dl

#———————————————————————————————— plot results and labels
plt.plot([1.53,1.59],[11.6,11.6],'k')
plt.text(1.6,11.5,'si=')
d='%7.3e'%(sl)
plt.text(1.65,11.5,d)
plt.text(1.83,11.5,'MW/m³')

plt.arrow(1.4,5,.085,0,head_width=.5,head_length=.01,linewidth=.2)
plt.arrow(1.6,5,-.085,0,head_width=.5,head_length=.01,linewidth=.2)

plt.text(1.15,5,'$\Delta \lambda$=')

dle='%7.3e'% (dl)
dls=str(dle)
plt.text(1.18,5,dls)
plt.text(1.35,5,'m')')

plt.text(1.145,4,'=')
dl=dl*10**6
dle='%7.3e'%(dl)
dls=str(dle)
plt.text(1.18,4,dls)
plt.text(1.35,4,'um')

plt.text(1.35,16.5,'s($\lambda$)')
plt.text(1.52,2.5,'power$_{i}$=')
pl='%7.3e'%(pl)
pl=str(pl)
plt.text(1.65,2.5,pl)
plt.text(1.823,2.5,'MW/m²')
plt.text(1.45,-1.1,'$\lambda_{i}$=1.5')

plt.show()

Listing 10-2Program BANDINTEGRAL

接下来我们来看一下如图 10-8 所示的马克斯普朗克的整个黑体谱。因为使用的温度是太阳的温度，大约 5800^{k，所以它被称为“马普太阳光谱”}

生成此图的程序(列表 10-3 )遵循前面程序(列表 10-2 )中的逻辑，但这里您对从λ= . 01 X10⁶到 10 . X10⁶米(. 01μm 到 10.μm)的各个频带功率求和，以获得整个(几乎)S(λ)曲线下的面积。你在清单 10-2 中看到的波段显示为λ=1.5 um。

这里使用的过程是简单地沿着波长前进，计算每个波长的 S(λ)值，乘以∏λ以获得该波段内的功率，然后根据等式 10-10 对每个波段产生的功率求和。这将给出所有波长发射的总功率。

为了更精确地测量 S(λ)曲线下的总功率，您将积分范围扩展至 10x 10⁶米。这将是每平方米太阳表面发出的总光谱能量。然后你把它乘以太阳的球面面积，就得到太阳发出的总功率，这就是所谓的太阳光度。在图 10-8 中，它被称为“太阳能总输出”如图所示，程序计算出的值为 3.816x10 ²⁶ 瓦。这与公布的数值非常一致。

许多研究人员使用 e=1.0 作为发射率，这是一种理想化的假设，即太阳是一个完美的辐射体(你可以假设它不是)。这里你使用 e = .984 的发射率。当你使用普朗克光谱计算太阳常数时，太阳常数已经由卫星测量(下一节)，你必须在你的计算中降低太阳的温度，或者将其发射率降低到小于 1.0，以便使结果与测量值一致。如果你选择保持 5800 K 的太阳温度，那么你必须将发射率降低到 0.984 才能获得一致。另一个选择，你会看到，是保持 e=1.0，把太阳的温度降到 5777 K。

"""
PLANCKSSOLARSPECTRUM
"""

import numpy as np
import matplotlib.pyplot as plt
#—————————————————————————————————————————————————— set up axes
ymax=100
plt.axis([0,3,0,ymax])
plt.xlabel('Wavelength _ (_m)')
plt.ylabel('S(λ) (MW/m^{3}) x 10^-6')

plt.grid(True)
plt.title('Max Planck's Solar Spectrum')
#———————————————————————————————————————— establish parameters
c=2.9979*(10.**8)           # speed of light in a vacuum m/s
h=6.63*(10.**-34)           # Planck's Constant J.s
kb=1.38*(10**-23)           # Boltzmann's Constant J/K
e=.984                      # emissivity
t=5800\.                     # K
lamin=.01*10**-6            # m
lamax=10.*10**-6            # m
dla=.01*10**-6              # m
st=0\.                       # set area under s curve to zero
#—————————————————————————————— plot s curve and calculate area
for la in np.arange(lamin,lamax,dla):
    a1=2.*np.pi*c*c*h/(la**5.)
    a2=h*c/(la*kb*t)
    sl=e*a1/(np.exp(a2)-1.)                # W/m³
    sl=sl*10**-6                           # MW/m³
    bandarea=sl*dls                        # band area MW/m²
    st=st+bandarea                         # sum band areas MW/m²
    slg=sl*10**-6                          # scale to plot
    lag=la*10**6                           # scale to plot
    plt.scatter(lag,slg,s=1,color='r')
#—————————————————————————————————— multiply the Sun's surface area
ds=1.39*10**9                # Sun's diameter m
spas=np.pi*ds**2\.            # Sun's spherical area m²
to=spas*st                   # Sun's total output MW
to=to*10**6                  # Sun's total output W
#———————————————————————————————————————————————— plot results
plt.text(.8,58.,'5800')
plt.text(1.05,58, '°K')
plt.plot([.39,.39],[-0.,100.],'b–')
plt.plot([.7,.7],[-0.,100.],'b–')

plt.text(.3,-10,'.390')
plt.text(.6,-10,'.700')
plt.text(.15,90.,'UV')
plt.text(.8,90.,'long wave infrared')
plt.arrow(1.75,91.,.8,0.,head_width=1.,head_length=.1,color='r')
plt.text(1.2,40.,'total solar output =')
so='dd=str(so)
plt.text(2.1,40,dd)
plt.text(2.7,40,'W')
plt.text(1.2,34,'emissivity =')
e=str(e)
plt.text(1.8,34,e)
plt.text(.5,75.,'v')
plt.text(.53,70.,'i')
plt.text(.5,65.,'s')
plt.text(.53,60.,'i')
plt.text(.5,55.,'b')
plt.text(.53,50.,'l')
plt.text(.5,45.,'e')
plt.plot([1.49,1.49],[0.,11.61],color='g')
plt.plot([1.5,1.5],[0.,11.61],color='g')
plt.plot([1.51,1.51],[0.,11.61],color='g')
#—————————————————— calculate s at la=1.5x10^-6 m and band power pband
laband=1.5*10**-6
a1=2.*np.pi*c*c*h/(laband**5.)
a2=h*c/(laband*kb*t)
sband=a1/(np.exp(a2)-1.)
sband=sband*10**-12
pband=sband*dla # MW/sq meter
pband=pband*10**6 # W/sq meter
#———————————————————————————————————————————————— plot band
plt.plot([1.55,1.7],[12.5,15.],color='k')
plt.text(1.72,14.,' p=')
pband='pband=str(pband)
plt.text(1.9,14,pband)
plt.text(2.4,14,'MW/m²')

plt.arrow(1.35,5,.1,0,head_width=1, head_length=.05, ec="k", fc="k")
plt.arrow(1.65,5,-.1,0,head_width=1, head_length=.05, ec="k", fc="k")
plt.text(.82,4.9,'Δλ = :01μm' )

plt.show()

Listing 10-3Program PLANCKSSOLARSPECTRUM

10.3 地球辐照度

图 10-10 显示了到达地球的太阳辐射光谱。图 10-11 显示地球围绕太阳旋转。这是用来计算地球截获的太阳总发电量的模型，即太阳常数。两个球体之间的距离平均为 1 天文单位，约为 93，000，000 英里。它在一个轨道上是变化的。标有 _p 的圆盘的面积等于地球的横截面积。A _p 截取的太阳能负责加热地球。

图 10-10

Svpectrum of solar power reaching Earth, the Earth’s solar irradiance, produced by Listing 10-3, which has been modified by inclusion of the inverse square law

注意这些值比图 10-8 中的值低多少。这是因为到达地球并被 A _p 拦截的太阳输出的功率强度根据

)

(10-12)的平方反比定律随着从太阳到地球的距离而减小，其中 p _s 是太阳表面的功率强度，p _p 是被 A _p 拦截的强度，r _s 是太阳的半径，r _{A p 截获的总功率因此为}

(10-13)

(10-14)

当等式 10-13 包含在列表 10-3 中时，频谱减少到图 10-10 。再次注意，由于平方反比定律，这些值比图 10-8 中的值低了多少。

P _p ，即到达地球大气层顶部的太阳能功率，称为太阳常数。卫星测得的数值约为 1361 瓦/米 ² 。其中大约 30%是通过反照率效应从地球表面和大气反射的，如雪、冰、云、水等。其余的被地球吸收了。其中大部分被重新辐射回太空，让地球达到热平衡。地球也是一个热(暖)体，它向太空展示自己的热辐射。但是一些应该再辐射的物质被包括二氧化碳在内的温室气体阻挡，导致全球变暖。众所周知，所有这些都是气候研究人员正在积极研究的。

10.3.1 地球太阳模型

图 10-11 显示地球绕太阳运行，清单 10-4 包含代码。

图 10-11

The Earth-Sun Model produced by Listing 10-4

"""
EARTHSUN
"""

import matplotlib.pyplot as plt
import numpy as np
from math import radians, sin, cos, sqrt

plt.axis([-100,150,-100,150])

plt.grid(False)
plt.axis('off')
sfx=2.5/3.8

#———————————————————————background
for x in range(-100,150,2):
    for y in range(-100,150,2):
        plt.scatter(x,y,s=40,color='midnightblue')

phimin=0.
phimax=2.*np.pi
dphi=phimax/100.

rs=40.
re=20.

ys=15.
ye=2.

xos=50.
yos=0.
zos=0.

#———————————————————Sun's core
plt.scatter(xos,yos,s=4300,color='yellow')

#———————————————————Sun horizontals
rx=radians(20)

for ys in np.arange(-rs,rs,5):
    for phi in np.arange(phimin,phimax,dphi):
        rp=np.sqrt(rs*rs-ys*ys)
        xp=rp*np.sin(phi)
        yp=ys
        zp=rp*np.cos(phi)
        px=xos +sfx*xp*1\. +yp*0\. +zp*0.
        py=yos +xp*0\. +yp*np.cos(rx) -zp*np.sin(rx)
        pz=zos +xp*0\. +yp*np.sin(rx) +zp*np.cos(rx)
        if pz > 0 :
            plt.scatter(px,py,s=1,color='red')

#————————————————————Sun verticals
alphamin=0.
alphamax=2.*np.pi
dalpha=alphamax/30.

for alpha in np.arange(alphamin,alphamax,dalpha):
    for phi in np.arange(phimin,phimax,dphi):
        xp=rs*np.sin(phi)*np.sin(alpha)
        yp=rs*np.cos(phi)
        zp=rs*np.sin(phi)*np.cos(alpha)
        px=xos +sfx*(xp*1\. +yp*0\. +zp*0.)
        py=yos +xp*0\. +yp*np.cos(rx) -zp*np.sin(rx)
        pz=zos +xp*0\. +yp*(np.sin(rx)) +zp*np.cos(rx)
        if pz > 0 :
            plt.scatter(px,py,s=1,color='red')

#——————————————————————Earth's clouds
xoe=-50.
yoe=20.
zoe=-10.

plt.scatter(xoe,yoe,s=800,color='white')

#———————————————————— Earth horizontals
rx=20.*np.pi/180.
dphi=phimax/100.

for ys in np.arange(-re,re,2):
    for phi in np.arange(phimin,phimax,dphi):
        rp=np.sqrt(re*re-ys*ys)
        xp=rp*np.sin(phi)
        yp=ys
        zp=rp*np.cos(phi)
        px=xoe +sfx*(+xp*1\. +yp*0\. +zp*0.)
        py=yoe +xp*0\. +yp*np.cos(rx) -zp*np.sin(rx)
        pz=zoe +xp*0\. +yp*(np.sin(rx)) +zp*np.cos(rx)
        if pz > 0 :
            plt.scatter(px,py,s=.1,color='#add8e6')

#—————————————————————Earth verticals
alphamin=0.
alphamax=2.*np.pi
dalpha=alphamax/30.

for alpha in np.arange(alphamin,alphamax,dalpha):
    for phi in np.arange(phimin,phimax,dphi):
        xp=re*np.sin(phi)*np.sin(alpha)
        yp=re*np.cos(phi)
        zp=re*np.sin(phi)*np.cos(alpha)
        px=xoe +sfx*(xp*1\. +yp*0\. +zp*0.)
        py=yoe +xp*0\. +yp*np.cos(rx) -zp*np.sin(rx)
        pz=zoe +xp*0\. +yp*(np.sin(rx)) +zp*np.cos(rx)
        if pz > 0 :
            plt.scatter(px,py,s=.1,color='#add8e6')

plt.arrow(xos-rs*sfx-3,yos+2,xoe-(xos-rs*sfx)+re+3,yoe-yos-6.2,color='r',
                              head_length=4.,head_width=3.)
plt.text(-14,16,'1 AU',color='white')
plt.text(80,-29,'Sun',color='white')
plt.text(-84,10,'Earth',color='white')

#———————————————————————front orbit
deltamin=0.*np.pi/180.
deltamax=195.*np.pi/180.
ddelta=deltamax/60.

for delta in np.arange(deltamin,deltamax,ddelta):
    r=108./sfx
    xp=r*np.cos(delta)
    yp=0.
    zp=r*np.sin(delta)
    px=xos +sfx*(xp*1\. +yp*0\. +zp*0.)
    py=yos +xp*0\. +yp*np.cos(rx) -zp*np.sin(rx)
    pz=zos +xp*0\. +yp*(np.sin(rx)) +zp*np.cos(rx)
    plt.scatter(px,py,s=1,color='white')

#———————————————————————back orbit
deltamin=220.*np.pi/180.
deltamax=360.*np.pi/180.

for delta in np.arange(deltamin,deltamax,ddelta):
    r=108./sfx
    xp=r*np.cos(delta)
    yp=0.
    zp=r*np.sin(delta)
    px=xos +sfx*xp*1\. +yp*0\. +zp*0.
    py=yos +xp*0\. +yp*np.cos(rx) -zp*np.sin(rx)
    pz=zos +xp*0\. +yp*(np.sin(rx)) +zp*np.cos(rx)
    plt.scatter(px,py,s=1,color='white')

#———————————————————Ap disc
xoc=xoe+re*sfx
yoc=yoe-2.5
zoc=zoe
rc=.83*re
phi1=0
phi2=2*np.pi
dphi=(phi2-phi1)/200
ry=-25*np.pi/180

for phi in np.arange(phi1,phi2,dphi):
    xc=xoc
    yc=rc*np.sin(phi)
    zc=rc*np.cos(phi)
    px=xoc+zc*np.sin(ry)
    py=yoc+yc
    pz=zoc+zc*np.cos(ry)
    plt.scatter(px,py,s=.03 ,color='white')

plt.scatter(xoe+re*sfx,yoe-2,s=6,color='white')

plt.arrow(-20,60,(xoe+re*sfx)+24,(yoe+re/2)-60-2,color='white',
                             linewidth=.5,head_width=2.,head_length=3)
plt.text(-18,60,'A
p
',color='white')

plt.show()

Listing 10-4Program EARTHSUN

10.4 总结

在本章中，你已经看到了 Python 图形编程的一些典型应用。在第一部分中，你看到了创建土星的图像是多么的容易。这颗行星的身体利用了前几章开发的球体和阴影算法；这些带是通过在 x，z 平面上构建同心环而形成的。阴影算法需要一点原始几何。然后整个东西在空间绕 x，y，z 轴旋转。在关于太阳辐射的章节中，您学习了太阳辐射的物理学，尤其是马克斯·普朗克的黑体辐射公式，以及 Python 构建技术插图的能力。特别值得注意的是用于绘制的缩放变量的技术。然后你学习了如何构建图像，如图 10-8 所示，显示了用于理解地球太阳辐照度的模型。

十一、从哪里获得 Python

互联网上有几个地方可以下载各种版本的 Python。我将 Anaconda 与 Spyder ² 和 Python 3.5 一起使用。这可以从 https://docs.continuum.io/anaconda/install 的 Continuum Analytics 下载。

很洒脱；请按照说明操作。虽然我使用 Python 3.5，但我建议使用最新版本。

桌面上会出现一个图标。如果没有，请查看您的已安装程序列表，并将其拖到桌面上。双击它以运行环境。您将在左侧窗格中输入 Python 脚本。输入程序代码后，单击顶部的运行按钮或按键盘上的 F5 键。您可能会被告知打开一个新的控制台。单击顶部的控制台按钮，然后选择“打开 IPython 控制台”选项。试着再运行一次。结果应该出现在右下方的窗格中。

右上方有一个显示变量状态的窗格。我从来不用它；事实上，我关闭它是为了给输出留出更多空间。如果我想知道一个特定的变量在做什么，我通常会在程序中放一个 print 语句。变量的历史记录将出现在输出窗格中。

如果您发现您的程序正在做意想不到的事情，打开一个新的控制台并重新运行程序有时会有所帮助。

十二、普朗克辐射定律和斯特凡-玻尔兹曼方程

在第十章中，向你介绍了马普著名的黑体辐射方程:

)

(B-1)

表面在带宽λ ₁ → λ ₂ 上发射的功率为

)

(B-2)

用方程 B-1，这就变成了

)

(B-3)

在第十章中，你对方程 B-3 进行了数值积分。在这里，您将对其进行数学积分，并表明它可用于推导黑体辐射的 Stefan-Boltzmann 定律

)

(B-4)其中 t 是表面的绝对温度，p 是单位面积辐射的功率，k _B 是 Boltzmann 常数，h 是 Planck 常数，c 是光速，ϵ是表面的发射率。从区域 A 的表面辐射的功率是

)

(B-5)，其中

)

(B-6) σ被称为 Stefan-Boltzmann 常数。等式 B-4 将表面辐射的功率强度与其温度 t 的四次方相关联。该等式通常用于科学和工程中。

为了进行导致方程 B-4 的积分，你从普朗克的辐射方程开始(也在上面的 12-1 中示出):

)

(B-7)你想要从λ=0 到λ=∞对其积分以获得每单位的总功率

所有波长辐射的区域 p。让 C ₁ =ϵ2πhc ² 和

)

，你得到

)

(B-8)

如果你做了下面的替换

)

(B-9)稍微忙乱一下你就有了

)

(B-10)

利用众所周知的:)关系

)

(B-11)，将 C ₁ 和 C ₂ 代入方程 B-10，得到

)

(B-12)，与上面的方程 B-4 相同。