溫馨提示×
您好，登錄后才能下訂單哦！
點(diǎn)擊 登錄注冊(cè) 即表示同意《億速云用戶服務(wù)條款》
您好，登錄后才能下訂單哦！
本篇文章為大家展示了如何用Python畫一個(gè)絕美土星環(huán)，內(nèi)容簡(jiǎn)明扼要并且容易理解，絕對(duì)能使你眼前一亮，通過這篇文章的詳細(xì)介紹希望你能有所收獲。
土星的行星環(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)角度。
運(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()上述內(nèi)容就是如何用Python畫一個(gè)絕美土星環(huán)，你們學(xué)到知識(shí)或技能了嗎？如果還想學(xué)到更多技能或者豐富自己的知識(shí)儲(chǔ)備，歡迎關(guān)注億速云行業(yè)資訊頻道。
免責(zé)聲明：本站發(fā)布的內(nèi)容（圖片、視頻和文字）以原創(chuàng)、轉(zhuǎn)載和分享為主，文章觀點(diǎn)不代表本網(wǎng)站立場(chǎng)，如果涉及侵權(quán)請(qǐng)聯(lián)系站長(zhǎng)郵箱：is@yisu.com進(jìn)行舉報(bào)，并提供相關(guān)證據(jù)，一經(jīng)查實(shí)，將立刻刪除涉嫌侵權(quán)內(nèi)容。
億速云公眾號(hào)
手機(jī)網(wǎng)站二維碼