需求
搜集到的方法和问题
- 直接用牛顿力学公式进行模拟的问题
- - 物理常数调整问题
  - 初始化参量不直观
- 实时求解动力学微分方程
插值方法
- 问题简化
- 求解 θ ( t ) \theta(t) θ(t) (?)
- - 常数项求解
  - - Note：开普勒方程
- 回到直角坐标
- 插值方案选择
- 行星椭圆轨道的三次样条插值
- - 对 θ \theta θ 进行插值的方法
  - 对 r ⃗ \vec r r 进行插值的方案
  - 插值取样方法
  - - 1. τ \tau τ 为负数的问题
    - 2. 取样点分布的问题

需求

给定椭圆轨道参数（轨道倾角、升交点经度、轨道周期、轨道半长轴、轨道偏心率、初始相位或初始位置），在游戏中实时演算椭圆轨道中行星的位置。

可以忽略轨道进动
节约计算性能
在大轨道半长轴时精度损失相对小

笔者在网络上找了很久，都没有找到一篇完整的适合游戏中计算任意时刻行星在椭圆轨道上位置的文章。于是决定自己研究一下。

搜集到的方法和问题

直接用牛顿力学公式进行模拟的问题

这是最直觉的方法，但问题在于给出合适的初始参数。

物理常数调整问题

在游戏中一般不会直接使用现实中各个物理量的值，例如恒星质量、行星质量、万有引力常数 G 等，游戏开发者往往需要根据需求对现实世界进行缩放。

初始化参量不直观

另外，利用牛顿力学公式，需要根据给定的轨道计算出行星的初始位置和初始速度。一旦参数调整出现一点点小问题，轨道便会偏离设定（例如巨大的轨道进动），更有甚者甚至无法形成稳定的椭圆轨道。
而这和常用描述轨道的参数例如 轨道倾角、升交点经度、轨道周期、轨道半长轴、轨道偏心率 没有直观的关系——我们并不能从初始位置和初始速度简单地推测出这些参数。

实时求解动力学微分方程

实时解微分方程是什么性能黑洞，更别提精度浪费的问题了。

插值方法

问题简化

我们先不考虑轨道倾角和升交点经度，把行星的轨道放到平面上。

根据高中物理书上的知识：行星的轨道呈椭圆形，其围绕的恒星在椭圆的一个焦点上。
那么，把所围绕的恒星放在坐标原点，用极坐标表示行星的椭圆轨道为：
r = p 1 + e cos ⁡ θ \begin{equation} r=\frac{p}{1+e \cos\theta} \end{equation} r=1+ecosθp
这里 p p p 是椭圆的通径、 e e e 是椭圆的离心率， p = b 2 / a , 0 < e < 1 p=b^2/a, 0<e<1 p=b2/a,0<e<1.
p = a ( 1 − e 2 ) p=a(1-e^2) p=a(1−e2)
这样一来，行星在任意时刻的位置问题就变成了求解 θ ( t ) {\theta(t)} θ(t) 的问题。

求解 θ ( t ) \theta(t) θ(t) (?)

这里复习一下大学物理知识：角动量的定义和角动量守恒（开普勒第二定律）
J : = r 2 θ ˙ = c o n s t . \begin{equation} J:=r^2 \.\theta=\mathrm{const.} \end{equation} J:=r2θ˙=const.
这里 θ ˙ = d θ d t \.\theta=\frac{\mathrm{d}\theta}{\mathrm{d}t} θ˙=dtdθ，转化为 θ ( t ) \theta(t) θ(t) 的方程，并设无量纲化时间 τ = t T \tau=\frac{t}{T} τ=Tt 使方程核心部分无量纲化，得到
p 2 J T θ ′ ( τ ) ( 1 + e cos ⁡ θ ) 2 = 1 \frac{p^2}{JT}\frac{\theta'(\tau)}{(1+e \cos\theta)^2}=1 JTp2(1+ecosθ)2θ′(τ)=1
显然这个微分方程中很难解出 θ ( τ ) \theta(\tau) θ(τ)，但没关系，我们的目的是获取 θ ( τ ) \theta(\tau) θ(τ) 的插值，而且 θ ( τ ) \theta(\tau) θ(τ) 显然是个单调递增的函数，因此也可以去求解他的反函数 τ ( θ ) \tau(\theta) τ(θ)，有
τ ′ ( θ ) = p 2 J T 1 ( 1 + e cos ⁡ θ ) 2 \begin{equation} \tau'(\theta)=\frac{p^2}{JT}\frac{1}{(1+e \cos\theta)^2} \end{equation} τ′(θ)=JTp2(1+ecosθ)21
i.e.
τ ( θ ) = p 2 J T ∫ 0 θ 1 ( 1 + e cos ⁡ ϕ ) 2 d ϕ \begin{equation} \tau(\theta)=\frac{p^2}{JT}\int_0^\theta \frac{1}{(1+e \cos\phi)^2}\mathrm{d}\phi \end{equation} τ(θ)=JTp2∫0θ(1+ecosϕ)21dϕ
由此，我们可以用积分的方法获得 θ ( τ ) \theta(\tau) θ(τ) 上的坐标点。现在只剩下一个问题：如何得到式子中的常数项。（虽然有物理直觉的同学可以看出：这个常数项只与偏心率 e e e 有关，而且大概率含一个 2 π 2\pi 2π）

常数项求解

因为 τ = t / T \tau=t/T τ=t/T 是无量纲化的时间，所以行星运行一周，即 θ = 2 π \theta=2\pi θ=2π 时， τ = 1 \tau=1 τ=1
这个活就交给 Mathematica 了，答案是 2 π ( 1 − e 2 ) 3 / 2 \frac{2 \pi }{\left(1-e^2\right)^{3/2}} (1−e2)3/22π，所以
τ ( θ ) = τ 0 + ( 1 − e 2 ) 3 / 2 2 π ∫ 0 θ 1 ( 1 + e cos ⁡ ϕ ) 2 d ϕ \begin{equation} \tau(\theta)=\tau_0+\frac{\left(1-e^2\right)^{3/2}}{2 \pi }\int_0^\theta \frac{1}{(1+e \cos\phi)^2}\mathrm{d}\phi \end{equation} τ(θ)=τ0+2π(1−e2)3/2∫0θ(1+ecosϕ)21dϕ
这里 τ 0 \tau_0 τ0 是行星的初始“相位”。从这个公式就可以直接计算 τ ( θ ) \tau(\theta) τ(θ) 曲线上坐标点了。

Note：开普勒方程

推完一遍发现这个积分有解析解，Mathematica 解出来是

τ ( θ ) = ( 1 − e 2 ) 3 / 2 2 π ( e sin ⁡ ( θ ) ( e 2 − 1 ) ( e cos ⁡ ( θ ) + 1 ) + 2 tanh ⁡ − 1 ( ( e − 1 ) tan ⁡ ( θ 2 ) e 2 − 1 ) ( e 2 − 1 ) 3 / 2 ) \begin{aligned} &\tau(\theta)\\ =&\frac{\left(1-e^2\right)^{3/2}}{2 \pi }\left(\frac{e \sin (\theta )}{\left(e^2-1\right) (e \cos (\theta )+1)}+\frac{2 \tanh ^{-1}\left(\frac{(e-1) \tan \left(\frac{\theta }{2}\right)}{\sqrt{e^2-1}}\right)}{\left(e^2-1\right)^{3/2}}\right) \end{aligned} =τ(θ)2π(1−e2)3/2 (e2−1)(ecos(θ)+1)esin(θ)+(e2−1)3/22tanh−1(e2−1 (e−1)tan(2θ))
可以化简为
τ ( θ ) = 1 2 π ( E − e sin ⁡ E ) , [ E = 2 arctan ⁡ 1 − e 1 + e tan ⁡ ( θ 2 ) ] \begin{equation} \tau(\theta)=\frac1{2\pi}(E-e\sin E),\quad \left[E=2\arctan \sqrt{\frac{1-e}{1+e}}\tan\left(\frac{\theta}{2}\right)\right] \end{equation} τ(θ)=2π1(E−esinE),[E=2arctan1+e1−e tan(2θ)]
然后这个东西叫做开普勒方程。

_{所以我是重新推了一遍开普勒方程。-_-||}

回到直角坐标

现在我们获得了计算 θ ( t ) \theta(t) θ(t) 上坐标点的方法，进行插值后，也就获取了计算平面内行星的实时位置的方法。
设轨道倾角为 α \alpha α，升交点经度为 Ω \Omega Ω，转换为旋转矩阵，得到

( x y z ) = R z ( Ω ) R x ( α ) ( ρ [ θ ( t ) ] cos ⁡ θ ( t ) ρ [ θ ( t ) ] sin ⁡ θ ( t ) 0 ) \begin{equation} \begin{pmatrix} x\\ y\\ z \end{pmatrix}=R_z(\Omega)R_x(\alpha)\begin{pmatrix} \rho[\theta(t)]\cos \theta(t)\\ \rho[\theta(t)]\sin \theta(t)\\ 0 \end{pmatrix} \end{equation} xyz =Rz(Ω)Rx(α) ρ[θ(t)]cosθ(t)ρ[θ(t)]sinθ(t)0
( x y z ) = ( cos ⁡ Ω − sin ⁡ Ω 0 sin ⁡ Ω cos ⁡ Ω 0 0 0 1 ) ( 1 0 0 0 cos ⁡ α − sin ⁡ α 0 sin ⁡ α cos ⁡ α ) ( ρ [ θ ( t ) ] cos ⁡ θ ( t ) ρ [ θ ( t ) ] sin ⁡ θ ( t ) 0 ) \begin{pmatrix} x\\ y\\ z \end{pmatrix}= \begin{pmatrix} \cos\Omega&-\sin\Omega&0\\ \sin\Omega&\cos\Omega&0\\ 0&0&1 \end{pmatrix} \begin{pmatrix} 1&0&0\\ 0&\cos\alpha&-\sin\alpha\\ 0&\sin\alpha&\cos\alpha \end{pmatrix} \begin{pmatrix} \rho[\theta(t)]\cos \theta(t)\\ \rho[\theta(t)]\sin \theta(t)\\ 0 \end{pmatrix} xyz = cosΩsinΩ0−sinΩcosΩ0001 1000cosαsinα0−sinαcosα ρ[θ(t)]cosθ(t)ρ[θ(t)]sinθ(t)0

这里往下有两套计算方案：

将 θ ( t ) \theta(t) θ(t) 对应的空间直角坐标点 ( x , y , z ) (x,y,z) (x,y,z) 算出，并分别对 x ( t ) , y ( t ) , z ( t ) x(t), y(t), z(t) x(t),y(t),z(t) 进行插值。
直接对 θ ( t ) \theta(t) θ(t) 进行插值，实时演算时求出对应的 θ \theta θ，再用坐标变换计算 ( x , y , z ) (x,y,z) (x,y,z)。

插值方案选择

线性插值
- 取样点多
- 计算简单
牛顿法 n 次样条
- 承袭性，加点方便
三次样条插值（全周期，周期边界条件）
- 取样点少
- 计算稍微复杂

_{不会有人用拉格朗日插值来做这个吧~~不会吧不会吧~~}

行星椭圆轨道的三次样条插值

对 θ \theta θ 进行插值的方法

三次样条插值的要求是在插值点处，函数的二阶导连续的三次多项式分段插值

那么，根据方程 (6) 和微分方程 (3)，我们可以给出每个取样点（给定 θ i \theta_i θi）的值以及对应的一阶导数
τ i = 1 2 π ( E i − e sin ⁡ E i ) , [ E i = 2 arctan ⁡ 1 − e 1 + e tan ⁡ ( θ i 2 ) ] θ i ′ = 2 π ( 1 − e 2 ) 3 / 2 ( 1 + e cos ⁡ θ i ) 2 \tau_i=\frac1{2\pi}(E_i-e\sin E_i),\quad \left[E_i=2\arctan \sqrt{\frac{1-e}{1+e}}\tan\left(\frac{\theta_i}{2}\right)\right]\\ \theta_i'=\frac{2 \pi }{\left(1-e^2\right)^{3/2}}{(1+e \cos\theta_i)^2} τi=2π1(Ei−esinEi),[Ei=2arctan1+e1−e tan(2θi)]θi′=(1−e2)3/22π(1+ecosθi)2

由此可以计算得到每一段上三次样条的函数，假设写成如下形式
θ i ( τ ) = ( τ − τ i ) ( ( τ − τ i ) ( a i ( τ − τ i ) + b i ) + c i ) + d i ( τ i ≤ τ < τ i + 1 ) \theta_i(\tau)=(\tau-\tau_i)((\tau-\tau_i)(a^i(\tau-\tau_i)+b^i)+c^i)+d^i \quad(\tau_i\le\tau<\tau_{i+1}) θi(τ)=(τ−τi)((τ−τi)(ai(τ−τi)+bi)+ci)+di(τi≤τ<τi+1)
那么有
{ a i = ( θ i ′ + θ i + 1 ′ ) h i − 2 ( θ i + 1 − θ i ) h i 3 b i = − h i ( 2 θ i ′ + θ i + 1 ′ ) − 3 ( θ i + 1 − θ i ) h i 2 c i = θ i ′ d i = θ i \left\{ \begin{aligned} a^i&=\frac{(\theta'_i+\theta'_{i+1})h_i-2(\theta_{i+1}-\theta_i)}{h_i^3}\\ b^i&=-\frac{h_i(2\theta'_i+\theta'_{i+1})-3(\theta_{i+1}-\theta_i)}{h_i^2}\\ c^i&=\theta'_i\\ d^i&=\theta_i\\ \end{aligned} \right. ⎩ ⎨ ⎧aibicidi=hi3(θi′+θi+1′)hi−2(θi+1−θi)=−hi2hi(2θi′+θi+1′)−3(θi+1−θi)=θi′=θi

对 r ⃗ \vec r r 进行插值的方案

既然已经有了 θ ′ \theta' θ′，那么也可以计算得到 x , y , z x, y, z x,y,z 的插值函数。星体最终的坐标是从轨道平面上的坐标线性变换得来的
( x ( t ) y ( t ) z ( t ) ) = ( cos ⁡ Ω − sin ⁡ Ω 0 sin ⁡ Ω cos ⁡ Ω 0 0 0 1 ) ( 1 0 0 0 cos ⁡ α − sin ⁡ α 0 sin ⁡ α cos ⁡ α ) ( ρ [ θ ( t ) ] cos ⁡ θ ( t ) ρ [ θ ( t ) ] sin ⁡ θ ( t ) 0 ) \begin{pmatrix} x(t)\\ y(t)\\ z(t) \end{pmatrix}= \begin{pmatrix} \cos\Omega&-\sin\Omega&0\\ \sin\Omega&\cos\Omega&0\\ 0&0&1 \end{pmatrix} \begin{pmatrix} 1&0&0\\ 0&\cos\alpha&-\sin\alpha\\ 0&\sin\alpha&\cos\alpha \end{pmatrix} \begin{pmatrix} \rho[\theta(t)]\cos \theta(t)\\ \rho[\theta(t)]\sin \theta(t)\\ 0 \end{pmatrix} x(t)y(t)z(t) = cosΩsinΩ0−sinΩcosΩ0001 1000cosαsinα0−sinαcosα ρ[θ(t)]cosθ(t)ρ[θ(t)]sinθ(t)0
所以有
( x ′ ( t ) y ′ ( t ) z ′ ( t ) ) = ( cos ⁡ Ω − sin ⁡ Ω 0 sin ⁡ Ω cos ⁡ Ω 0 0 0 1 ) ( 1 0 0 0 cos ⁡ α − sin ⁡ α 0 sin ⁡ α cos ⁡ α ) ( [ ρ [ θ ( t ) ] cos ⁡ θ ( t ) ] ′ [ ρ [ θ ( t ) ] sin ⁡ θ ( t ) ] ′ 0 ) \begin{pmatrix} x'(t)\\ y'(t)\\ z'(t) \end{pmatrix}= \begin{pmatrix} \cos\Omega&-\sin\Omega&0\\ \sin\Omega&\cos\Omega&0\\ 0&0&1 \end{pmatrix} \begin{pmatrix} 1&0&0\\ 0&\cos\alpha&-\sin\alpha\\ 0&\sin\alpha&\cos\alpha \end{pmatrix} \begin{pmatrix} [\rho[\theta(t)]\cos \theta(t)]'\\ [\rho[\theta(t)]\sin \theta(t)]'\\ 0 \end{pmatrix} x′(t)y′(t)z′(t) = cosΩsinΩ0−sinΩcosΩ0001 1000cosαsinα0−sinαcosα [ρ[θ(t)]cosθ(t)]′[ρ[θ(t)]sinθ(t)]′0
而
[ ρ [ θ ( t ) ] cos ⁡ θ ( t ) ] ′ = − sin ⁡ θ ( t ) ⋅ f [ ρ [ θ ( t ) ] sin ⁡ θ ( t ) ] ′ = ( e + cos ⁡ θ ( t ) ) ⋅ f ( f = 1 p ρ 2 ( t ) θ ′ ( t ) = J / p = 2 π ( 1 − e 2 ) 3 / 2 p T = const . ) [\rho[\theta(t)]\cos \theta(t)]'=-\sin \theta(t)\cdot f\\ [\rho[\theta(t)]\sin \theta(t)]'=(e+\cos\theta(t)) \cdot f\\ \left(f = \frac1p\rho^2(t)\theta'(t)=J/p=\frac{2 \pi }{\left(1-e^2\right)^{3/2}}\frac pT=\text{const}.\right) [ρ[θ(t)]cosθ(t)]′=−sinθ(t)⋅f[ρ[θ(t)]sinθ(t)]′=(e+cosθ(t))⋅f(f=p1ρ2(t)θ′(t)=J/p=(1−e2)3/22πTp=const.)
至此，便得到了行星轨道直角坐标关于 t t t 的一阶导，代入上面插值系数的公式即可。

插值取样方法

上面的取样公式在具体实现时有两个问题：

θ \theta θ 取到 ( π , 2 π ) (\pi,2\pi) (π,2π) 区间时，按照公式计算 τ \tau τ 会得到负数
在离心率很大时， τ ( θ ) \tau(\theta) τ(θ) 很接近于阶跃函数，因此需要根据离心率来动态地获取采样点。

1. τ \tau τ 为负数的问题

直接计算三角函数得到的 τ ( θ ) \tau(\theta) τ(θ) 的函数图像

产生该现象的原因：在计算 E = 2 arctan ⁡ ( 1 − e 1 + e tan ⁡ θ 2 ) E=2\arctan\left(\sqrt{\frac{1-e}{1+e}}\tan\frac\theta2\right) E=2arctan(1+e1−e tan2θ) 时，默认的 arctan ⁡ \arctan arctan 函数的值域是 ( − π , π ) (-\pi,\pi) (−π,π) 区间，对应到无量纲化时间 τ \tau τ 上是 ( − 1 2 , 1 2 ) \left(-\frac12,\frac12\right) (−21,21)

两种解决方案

修改插值的区间为 τ ∈ [ − 1 2 , 1 2 ] \tau\in\left[-\frac12,\frac12\right] τ∈[−21,21], 对应采样的区间为 θ ∈ [ − π , π ] \theta\in[-\pi,\pi] θ∈[−π,π]
采样区间依然为 θ ∈ [ 0 , 2 π ] \theta\in[0,2\pi] θ∈[0,2π]，但在 ( π , 2 π ] (\pi,2\pi] (π,2π] 上对计算得到的 τ \tau τ 加 1，得到的 τ ( θ ) \tau(\theta) τ(θ) 曲线如图

2. 取样点分布的问题

当轨道离心率比较大时（例如哈雷彗星 0.967），如果仅仅均匀地取多个 θ \theta θ 值，那么在近日点附近的时间分辨率会非常小，轨道将严重失真。为了避免这种情况，可以通过加取样点的方法，保证取样点在时间上的间隔不至于太大。

具体做法：

均匀取 θ \theta θ，计算出对应的 τ \tau τ 值
按 τ \tau τ 排序后，检查两个相邻取样点之间的 Δ τ i = τ i + 1 − τ i \Delta\tau_i=\tau_{i+1}-\tau_i Δτi=τi+1−τi，如果大于我们设定的阈值，那就在这两个取样点的中点处（ θ = 0.5 ( θ i + θ i + 1 \theta=0.5(\theta_i+\theta_{i+1} θ=0.5(θi+θi+1) 取样并加入取样点列表
反复执行第 2 步，直到所有相邻两点间 τ \tau τ 的差值都小于我们设定的阈值

计算任意时刻椭圆轨道上行星的位置的插值方法相关推荐

计算任意时刻格林尼治视恒星时角
近期学习卫星轨道方面的一些知识,遇到计算任意时刻格林尼治视恒星时角的问题,在网上搜了好久也没有一个完整的解决方案,后来通过,网上的一些零碎的信息,终于完成了计算格林尼治视恒星时角的程序,先整理如下. ...
Hologres 弹性计算在 OLAP 分析上的实践和探索
简介: 1.本文介绍了OLAP分析在大数据分析中的位置 2.分析并介绍目前大数据OLAP遇到的分析性能.资源隔离.高可用.弹性扩缩容等核心问题 3.解析阿里云Hologres是如何解决极致性能.弹性. ...
设置android启动器,教程：在任意 Android 设备上安装 HTC 专属桌面启动器
HTC 手机的标志性 Sense UI 几乎是伴随着 Android 系统成长起来的,BlinkFeed 作为 Sense 5 的新特性进入了人们的视野之中.经历了几代更新,BlinkFeed 受到更 ...
T-SQL 2 Tips: 1.计算任意两日期之间的周一到周日分别各有几个! 2.根据出生日期计算精确年龄!...
这两个小技巧,不写不知道,一写吓一跳! 都是看似简单,实际做起来就懵,得仔细想一想,才能写对! 凡是有日期运算的程序都要细心哦! 先说第二个: 2.根据出生日期精确计算年龄! 所谓计算精确年龄就是 ...
android分辨率 x y,Android - 确定不同分辨率的位图上的特定位置（X，Y坐标）？
我正在尝试创建的应用程序是棋盘游戏.它将有一个位图作为棋盘和棋子,移动到棋盘上的不同位置.该板的一般设计是方形的,具有一定数量的行和列,并具有用于外观的边框.想象一下国际象棋棋盘或拼字游戏板.Andr ...
NSURLSessionDataTask与NSOperationQueue实现多文件断点下载（任意时刻终止进程，重启应用，自动重启下载）...
效果展示 gif有点大,直接连接:7qnbrb.com1.z0.glb.clouddn.com/download.gi- 知识要点 NSOperationQueue线程队列的管理 NSURLSessi ...
凸多边形面积_C++计算任意多边形的面积
任意多边形的面积计算_拾忆楓灵的博客-CSDN博客blog.csdn.net 计算任意多边形的面积 - tenos - 博客园www.cnblogs.com 完美解决计算3D空间任意多边形面积_S ...
DSSM、CNN-DSSM、LSTM-DSSM等深度学习模型在计算语义相似度上的应用+距离运算
在NLP领域,语义相似度的计算一直是个难题:搜索场景下query和Doc的语义相似度.feeds场景下Doc和Doc的语义相似度.机器翻译场景下A句子和B句子的语义相似度等等.本文通过介绍DSSM.C ...
【perl计算基因在基因组上的距离】
perl计算基因在基因组上的距离介绍根据输入注释文件及目的基因文件,找出目的基因在基因组上的距离并输出到文件内一.文件格式基因注释文件:为细菌全基因组利用prokka注释文件,格式如下 gen ...

计算任意时刻椭圆轨道上行星的位置的插值方法

目录

需求