手把手教你用Python画一个绝美土星环

开发 后端
土星的行星环非常出名。虽然木星、土星、天王星和海王星也有环,但土星环是我们太阳系中最大、最亮、最广为人知的行星环。

 土星的行星环非常出名。虽然木星、土星、天王星和海王星也有环,但土星环是我们太阳系中最大、最亮、最广为人知的行星环。

它由小到灰尘的颗粒,大到巨石的物体组成。这些物体的成分主要是冰,一般认为是彗星或较大的小行星与土星的一颗卫星相撞时产生的,两者都撞成了小碎块。在远古时代,土星就已为人所知,但直到1610年,伽利略才首次用望远镜对它进行观测。

这个行星以罗马农业之神土星Saturn命名为,也就是我们每个星期的第六天Saturday。

图1至图5中的图像是由本文文末的代码生成的。每张图都呈现不同的方向角,图标题中有相应的说明。

图标题中还列出了入射光线的单位矢量分量,例如lx=+0.707,ly=+0.707,lz=0 表示左上象限中的光源;lx=-1,ly=0,lz=0表示来自右侧的光源。在图像中请注意行星在环上投射的阴影,尤其是在图5中能够看到行星轮廓的曲率。

 

手把手教你用Python画一个绝美土星环

 

▲图1 包含土星环和阴影的土星1:Rx=-20°, Ry=0, Rz=-10°, lx=1, ly=0, lz=0

 

手把手教你用Python画一个绝美土星环

 

▲图2 包含土星环和阴影的土星2:Rx=-8°, Ry=0, Rz=-30°, lx=0.707, ly=.707, lz=0

 

手把手教你用Python画一个绝美土星环

 

▲图3 包含土星环和阴影的土星3:Rx=20°, Ry=0, Rz=25°, lx=0.707, ly=0.707, lz=0

 

手把手教你用Python画一个绝美土星环

 

▲图4 包含土星环和阴影的土星4:Rx=-10°, Ry=0, Rz=25°, lx=0.707, ly=-0.707, lz=0

 

手把手教你用Python画一个绝美土星环

 

▲图5 包含土星环和阴影的土星5:Rx=20°, Ry=0, Rz=30°, lx=-1, ly=0, lz=0

为了进行比较,你可以在下面网址找到土星的摄像图:

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

图6所示为构建土星环所用的数学模型。这里介绍一种实现球体着色的算法。首先创建一个直立球体,也就是说,经度是垂直的,纬度是水平的(即平行于XZ平面),然后从初始方向开始,围绕x,y和z轴对球体进行旋转。

 

手把手教你用Python画一个绝美土星环

 

▲图6 土星环模型:行星俯视图和从XZ平面上Rx=0, Ry=0, Rz=0向下看环

我们对土星环也进行同样的操作。我们可以创建平行于XZ平面的水平环,然后将它和土星一起旋转相同的角度。土星环所处的平面穿过土星的球心,因此土星和环具有相同的旋转中心。

土星环绘制为一系列相邻的同心圆,每个同心圆由短线段组成。参考图6和文末代码,程序第42和43行设置了土星环的内半径和外半径,第44行设置同心圆的间距。土星环被分成七个同心环形带(图6中未画出)且具有不同的颜色,第45行的deltar是它们的宽度。

构成同心圆的每个线段都单独绘制。第48行从r1向r2进行绘制,通过径向循环绘制圆弧段。第49行是绕圆周方向绘制的循环。第50-61行执行旋转操作产生第62和63行中的全局绘图坐标xpg和ypg,旋转函数与先前程序中的相同。

接下来在第66-75行中设置线段的颜色。土星环是由不同颜色的条带构成的,这和NASA观测图像中看到的物理组成结果一致。从r = r1到r1 + deltar的第一个条带具有颜色clr=(.63,.54,.18),剩余的条带也是如此。

第五个条带省略掉了,因为它是空的,背景颜色能显示出来。第六个条带的宽度是其他条带的两倍,并且为第七个条带提供了颜色。

对于给定的光方向,从大多数角度上,行星体本身都会在环上投下阴影。参考图7,我们的目标是确定点p到底位于行星阴影区域内部还是外部。

球状的行星将产生圆形的阴影,阴影的直径与行星的尺寸相等,或者更准确地说,是球体的“大圆”。它是用通过圆心的平面切割球体而得到的最大圆,就像把橙子切成两半,你看到的是一个橙子的最大圆。

在图7中,这种阴影可能是由相同大小的圆盘投影产生的,也可能是由球状行星投影产生的,两种情况下,阴影的大小都是一样的。在土星的侧面图中,大圆显示为是一条通过平面中心的加粗线。

从图7的几何图形中可以看出,如果p位于|B| > rs的位置,则它位于阴影区之外,其中rs是土星的半径;如果|B| < rs,则p位于在阴影区之中。在绘制条带的时候,如果我们确定了p的位置在阴影区中,我们就把这个点涂成灰色,如果它在阴影区之外,我们就用第66-75行设置的条带颜色给它着色。

 

手把手教你用Python画一个绝美土星环

 

▲图7 阴影模型

我们的目标是求出给定位置p时的|B|值,由图7可看出:

|B|=|V|sin(φ)

并且我们知道:

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

其中û=-î 。将上述方程与|û|=1合并得:

B=V×û

|B|=|V×û|

在代码第78行中确定了入射光矢量î 的长度为1,但如果在第23-25行中输入的分量不计算为1,即

 

手把手教你用Python画一个绝美土星环

 

则入射光矢量可能不等于1。有必要的话,需要在第79-81行进行重构。第82-84行建立矢量V的分量。第85-87行计算B的分量。第88行给出其幅度magB = |B|。

第89行确定p是否位于阴影区内,如果是,则执行第90行V与î的点积。这是用来确定p是否位于行星朝向光源的一侧,在这种情况下,它与行星的暗侧相对,而不在阴影区。因为在第78-89行的阴影算法中并未区分p的位置,所以此处必须进行明确。如果p确实位于阴影区域内的暗侧,则在第91行中将p设置为中等灰色。

相信你一定注意到,图6的土星环上有一个暗色条带,这是因为土星环在该位置上是空的:那里没有固体颗粒物,不能反射阳光。因此我们透过条带看到了背景颜色'midnightblue'。但这会产生一个问题,即阴影颜色的绘制会覆盖该空白处的背景颜色,因此在第93和94行将其重置为'midnightblue'。

既然条带的颜色都建立好了,我们就可以通过一个个短线段来绘制土星环了。第97-100行计算第一个线段的起始位置。参考图6,第100-101行确定该线段与土星的遮挡关系,线段在前就会被绘制。

103-108行确定线段是否在土星后面,位于后面就不会被绘制。这种遮挡关系是通过计算点的全局坐标与土星中心之间的距离c来完成的。

第107行的意思是,如果c大于球体半径的1.075倍,则绘制该段。因子1.075的作用是防止线段与球体的表面重合,这是有必要的,不然小于球体半径的可见区段将不会被绘制。

本文代码生成的图像有两点需要注意:首先是颜色。美国宇航局的摄影图像呈现的是一种几乎没有颜色灰色,但是许多土星观察者都将它描述为金黄色,因此我们选择了金色。

所有摄影师都知道,在摄影图像中呈现物体的真实颜色是十分困难的,颜色取决于入射光和物体本身的颜色,或许最好的方法是依靠肉眼观察。

如果你不赞同本文代码所生成的图像中的颜色,可以通过更改程序中clr的定义来修改它们。需要注意的第二点,是图5中土星表面阴影的曲率,它表示着色算法是否按预期工作。

在程序的使用上,你可以自行更改第24-26行中的入射光的方向和第32-34行中的旋转角度。

本文篇幅有限,更多更详细的讲解请参阅《Python图形编程:2D和3D图像的创建》一书。

  • 土星代码

运行代码也需要一段时间,请耐心等待。

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() 
  • 1.
  • 2.
  • 3.
  • 4.
  • 5.
  • 6.
  • 7.
  • 8.
  • 9.
  • 10.
  • 11.
  • 12.
  • 13.
  • 14.
  • 15.
  • 16.
  • 17.
  • 18.
  • 19.
  • 20.
  • 21.
  • 22.
  • 23.
  • 24.
  • 25.
  • 26.
  • 27.
  • 28.
  • 29.
  • 30.
  • 31.
  • 32.
  • 33.
  • 34.
  • 35.
  • 36.
  • 37.
  • 38.
  • 39.
  • 40.
  • 41.
  • 42.
  • 43.
  • 44.
  • 45.
  • 46.
  • 47.
  • 48.
  • 49.
  • 50.
  • 51.
  • 52.
  • 53.
  • 54.
  • 55.
  • 56.
  • 57.
  • 58.
  • 59.
  • 60.
  • 61.
  • 62.
  • 63.
  • 64.
  • 65.
  • 66.
  • 67.
  • 68.
  • 69.
  • 70.
  • 71.
  • 72.
  • 73.
  • 74.
  • 75.
  • 76.
  • 77.
  • 78.
  • 79.
  • 80.
  • 81.
  • 82.
  • 83.
  • 84.
  • 85.
  • 86.
  • 87.
  • 88.
  • 89.
  • 90.
  • 91.
  • 92.
  • 93.
  • 94.
  • 95.
  • 96.
  • 97.
  • 98.
  • 99.
  • 100.
  • 101.
  • 102.
  • 103.
  • 104.
  • 105.
  • 106.
  • 107.
  • 108.
  • 109.
  • 110.
  • 111.
  • 112.

 

责任编辑:华轩 来源: 今日头条
相关推荐

2021-08-09 13:31:25

PythonExcel代码

2021-07-12 09:03:50

Python任务管理器cmd命令

2022-10-19 14:30:59

2021-07-12 14:35:26

代码架构云原生

2021-05-10 06:48:11

Python腾讯招聘

2021-12-11 20:20:19

Python算法线性

2021-02-02 13:31:35

Pycharm系统技巧Python

2020-05-26 10:20:56

Python开发工具

2011-03-28 16:14:38

jQuery

2021-02-06 14:55:05

大数据pandas数据分析

2021-02-04 09:00:57

SQLDjango原生

2022-02-25 09:41:05

python搜索引擎

2021-08-24 10:02:21

JavaScript网页搜索 前端

2024-11-05 16:40:24

JavaScript搜索引擎

2022-08-04 10:39:23

Jenkins集成CD

2021-05-17 21:30:06

Python求均值中值

2021-01-21 09:10:29

ECharts柱状图大数据

2021-01-08 10:32:24

Charts折线图数据可视化

2009-04-22 09:17:19

LINQSQL基础

2012-01-11 13:40:35

移动应用云服务
点赞
收藏

51CTO技术栈公众号