layout: post
title: 利用PyMOL制作反应势能面示意动画
categories:
- 科
tags:
- 3d
- pymol
- python
math: true
- 2013-04-18 10:43:35 发布图片
- 2013-11-14 21:18:18 增加教程
- 2016-09-07 20:17:18 颜色映射
利用PyMOL内置的OpenGL制作了势能面反应的示意图, 具体制作过程容后再讲.
PyMOL是基于Python的分子可视化软件, 在结构生物学中使用十分广泛. 它的主要特色是对生物大分子显示效果好, 并自带一个高效的光线追踪渲染器, 能渲染出逼真的效果. 此外, PyMOL还支持脚本, 可用于精确控制显示, 并提供了一个基本的OpenGL的接口, 以便用户进行一些简单的三维动画设计.
三维建模软件有很多, 不同的软件适用于不同的领域. 最基础的如OpenGL之类, 需要你利用最基本的图元和场景选项一点一点地构建出整个场景. POV-Ray之类则侧重于光线追踪渲染, 并集成了许多可用的模型与场景, 更高级一些. Blender之类则是更通用更专业的三维建模软件, 并非专用于分子建模与可视化. 这样看起来, PyMOL算是处于OpenGL和Blender中间, 集成了POV-Ray的部分功能, 专门用于分子可视化的. 建模时它可以实时可视化, 调试很分方便, 因此对于一些简单的应用很合适.
PyMOL中, 三维模型被称作CGO(Compiled Graphics Objects), 可在其中引用一些OpenGL基本图元. 使用方法和OpenGL很类似, 但由于调用是基于Python的, 所以比直接使用OpenGL简单一些.
下面说说上面两个动画的制作方法.
首先我们需要一个二维势能面的模型. 这个模型最好能够具有一般势能面的特点, 并且包含各类驻点, 如极大点, 极小点, 鞍点. 这就要求势能面上至少有两个极大点. 经过比较, 我发现二元函数 $F(x,y)={\sin x \over x}+{\sin y \over y}$ 满足要求. 但此函数为周期函数, 所以最好加上一个线性项去掉周期性, 并适当调整大小. 最终我用的函数是 $F(x,y)=4({\sin x \over x}+{\sin y \over y})+0.1x$.
确定了势能函数的解析式就可以创建势能曲面了. 方法是最基本的剖分, 用三角面片对整个区间进行剖分, 再将剖分所得的三角面片组合起来. 值得注意的是剖分的方向和三角面片的法向.
至于曲面上球体运动的模拟, 可利用简单的线性步长方法. 若需要更加真实的效果, 则可利用分子动力学中的Verlet积分方法进行计算.
使用方法
- 运行代码:
run PESdemo.py
- 背景设为白色:
bg white
- 设定光线追踪:
set ray_trace_frames=1
- 输出png图片:
mpng PESdemo
渲染好的图片将输出为PESdemoXXXX.png, XXXX为编号. 利用这些图片就可以制作成动画.
支持动画的图片格式目前主要两种:
gif 最常用的, 大家也最熟悉. 各种软件支持最好, 算是通用格式. 可惜只有256色, 失真有时很严重, 特别是对光线追踪渲染过的图片.
apng 基于png的动画图片, 效果等同于png. 可惜目前浏览器支持不广, 尚未得到png官方承认.
下面两个图就是这两种格式的对比
更具体的信息可参考下面的网文:
- 小牛犊APNG力挫老古董MNG
- APNG编辑制作工具
如果你要使用其他颜色映射方案的话, 可参考我的另一篇博文几种颜色映射方案的解析式.
代码
<table class="highlighttable"><th colspan="2">PESdemo.py</th><tr><td><div class="linenodiv" style="background-color: #f0f0f0; padding-right: 10px"><pre style="line-height: 125%"> 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
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188</pre></div></td><td class="code"><div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span style="color: #008800; font-style: italic"># coding: utf-8</span>
<span style="color: #008800; font-style: italic"># ##############################################################################</span>
<span style="color: #008800; font-style: italic"># 2013-11-14 10:25:11 简单示例</span>
<span style="color: #008800; font-style: italic"># 2016-09-07 20:23:18 注释, 颜色映射</span>
<span style="color: #008800; font-style: italic"># ##############################################################################</span>
<span style="color: #AA22FF; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">math</span>
<span style="color: #AA22FF; font-weight: bold">from</span> <span style="color: #0000FF; font-weight: bold">pymol</span> <span style="color: #AA22FF; font-weight: bold">import</span> cmd
<span style="color: #AA22FF; font-weight: bold">from</span> <span style="color: #0000FF; font-weight: bold">pymol.cgo</span> <span style="color: #AA22FF; font-weight: bold">import</span> <span style="color: #666666">*</span>
<span style="color: #008800; font-style: italic"># 势能面函数</span>
<span style="color: #AA22FF; font-weight: bold">def</span> <span style="color: #00A000">Fxy</span>(x,y):
<span style="color: #AA22FF; font-weight: bold">if</span> x<span style="color: #666666">==0</span>: Fx <span style="color: #666666">=</span> <span style="color: #666666">4</span>
<span style="color: #AA22FF; font-weight: bold">else</span>: Fx <span style="color: #666666">=</span> <span style="color: #666666">4*</span>math<span style="color: #666666">.</span>sin(x)<span style="color: #666666">/</span>x<span style="color: #666666">+.1*</span>x
<span style="color: #AA22FF; font-weight: bold">if</span> y<span style="color: #666666">==0</span>: Fy <span style="color: #666666">=</span> <span style="color: #666666">4</span>
<span style="color: #AA22FF; font-weight: bold">else</span>: Fy <span style="color: #666666">=</span> <span style="color: #666666">4*</span>math<span style="color: #666666">.</span>sin(y)<span style="color: #666666">/</span>y
<span style="color: #AA22FF; font-weight: bold">return</span> Fx<span style="color: #666666">+</span>Fy
<span style="color: #008800; font-style: italic"># 势能面函数法向</span>
<span style="color: #AA22FF; font-weight: bold">def</span> <span style="color: #00A000">dFxy</span>(x,y):
<span style="color: #AA22FF; font-weight: bold">if</span> x<span style="color: #666666">==0</span>: dFx <span style="color: #666666">=</span> <span style="color: #666666">0</span>
<span style="color: #AA22FF; font-weight: bold">else</span>: dFx <span style="color: #666666">=</span> <span style="color: #666666">4*</span>(x<span style="color: #666666">*</span>math<span style="color: #666666">.</span>cos(x)<span style="color: #666666">-</span>math<span style="color: #666666">.</span>sin(x))<span style="color: #666666">/</span>x<span style="color: #666666">**2+0.1</span>
<span style="color: #AA22FF; font-weight: bold">if</span> y<span style="color: #666666">==0</span>: dFy <span style="color: #666666">=</span> <span style="color: #666666">0</span>
<span style="color: #AA22FF; font-weight: bold">else</span>: dFy <span style="color: #666666">=</span> <span style="color: #666666">4*</span>(y<span style="color: #666666">*</span>math<span style="color: #666666">.</span>cos(y)<span style="color: #666666">-</span>math<span style="color: #666666">.</span>sin(y))<span style="color: #666666">/</span>y<span style="color: #666666">**2</span>
Rtmp <span style="color: #666666">=</span> <span style="color: #666666">1/</span>math<span style="color: #666666">.</span>sqrt(dFx<span style="color: #666666">*</span>dFx<span style="color: #666666">+</span>dFy<span style="color: #666666">*</span>dFy<span style="color: #666666">+1</span>)
<span style="color: #AA22FF; font-weight: bold">return</span> <span style="color: #666666">-</span>dFx<span style="color: #666666">*</span>Rtmp, <span style="color: #666666">-</span>dFy<span style="color: #666666">*</span>Rtmp, Rtmp
<span style="color: #AA22FF; font-weight: bold">def</span> <span style="color: #00A000">RGB</span>(V, Vmin, Vmax):
dV<span style="color: #666666">=</span>Vmax<span style="color: #666666">-</span>Vmin; x<span style="color: #666666">=</span>(V<span style="color: #666666">-</span>Vmin)<span style="color: #666666">/</span>dV
r<span style="color: #666666">=1</span>; g<span style="color: #666666">=1</span>; b<span style="color: #666666">=1</span>
<span style="color: #AA22FF; font-weight: bold">if</span> x<span style="color: #666666"><0.25</span>: r <span style="color: #666666">=</span> <span style="color: #666666">0</span>; g <span style="color: #666666">=</span> <span style="color: #666666">4*</span>x
<span style="color: #AA22FF; font-weight: bold">elif</span> x<span style="color: #666666"><0.50</span>: r <span style="color: #666666">=</span> <span style="color: #666666">0</span>; b <span style="color: #666666">=</span> <span style="color: #666666">2-4*</span>x
<span style="color: #AA22FF; font-weight: bold">elif</span> x<span style="color: #666666"><0.75</span>: r <span style="color: #666666">=</span> <span style="color: #666666">4*</span>x<span style="color: #666666">-2</span>; b <span style="color: #666666">=</span> <span style="color: #666666">0</span>
<span style="color: #AA22FF; font-weight: bold">else</span>: g <span style="color: #666666">=</span> <span style="color: #666666">4-4*</span>x; b <span style="color: #666666">=</span> <span style="color: #666666">0</span>
r<span style="color: #666666">=</span><span style="color: #AA22FF">min</span>(<span style="color: #AA22FF">max</span>(r,<span style="color: #666666">0</span>), <span style="color: #666666">1</span>)
g<span style="color: #666666">=</span><span style="color: #AA22FF">min</span>(<span style="color: #AA22FF">max</span>(g,<span style="color: #666666">0</span>) ,<span style="color: #666666">1</span>)
b<span style="color: #666666">=</span><span style="color: #AA22FF">min</span>(<span style="color: #AA22FF">max</span>(b,<span style="color: #666666">0</span>) ,<span style="color: #666666">1</span>)
<span style="color: #AA22FF; font-weight: bold">return</span> r, g, b
<span style="color: #008800; font-style: italic"># 是否使用法向, 显示网格, 颜色映射, 使用Verlet方法计算轨迹</span>
YesNorm<span style="color: #666666">=1</span>; YesGrid<span style="color: #666666">=1</span>; YesMap<span style="color: #666666">=0</span>; YesTrj<span style="color: #666666">=1</span>
Xini <span style="color: #666666">=</span> <span style="color: #666666">4.5</span> <span style="color: #008800; font-style: italic"># 小球初始位置</span>
Xmin <span style="color: #666666">=</span> <span style="color: #666666">-7</span>; Xmax <span style="color: #666666">=</span> <span style="color: #666666">5.5</span>; dX <span style="color: #666666">=</span> <span style="color: #666666">.5</span>; Nx <span style="color: #666666">=</span> <span style="color: #AA22FF">int</span>((Xmax<span style="color: #666666">-</span>Xmin)<span style="color: #666666">/</span>dX)
Ymin <span style="color: #666666">=</span> <span style="color: #666666">-7</span>; Ymax <span style="color: #666666">=</span> <span style="color: #666666">7</span>; dY <span style="color: #666666">=</span> <span style="color: #666666">.5</span>; Ny <span style="color: #666666">=</span> <span style="color: #AA22FF">int</span>((Ymax<span style="color: #666666">-</span>Ymin)<span style="color: #666666">/</span>dY)
X <span style="color: #666666">=</span> [ <span style="color: #666666">0</span> <span style="color: #AA22FF; font-weight: bold">for</span> i <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #AA22FF">range</span>(Nx) ]
Y <span style="color: #666666">=</span> [ <span style="color: #666666">0</span> <span style="color: #AA22FF; font-weight: bold">for</span> j <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #AA22FF">range</span>(Ny) ]
Z <span style="color: #666666">=</span> [ [<span style="color: #666666">0</span>]<span style="color: #666666">*</span>Ny <span style="color: #AA22FF; font-weight: bold">for</span> i <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #AA22FF">range</span>(Nx) ]
Zx<span style="color: #666666">=</span> [ [<span style="color: #666666">0</span>]<span style="color: #666666">*</span>Ny <span style="color: #AA22FF; font-weight: bold">for</span> i <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #AA22FF">range</span>(Nx) ]
Zy<span style="color: #666666">=</span> [ [<span style="color: #666666">0</span>]<span style="color: #666666">*</span>Ny <span style="color: #AA22FF; font-weight: bold">for</span> i <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #AA22FF">range</span>(Nx) ]
Zz<span style="color: #666666">=</span> [ [<span style="color: #666666">1</span>]<span style="color: #666666">*</span>Ny <span style="color: #AA22FF; font-weight: bold">for</span> i <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #AA22FF">range</span>(Nx) ]
<span style="color: #AA22FF; font-weight: bold">for</span> i <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #AA22FF">range</span>(Nx): X[i] <span style="color: #666666">=</span> Xmin<span style="color: #666666">+</span>dX<span style="color: #666666">*</span>i
<span style="color: #AA22FF; font-weight: bold">for</span> j <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #AA22FF">range</span>(Ny): Y[j] <span style="color: #666666">=</span> Ymin<span style="color: #666666">+</span>dY<span style="color: #666666">*</span>j
<span style="color: #AA22FF; font-weight: bold">for</span> i <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #AA22FF">range</span>(Nx):
x<span style="color: #666666">=</span>X[i]
<span style="color: #AA22FF; font-weight: bold">for</span> j <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #AA22FF">range</span>(Ny):
y<span style="color: #666666">=</span>Y[j]
Z[i][j] <span style="color: #666666">=</span> Fxy(x,y)
<span style="color: #AA22FF; font-weight: bold">if</span> YesNorm: Zx[i][j], Zy[i][j], Zz[i][j] <span style="color: #666666">=</span> dFxy(x,y)
<span style="color: #008800; font-style: italic"># 获取极值, 用于颜色映射</span>
Zmin<span style="color: #666666">=</span><span style="color: #AA22FF">min</span>(<span style="color: #AA22FF">min</span>(Z)); Zmax<span style="color: #666666">=</span><span style="color: #AA22FF">max</span>(<span style="color: #AA22FF">max</span>(Z))
PES <span style="color: #666666">=</span> []
<span style="color: #008800; font-style: italic"># 绘制网格</span>
<span style="color: #AA22FF; font-weight: bold">if</span> YesGrid:
PES<span style="color: #666666">.</span>extend( [ COLOR, <span style="color: #666666">0</span>, <span style="color: #666666">0</span>, <span style="color: #666666">0</span> ] )
<span style="color: #AA22FF; font-weight: bold">for</span> i <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #AA22FF">range</span>(Nx):
PES<span style="color: #666666">.</span>extend( [ BEGIN, LINE_STRIP ] )
<span style="color: #AA22FF; font-weight: bold">for</span> j <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #AA22FF">range</span>(<span style="color: #666666">0</span>,Ny):
PES<span style="color: #666666">.</span>extend( [ NORMAL, Zx[i][j], Zy[i][j], Zz[i][j] ] )
PES<span style="color: #666666">.</span>extend( [ VERTEX, X[i], Y[j], Z[i][j] ] )
PES<span style="color: #666666">.</span>append( END )
<span style="color: #AA22FF; font-weight: bold">for</span> j <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #AA22FF">range</span>(Ny):
PES<span style="color: #666666">.</span>extend( [ BEGIN, LINE_STRIP ] )
<span style="color: #AA22FF; font-weight: bold">for</span> i <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #AA22FF">range</span>(Nx):
PES<span style="color: #666666">.</span>extend( [ NORMAL, Zx[i][j], Zy[i][j], Zz[i][j] ] )
PES<span style="color: #666666">.</span>extend( [ VERTEX, X[i], Y[j], Z[i][j] ] )
PES<span style="color: #666666">.</span>append( END )
<span style="color: #008800; font-style: italic"># 绘制表面</span>
PES<span style="color: #666666">.</span>extend( [ COLOR, <span style="color: #666666">.19</span>, <span style="color: #666666">.6</span>, <span style="color: #666666">.83</span> ] )
<span style="color: #AA22FF; font-weight: bold">for</span> j <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #AA22FF">range</span>(Ny<span style="color: #666666">-1</span>):
PES<span style="color: #666666">.</span>extend( [ BEGIN, TRIANGLE_STRIP ] )
<span style="color: #AA22FF; font-weight: bold">for</span> i <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #AA22FF">range</span>(Nx):
<span style="color: #AA22FF; font-weight: bold">if</span> YesMap:
r, g, b<span style="color: #666666">=</span>RGB(Z[i][j<span style="color: #666666">+1</span>], Zmin, Zmax)
PES<span style="color: #666666">.</span>extend( [ COLOR, r, g, b ] )
PES<span style="color: #666666">.</span>extend( [ NORMAL, Zx[i][j<span style="color: #666666">+1</span>], Zy[i][j<span style="color: #666666">+1</span>], Zz[i][j<span style="color: #666666">+1</span>] ] )
PES<span style="color: #666666">.</span>extend( [ VERTEX, X[i], Y[j<span style="color: #666666">+1</span>], Z[i][j<span style="color: #666666">+1</span>] ] )
<span style="color: #AA22FF; font-weight: bold">if</span> YesMap:
r, g, b<span style="color: #666666">=</span>RGB(Z[i][j], Zmin, Zmax)
PES<span style="color: #666666">.</span>extend( [ COLOR, r, g, b ] )
PES<span style="color: #666666">.</span>extend( [ NORMAL, Zx[i][j], Zy[i][j], Zz[i][j] ] )
PES<span style="color: #666666">.</span>extend( [ VERTEX, X[i], Y[j], Z[i][j] ] )
PES<span style="color: #666666">.</span>append( END )
<span style="color: #008800; font-style: italic"># 绘制路径</span>
PES<span style="color: #666666">.</span>extend( [ LINEWIDTH, <span style="color: #666666">5</span> ] )
PES<span style="color: #666666">.</span>extend( [ BEGIN, LINE_STRIP ] )
PES<span style="color: #666666">.</span>extend( [ COLOR, <span style="color: #666666">1.</span>, <span style="color: #666666">0.</span>, <span style="color: #666666">.0</span> ] )
x <span style="color: #666666">=</span> <span style="color: #666666">-</span>Xini; y <span style="color: #666666">=</span> <span style="color: #666666">-</span>Xini
<span style="color: #AA22FF; font-weight: bold">while</span> y<span style="color: #666666"><=</span>Xini:
z <span style="color: #666666">=</span> Fxy(x,y)
dFx, dFy, dFz <span style="color: #666666">=</span> dFxy(x,y)
PES<span style="color: #666666">.</span>extend( [ NORMAL, dFx, dFy, dFz ] )
PES<span style="color: #666666">.</span>extend( [ VERTEX, x, y, z ] )
y <span style="color: #666666">=</span> y<span style="color: #666666">+0.1</span>
PES<span style="color: #666666">.</span>append( END )
PES<span style="color: #666666">.</span>extend( [ BEGIN, LINE_STRIP] )
PES<span style="color: #666666">.</span>extend( [ COLOR, <span style="color: #666666">0</span>, <span style="color: #666666">1.</span>, <span style="color: #666666">.0</span> ] )
x <span style="color: #666666">=</span> <span style="color: #666666">-</span>Xini; y <span style="color: #666666">=</span> <span style="color: #666666">-</span>Xini
<span style="color: #AA22FF; font-weight: bold">while</span> x<span style="color: #666666"><=</span>Xini:
z <span style="color: #666666">=</span> Fxy(x,y)
dFx, dFy, dFz <span style="color: #666666">=</span> dFxy(x,y)
PES<span style="color: #666666">.</span>extend( [ NORMAL, dFx, dFy, dFz ] )
PES<span style="color: #666666">.</span>extend( [ VERTEX, x, y, z ] )
x <span style="color: #666666">=</span> x<span style="color: #666666">+0.1</span>
PES<span style="color: #666666">.</span>append( END )
<span style="color: #008800; font-style: italic"># 绘制小球</span>
x <span style="color: #666666">=</span> <span style="color: #666666">-</span>Xini; y <span style="color: #666666">=</span> <span style="color: #666666">-</span>Xini; z <span style="color: #666666">=</span> Fxy(x,y)
PES<span style="color: #666666">.</span>extend ( [ COLOR, <span style="color: #666666">1</span>, <span style="color: #666666">1</span>, <span style="color: #666666">0</span> ] )
PES<span style="color: #666666">.</span>extend ( [ SPHERE, x, y, z, <span style="color: #666666">.2</span> ] )
x <span style="color: #666666">=</span> <span style="color: #666666">-</span>Xini; y <span style="color: #666666">=</span> Xini; z <span style="color: #666666">=</span> Fxy(x,y)
PES<span style="color: #666666">.</span>extend ( [ COLOR, <span style="color: #666666">1</span>, <span style="color: #666666">0</span>, <span style="color: #666666">0</span> ] )
PES<span style="color: #666666">.</span>extend ( [ SPHERE, x, y, z, <span style="color: #666666">.2</span> ] )
x <span style="color: #666666">=</span> Xini; y <span style="color: #666666">=</span> <span style="color: #666666">-</span>Xini; z <span style="color: #666666">=</span> Fxy(x,y)
PES<span style="color: #666666">.</span>extend ( [ COLOR, <span style="color: #666666">0</span>, <span style="color: #666666">1</span>, <span style="color: #666666">0</span> ] )
PES<span style="color: #666666">.</span>extend ( [ SPHERE, x, y, z, <span style="color: #666666">.2</span> ] )
cmd<span style="color: #666666">.</span>load_cgo(PES, <span style="color: #BB4444">'</span>PES<span style="color: #BB4444">'</span>)
<span style="color: #008800; font-style: italic"># 模拟小球运行</span>
Rsph <span style="color: #666666">=</span> <span style="color: #666666">0.5</span>
<span style="color: #AA22FF; font-weight: bold">if</span> <span style="color: #AA22FF; font-weight: bold">not</span> YesTrj:
x <span style="color: #666666">=</span> <span style="color: #666666">-</span>Xini
<span style="color: #AA22FF; font-weight: bold">for</span> Ifrm <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #AA22FF">range</span>(<span style="color: #666666">30</span>):
y <span style="color: #666666">=</span> <span style="color: #666666">-</span>Xini<span style="color: #666666">+</span>Ifrm<span style="color: #666666">*0.3</span>; z <span style="color: #666666">=</span> Fxy(x,y)
dFx, dFy, dFz <span style="color: #666666">=</span> dFxy(x,y)
SYS <span style="color: #666666">=</span> [ COLOR, <span style="color: #666666">1</span>, <span style="color: #666666">1-</span>Ifrm<span style="color: #666666">/30.</span>, <span style="color: #666666">0</span> ]
SYS<span style="color: #666666">.</span>extend( [ SPHERE, x<span style="color: #666666">+</span>Rsph<span style="color: #666666">*</span>dFx, y<span style="color: #666666">+</span>Rsph<span style="color: #666666">*</span>dFy, z<span style="color: #666666">+</span>Rsph<span style="color: #666666">*</span>dFz, Rsph ] )
cmd<span style="color: #666666">.</span>load_cgo(SYS, <span style="color: #BB4444">'</span>SYS<span style="color: #BB4444">'</span>, Ifrm)
y <span style="color: #666666">=</span> <span style="color: #666666">-</span>Xini
<span style="color: #AA22FF; font-weight: bold">for</span> Ifrm <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #AA22FF">range</span>(<span style="color: #666666">30</span>):
x <span style="color: #666666">=</span> <span style="color: #666666">-</span>Xini<span style="color: #666666">+</span>Ifrm<span style="color: #666666">*0.3</span>; z <span style="color: #666666">=</span> Fxy(x,y)
dFx, dFy, dFz <span style="color: #666666">=</span> dFxy(x,y)
SYS <span style="color: #666666">=</span> [ COLOR, <span style="color: #666666">1-</span>Ifrm<span style="color: #666666">/30.</span>, <span style="color: #666666">1</span>, <span style="color: #666666">0</span> ]
SYS<span style="color: #666666">.</span>extend( [ SPHERE, x<span style="color: #666666">+</span>Rsph<span style="color: #666666">*</span>dFx, y<span style="color: #666666">+</span>Rsph<span style="color: #666666">*</span>dFy, z<span style="color: #666666">+</span>Rsph<span style="color: #666666">*</span>dFz, Rsph ] )
cmd<span style="color: #666666">.</span>load_cgo(SYS, <span style="color: #BB4444">'</span>SYS<span style="color: #BB4444">'</span>, Ifrm<span style="color: #666666">+30</span>)
<span style="color: #AA22FF; font-weight: bold">else</span>:
<span style="color: #008800; font-style: italic"># 初始位置和速度</span>
x <span style="color: #666666">=</span> <span style="color: #666666">-</span>Xini; y<span style="color: #666666">=-</span>Xini
vx <span style="color: #666666">=</span> <span style="color: #666666">3.28</span>; vy <span style="color: #666666">=</span> <span style="color: #666666">0</span>; Nfrm<span style="color: #666666">=130</span>
vx <span style="color: #666666">=</span> <span style="color: #666666">0</span>; vy <span style="color: #666666">=</span> <span style="color: #666666">3.15</span>; Nfrm<span style="color: #666666">=120</span>
dt <span style="color: #666666">=</span> <span style="color: #666666">0.05</span>; E0 <span style="color: #666666">=</span> <span style="color: #666666">0.5*</span>(vx<span style="color: #666666">*</span>vx<span style="color: #666666">+</span>vy<span style="color: #666666">*</span>vy)<span style="color: #666666">+</span>Fxy(x,y)
dFx0, dFy0, dFz0 <span style="color: #666666">=</span> dFxy(x,y)
<span style="color: #AA22FF; font-weight: bold">for</span> Ifrm <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #AA22FF">range</span>(Nfrm):
SYS <span style="color: #666666">=</span> [ COLOR, <span style="color: #666666">1</span>, <span style="color: #666666">1</span>, <span style="color: #666666">1</span> ]
x <span style="color: #666666">=</span> x <span style="color: #666666">+</span> (vx <span style="color: #666666">+</span> <span style="color: #666666">0.5*</span>dFx0<span style="color: #666666">*</span>dt)<span style="color: #666666">*</span>dt;
y <span style="color: #666666">=</span> y <span style="color: #666666">+</span> (vy <span style="color: #666666">+</span> <span style="color: #666666">0.5*</span>dFy0<span style="color: #666666">*</span>dt)<span style="color: #666666">*</span>dt;
z <span style="color: #666666">=</span> Fxy(x,y)
dFx, dFy, dFz <span style="color: #666666">=</span> dFxy(x,y)
SYS<span style="color: #666666">.</span>extend( [ SPHERE, x<span style="color: #666666">+</span>Rsph<span style="color: #666666">*</span>dFx, y<span style="color: #666666">+</span>Rsph<span style="color: #666666">*</span>dFy, z<span style="color: #666666">+</span>Rsph<span style="color: #666666">*</span>dFz, Rsph ] )
cmd<span style="color: #666666">.</span>load_cgo(SYS, <span style="color: #BB4444">'</span>SYS<span style="color: #BB4444">'</span>, Ifrm)
vx <span style="color: #666666">=</span> vx<span style="color: #666666">+0.5*</span>(dFx0<span style="color: #666666">+</span>dFx)<span style="color: #666666">*</span>dt
vy <span style="color: #666666">=</span> vy<span style="color: #666666">+0.5*</span>(dFy0<span style="color: #666666">+</span>dFy)<span style="color: #666666">*</span>dt
E <span style="color: #666666">=</span> <span style="color: #666666">0.5*</span>(vx<span style="color: #666666">*</span>vx<span style="color: #666666">+</span>vy<span style="color: #666666">*</span>vy)
Rtmp <span style="color: #666666">=</span> math<span style="color: #666666">.</span>sqrt(<span style="color: #666666">2.*</span><span style="color: #AA22FF">abs</span>(E0<span style="color: #666666">-</span>z)<span style="color: #666666">/</span>(vx<span style="color: #666666">*</span>vx<span style="color: #666666">+</span>vy<span style="color: #666666">*</span>vy))
vx <span style="color: #666666">=</span> vx<span style="color: #666666">*</span>Rtmp; vy <span style="color: #666666">=</span> vy<span style="color: #666666">*</span>Rtmp
dFx0, dFy0 <span style="color: #666666">=</span> dFx, dFy
cmd<span style="color: #666666">.</span>reset()
cmd<span style="color: #666666">.</span>zoom(<span style="color: #BB4444">'</span>PES<span style="color: #BB4444">'</span>, <span style="color: #666666">1.0</span>)
cmd<span style="color: #666666">.</span>clip(<span style="color: #BB4444">'</span>far<span style="color: #BB4444">'</span>, <span style="color: #666666">-10.0</span>)
cmd<span style="color: #666666">.</span>turn(<span style="color: #BB4444">'</span>z<span style="color: #BB4444">'</span>, <span style="color: #666666">30</span>)
cmd<span style="color: #666666">.</span>turn(<span style="color: #BB4444">'</span>x<span style="color: #BB4444">'</span>, <span style="color: #666666">-60</span>)
cmd<span style="color: #666666">.</span>mplay()
</pre></div>
</td></tr></table>