土星的行星環(huán)非常出名。雖然木星、土星、天王星和海王星也有環(huán)，但土星環(huán)是我們太陽系中最大、最亮、最廣為人知的行星環(huán)。

它由小到灰塵的顆粒，大到巨石的物體組成。這些物體的成分主要是冰，一般認(rèn)為是彗星或較大的小行星與土星的一顆衛(wèi)星相撞時(shí)產(chǎn)生的，兩者都撞成了小碎塊。在遠(yuǎn)古時(shí)代，土星就已為人所知，但直到1610年，伽利略才首次用望遠(yuǎn)鏡對(duì)它進(jìn)行觀測(cè)。

這個(gè)行星以羅馬農(nóng)業(yè)之神土星Saturn命名為，也就是我們每個(gè)星期的第六天Saturday。

圖1至圖5中的圖像是由本文文末的代碼生成的。每張圖都呈現(xiàn)不同的方向角，圖標(biāo)題中有相應(yīng)的說明。

圖標(biāo)題中還列出了入射光線的單位矢量分量，例如lx=+0.707，ly=+0.707，lz=0 表示左上象限中的光源;lx=-1，ly=0，lz=0表示來自右側(cè)的光源。在圖像中請(qǐng)注意行星在環(huán)上投射的陰影，尤其是在圖5中能夠看到行星輪廓的曲率。

▲圖1 包含土星環(huán)和陰影的土星1：Rx=-20°, Ry=0, Rz=-10°, lx=1, ly=0, lz=0

▲圖2 包含土星環(huán)和陰影的土星2：Rx=-8°, Ry=0, Rz=-30°, lx=0.707, ly=.707, lz=0

▲圖3 包含土星環(huán)和陰影的土星3：Rx=20°, Ry=0, Rz=25°, lx=0.707, ly=0.707, lz=0

▲圖4 包含土星環(huán)和陰影的土星4：Rx=-10°, Ry=0, Rz=25°, lx=0.707, ly=-0.707, lz=0

▲圖5 包含土星環(huán)和陰影的土星5：Rx=20°, Ry=0, Rz=30°, lx=-1, ly=0, lz=0

為了進(jìn)行比較，你可以在下面網(wǎng)址找到土星的攝像圖：

https://www.jpl.nasa.gov/spaceimages/?search=saturn&category=#submit

圖6所示為構(gòu)建土星環(huán)所用的數(shù)學(xué)模型。這里介紹一種實(shí)現(xiàn)球體著色的算法。首先創(chuàng)建一個(gè)直立球體，也就是說，經(jīng)度是垂直的，緯度是水平的(即平行于XZ平面)，然后從初始方向開始，圍繞x，y和z軸對(duì)球體進(jìn)行旋轉(zhuǎn)。

▲圖6 土星環(huán)模型：行星俯視圖和從XZ平面上Rx=0, Ry=0, Rz=0向下看環(huán)

我們對(duì)土星環(huán)也進(jìn)行同樣的操作。我們可以創(chuàng)建平行于XZ平面的水平環(huán)，然后將它和土星一起旋轉(zhuǎn)相同的角度。土星環(huán)所處的平面穿過土星的球心，因此土星和環(huán)具有相同的旋轉(zhuǎn)中心。

土星環(huán)繪制為一系列相鄰的同心圓，每個(gè)同心圓由短線段組成。參考圖6和文末代碼，程序第42和43行設(shè)置了土星環(huán)的內(nèi)半徑和外半徑，第44行設(shè)置同心圓的間距。土星環(huán)被分成七個(gè)同心環(huán)形帶(圖6中未畫出)且具有不同的顏色，第45行的deltar是它們的寬度。

構(gòu)成同心圓的每個(gè)線段都單獨(dú)繪制。第48行從r1向r2進(jìn)行繪制，通過徑向循環(huán)繪制圓弧段。第49行是繞圓周方向繪制的循環(huán)。第50-61行執(zhí)行旋轉(zhuǎn)操作產(chǎn)生第62和63行中的全局繪圖坐標(biāo)xpg和ypg，旋轉(zhuǎn)函數(shù)與先前程序中的相同。

接下來在第66-75行中設(shè)置線段的顏色。土星環(huán)是由不同顏色的條帶構(gòu)成的，這和NASA觀測(cè)圖像中看到的物理組成結(jié)果一致。從r = r1到r1 + deltar的第一個(gè)條帶具有顏色clr=(.63，.54，.18)，剩余的條帶也是如此。

第五個(gè)條帶省略掉了，因?yàn)樗强盏?，背景顏色能顯示出來。第六個(gè)條帶的寬度是其他條帶的兩倍，并且為第七個(gè)條帶提供了顏色。

對(duì)于給定的光方向，從大多數(shù)角度上，行星體本身都會(huì)在環(huán)上投下陰影。參考圖7，我們的目標(biāo)是確定點(diǎn)p到底位于行星陰影區(qū)域內(nèi)部還是外部。

球狀的行星將產(chǎn)生圓形的陰影，陰影的直徑與行星的尺寸相等，或者更準(zhǔn)確地說，是球體的“大圓”。它是用通過圓心的平面切割球體而得到的最大圓，就像把橙子切成兩半，你看到的是一個(gè)橙子的最大圓。

在圖7中，這種陰影可能是由相同大小的圓盤投影產(chǎn)生的，也可能是由球狀行星投影產(chǎn)生的，兩種情況下，陰影的大小都是一樣的。在土星的側(cè)面圖中，大圓顯示為是一條通過平面中心的加粗線。

從圖7的幾何圖形中可以看出，如果p位于|B| > rs的位置，則它位于陰影區(qū)之外，其中rs是土星的半徑;如果|B| < rs，則p位于在陰影區(qū)之中。在繪制條帶的時(shí)候，如果我們確定了p的位置在陰影區(qū)中，我們就把這個(gè)點(diǎn)涂成灰色，如果它在陰影區(qū)之外，我們就用第66-75行設(shè)置的條帶顏色給它著色。

▲圖7 陰影模型

我們的目標(biāo)是求出給定位置p時(shí)的|B|值，由圖7可看出：

|B|=|V|sin(φ)

并且我們知道：

V×û=|V||û|sin(φ)

其中û=-î 。將上述方程與|û|=1合并得：

B=V×û

|B|=|V×û|

在代碼第78行中確定了入射光矢量î 的長(zhǎng)度為1，但如果在第23-25行中輸入的分量不計(jì)算為1，即

則入射光矢量可能不等于1。有必要的話，需要在第79-81行進(jìn)行重構(gòu)。第82-84行建立矢量V的分量。第85-87行計(jì)算B的分量。第88行給出其幅度magB = |B|。

第89行確定p是否位于陰影區(qū)內(nèi)，如果是，則執(zhí)行第90行V與î的點(diǎn)積。這是用來確定p是否位于行星朝向光源的一側(cè)，在這種情況下，它與行星的暗側(cè)相對(duì)，而不在陰影區(qū)。因?yàn)樵诘?8-89行的陰影算法中并未區(qū)分p的位置，所以此處必須進(jìn)行明確。如果p確實(shí)位于陰影區(qū)域內(nèi)的暗側(cè)，則在第91行中將p設(shè)置為中等灰色。

相信你一定注意到，圖6的土星環(huán)上有一個(gè)暗色條帶，這是因?yàn)橥列黔h(huán)在該位置上是空的：那里沒有固體顆粒物，不能反射陽光。因此我們透過條帶看到了背景顏色'midnightblue'。但這會(huì)產(chǎn)生一個(gè)問題，即陰影顏色的繪制會(huì)覆蓋該空白處的背景顏色，因此在第93和94行將其重置為'midnightblue'。

既然條帶的顏色都建立好了，我們就可以通過一個(gè)個(gè)短線段來繪制土星環(huán)了。第97-100行計(jì)算第一個(gè)線段的起始位置。參考圖6，第100-101行確定該線段與土星的遮擋關(guān)系，線段在前就會(huì)被繪制。

103-108行確定線段是否在土星后面，位于后面就不會(huì)被繪制。這種遮擋關(guān)系是通過計(jì)算點(diǎn)的全局坐標(biāo)與土星中心之間的距離c來完成的。

第107行的意思是，如果c大于球體半徑的1.075倍，則繪制該段。因子1.075的作用是防止線段與球體的表面重合，這是有必要的，不然小于球體半徑的可見區(qū)段將不會(huì)被繪制。

本文代碼生成的圖像有兩點(diǎn)需要注意：首先是顏色。美國(guó)宇航局的攝影圖像呈現(xiàn)的是一種幾乎沒有顏色灰色，但是許多土星觀察者都將它描述為金黃色，因此我們選擇了金色。

所有攝影師都知道，在攝影圖像中呈現(xiàn)物體的真實(shí)顏色是十分困難的，顏色取決于入射光和物體本身的顏色，或許最好的方法是依靠肉眼觀察。

如果你不贊同本文代碼所生成的圖像中的顏色，可以通過更改程序中clr的定義來修改它們。需要注意的第二點(diǎn)，是圖5中土星表面陰影的曲率，它表示著色算法是否按預(yù)期工作。

在程序的使用上，你可以自行更改第24-26行中的入射光的方向和第32-34行中的旋轉(zhuǎn)角度。

本文篇幅有限，更多更詳細(xì)的講解請(qǐng)參閱《Python圖形編程：2D和3D圖像的創(chuàng)建》一書。

• 土星代碼

運(yùn)行代碼也需要一段時(shí)間，請(qǐng)耐心等待。

1""" 
  2SATURN 
  3""" 
  4 
  5import numpy as np 
  6import matplotlib.pyplot as plt 
  7from math import sin, cos, radians, sqrt 
  8 
  9plt.axis([0,150,100,0]) 
 10plt.axis('off') 
 11plt.grid(False) 
 12 
 13print('running') 
 14#—————————————————parameters 
 15g=[0]*3 
 16 
 17xc=80 #———sphere center 
 18yc=50 
 19zc=0 
 20 
 21rs=25 #———sphere radius 
 22 
 23lx=-1 #———light ray unit vector components 
 24ly=0 
 25lz=0 
 26 
 27IA=0 
 28IB=.8 
 29+n=2 
 30 
 31Rx=radians(-20) 
 32Ry=radians(0) 
 33Rz=radians(30) 
 34 
 35#————————same as SHADESPHERE—————– 
 36 
 37#———————————————————rings 
 38alpha1=radians(-10) 
 39alpha2=radians(370) 
 40dalpha=radians(.5) 
 41 
 42r1=rs*1.5 
 43r2=rs*2.2 
 44dr=rs*.02 
 45deltar=(r2-r1)/7 #———ring band width 
 46 
 47#—————————————rotate ring point p which is at r, alpha 
 48for 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 
112plt.show()