65.9K
CodeProject 正在变化。 阅读更多。
Home

轨道力学入门

charles922
starIconstarIconstarIconstarIconstarIcon

5.00/5 (13投票s)

2018 年 12 月 6 日

CPOL

11分钟阅读

viewsIcon

20977

downloadIcon

1588

轨道力学入门 - 二体问题

Orbital Mechanics

引言

这个主要项目叫做 OrbitMech。另外两个 Visual Studio 项目在 下方 有详细描述。OrbitMech 是一个使用 C#、Windows Forms 和 OpenGL 编写的 Windows Forms 项目。它允许您通过更改轨道参数来可视化一个小天体绕一个大得多的天体运行的轨道，例如卫星绕地球的轨道。您可以修改的参数如下所述，包括：

  1. 半长轴
  2. 偏心率
  3. 倾角
  4. 升交点经度（Omega）

背景

轨道力学是一个非常复杂的课题，但一个简化的子集被称为“二体问题”，其中一个物体绕一个质量大得多的物体运行，涉及到相当简单的数学。例如，一颗卫星绕地球运行，一艘宇宙飞船绕月球运行，或者一颗行星绕太阳运行。也就是说，我们假设小天体的质量相对于大天体的质量可以忽略不计。行星运动第一定律，由 开普勒 发现，指出轨道的形状是一个椭圆，其中质量较大的天体位于椭圆的一个焦点上。牛顿后来发展了 椭圆轨道 的数学。质量较大的天体将被称为“行星”。

ellipse2

正如您在几何课上所记得的，一个 *椭圆* 是一个平面点集，其中与两个固定点（称为焦点）的距离之和（显示为虚线）是一个常数。椭圆由 *半长轴* 和 *偏心率* 定义。椭圆的 *长轴* 是其最长的直径，是一条穿过中心和两个焦点的线段；半长轴是长轴的一半。偏心率的值从 0 到 1，用于衡量椭圆的“非圆形”程度；圆是偏心率为零的椭圆。在左侧的图像中，椭圆的偏心率约为 0.6。还显示了较短的轴，即 *短轴*。

有六个参数可以完全定义一个轨道，其中四个将在本项目中介绍。前两个轨道参数定义了椭圆的形状，因此也定义了轨道的形状。它们是椭圆的 *半长轴* 和 *偏心率*。

在上面的 图像 中，行星（质量较大的天体，本例中为地球）以及轨道（蓝色）在一个 Windows 窗体中显示。XYZ 轴分别与零度经度、九十度经度和北极对齐。窗体底部有四个滑块。通过使用相应的滑块可以更改半长轴和偏心率的值。请注意，当您更改偏心率时，椭圆的半短轴会显示在半长轴值旁边（以括号表示）。包含轨道的平面，即 *轨道平面*，显示为灰色。长轴显示为青色，短轴显示为品红色。

接下来的两个轨道参数定义了轨道平面的方向。它们是 *倾角*，它是轨道平面相对于升交点的垂直倾斜度（在轨道向上穿过轨道平面的地方测量），以及 *升交点经度*，它在水平上定向椭圆的升交点。升交点经度通常用希腊字母 omega (Ω) 来表示。这两个参数可以通过使用相应的滑块来更改。（为简单起见，我选择不使用最后两个轨道参数：*近拱角* 和 *平近点角*。）轨道周期（分钟）、远地点（轨道的最高点）和远地点的高度显示在窗体的左下角。

使用菜单栏上的“显示”，您可以切换坐标轴、行星和轨道平面的绘制。可以使用左、右、上、下、加号和减号键来移动“摄像机”。

  • 上/下箭头：以 (F) 度倾斜摄像机
  • 右/左箭头：以 (?) 度旋转摄像机
  • 加号/减号：以 R 单位缩放摄像机 内/外

三个摄像机设置显示在窗体的右下角。phi、theta 和 R 的计算包含在 Camera 类中，该类使用 OpenGL 实用方法 gluLookAt()。

OrbitMech 配置文件

4 个轨道参数（偏心率、倾角、半长轴和升交点经度或 Omega），以及摄像机设置，可以在 App.config 文件中指定。例如，这是一个 App.config 的片段，用于指定偏心率为 0.05，倾角为 50 度，半长轴为 6556 公里，Omega 为 30 度。

<setting name="Eccentricity" serializeAs="String">
	<value>0.05</value>
  </setting>
<setting name="Inclination" serializeAs="String">
	<value>50</value>
</setting>
<setting name="SemiMajorAxis" serializeAs="String">
	<value>6556</value>
 </setting>
<setting name="Omega" serializeAs="String">
	<value>30</value>
   </setting>

这是一个 App.config 的片段，用于指定摄像机倾斜（phi, F）为 20 度，旋转（theta, ?）为 60 度，以及缩放为 27。

<setting name="CameraPhi" serializeAs="String">
	<value>20</value>
</setting>
<setting name="CameraTheta" serializeAs="String">
	<value>60</value>
</setting>
<setting name="CameraR" serializeAs="String">
	<value>27</value>
</setting>

背景颜色和轨道平面的颜色也可以在 App.config 文件中指定，使用 ARGB 值或颜色名称，例如：

<setting name="OrbitPlaneColor" serializeAs="String">
	<value>26, 224, 224, 224</value>
</setting>
<setting name="BackgroundColor" serializeAs="String">
	<value>Gray</value>
</setting>

当您退出 OrbitMech 时，系统会询问您是否要将您的轨道参数和摄像机设置保存在应用程序设置中（在 Windows 7 的 AppData\Local 中），并在启动 OrbitMech 时，系统会询问您是否要恢复这些设置。有关详细信息，请参阅 OrbitMechMain.cs 中的 SaveValues() 和 RetrieveValues()。或者，轨道参数和摄像机设置可以通过菜单栏上的“文件->另存为...”保存在文本文件中，并通过“文件->打开...”恢复。

设置行星

计算位置和速度的关键常数是引力常数与行星质量的乘积，即 *标准引力参数*，通常用希腊字母 Mu (μ) 表示。您可以使用 App.config 文件中的配置项 Mu 来更改此值。例如，火星的 μ = 42828。您还可以通过 Textures 目录中存储的位图来更改应用于球体的纹理。Textures 目录及其内容必须复制到 bin\Debug 和 bin\Release 目录中。行星纹理通过 App.config 文件中的配置项 PlanetTexture 设置。例如，将其设置为“mars.jpg”将用火星的位图绘制球体。（感谢 http://vasilydev.blogspot.com 提供将位图绘制到球体上的代码。）您可以使用 App.config 文件中的配置项 PlanetRadiusKM 来设置行星的半径（以公里为单位），天体绕其运行。

这是 App.config 的一个片段，展示了如何将 Earth 设置为 planet。

<setting name="PlanetTexture" serializeAs="String">
	<value>earth.jpg</value>
</setting>
<setting name="PlanetRadiusKM" serializeAs="String">
	<value>6371.1</value>
</setting>
<setting name="Mu" serializeAs="String">
	<value>398600.0</value>
</setting>

这是 App.config 的一个片段，展示了如何将 Mars 设置为 planet。

<setting name="PlanetTexture" serializeAs="String">
	<value>mars.jpg</value>
</setting>
<setting name="PlanetRadiusKM" serializeAs="String">
	<value>3389.1</value>
</setting>
<setting name="Mu" serializeAs="String">
	<value>42828.0</value>
</setting>

这是 App.config 的一个片段，展示了如何将 Moon 设置为“planet”。

<setting name="PlanetTexture" serializeAs="String">
	<value>moon.jpg</value>
</setting>
<setting name="PlanetRadiusKM" serializeAs="String">
	<value>1737</value>
</setting>
<setting name="Mu" serializeAs="String">
	<value>4902.8</value>
</setting>

设置高度

如果您想输入轨道的离地高度，并让 OrbitMech 计算偏心率和半长轴，请将偏心率和半长轴都设置为 -1（负一），并在 App.config 的以下片段中指定离地高度（以公里为单位）：

<setting name="Eccentricity" serializeAs="String">
	<value>-1</value>
</setting>
<setting name="SemiMajorAxis" serializeAs="String">
	<value>-1</value>
</setting>
<setting name="MaxAltitude" serializeAs="String">
	<value>297</value>
</setting>
<setting name="MinAltitude" serializeAs="String">
	<value>111</value>
</setting>

阿波罗 8 号，第一次载人登月任务，其初始轨道倾角为 32 度，最小和最大离地高度分别为 297 和 111 公里。使用上述 App.config 片段中所示的值，OrbitMech 计算出半长轴为 1941 公里，偏心率为 0.048，如下图所示。（请注意，我如何通过“显示->关闭坐标轴”关闭了坐标轴，并在 App.config 中设置了背景和轨道平面颜色。）

Apollo8Orbit

Using the Code

下载并解压项目，然后用 Visual Studio 打开 OrbitMech.sln。它包含三个项目：

  1. OrbitMech - 一个 Windows Forms 项目，使用 OpenGL 以 3D 方式显示轨道。
  2. KeplerForm - 一个 Windows Forms 项目，使用 OpenGL 以 3D 方式模拟一个天体绕行星运行的动画。它使用 KeplerLib 进行计算。
  3. KeplerLib - 一个用于数值求解开普勒方程的库（有关数值方法的介绍，请参见 我的数值方法项目）。

使用 Visual Studio 主菜单中的“生成”->“重新生成解决方案”来生成解决方案。按 Ctrl-F5 启动而不调试，并显示上面图像中所示的 Windows 窗体。

OpenGL

虽然本文不是 OpenGL 教程，但这里为不熟悉 OpenGL 的人提供了一些解释。“OpenGL 编程指南”，也称为“红宝书”，是学习 OpenGL 的绝佳资源。

例如，这是 OrbitalBody.cs 中 DrawPlane() 用于绘制轨道平面的代码。在 OpenGL 中，所有几何对象最终都描述为一对 glBegin() 和 glEnd() 之间的顶点有序集（请注意，我使用了 C# 的花括号，以便您可以轻松地使用 Visual Studio 的 CTRL-] 匹配 glBegin() 和 glEnd() 调用 - 这只是一个便利，OpenGL 本身不需要）。在此示例中，轨道平面是一个带有四个顶点的矩形。矩形的尺寸由椭圆的半长轴和半短轴定义。矩形的角由四个 glVertex3f() 调用定义，并且矩形位于 XZ 平面上，因为 glVertex3f() 调用中的 Y 坐标（第二个参数）为零。

Gl.glBegin(Gl.GL_QUADS);  // Each set of 4 vertices form a quad
{
...
Gl.glVertex3f(-semiMajorAxis + focus, 0.0f, -semiMinorAxis);
Gl.glVertex3f(+semiMajorAxis + focus, 0.0f, -semiMinorAxis);
Gl.glVertex3f(+semiMajorAxis + focus, 0.0f, +semiMinorAxis);
Gl.glVertex3f(-semiMajorAxis + focus, 0.0f, +semiMinorAxis);
}
Gl.glEnd();

OpenGL 的一个关键概念是模型视图矩阵，用于应用 *模型变换*，例如旋转或平移场景的一部分。例如，这是 OrbitalBody.cs 中 DrawEllipse() 用于绘制实际轨道中的代码。当通过移动倾角滑块更改轨道平面的倾角时，从滑块控件获得的倾角值（以度为单位）被用作 glRotatef()（一个 *模型变换* 方法）的第一个参数。接下来的三个参数指定了旋转发生的向量的 x、y 和 z 分量，在本例中为 X 轴。通过 omega 进行旋转的代码是类似的 glRotatef() 调用。

Gl.glPushMatrix();
{
Gl.glRotatef(orbitalParameters.Inclination, 1.0f, 0.0f, 0.0f);
		...
DrawPlane(semiMajorAxis, semiMinorAxis, (float)focus);
	...
}
Gl.glPopMatrix();

正如红宝书所说：

实际上，glPushMatrix() 意味着“记住你现在的位置”，而 glPopMatrix() 意味着“回到你原来的位置”。

- OpenGL 红宝书

因此，上面的代码的效果是旋转轨道平面——即改变它的倾角——而不旋转场景中的任何其他东西。同样，我使用了 C# 的花括号，以便您可以使用 CTRL-] 匹配 glPushMatrix() 和 glPopMatrix() 调用，这只是一个便利，OpenGL 本身不需要。

Tao Framework

OrbitMech 使用 Tao Framework OpenGL 实现。这些文件包含在 OrbitMech 项目中，并且必须复制到 bin\Debug 和 bin\Release 目录中。

  1. Tao.OpenGL.dll
  2. ShadowEngine.dll
  3. glut32.dll

开普勒方程的数值解

Kepler

为了计算天体运行时位置和速度，有必要数值求解开普勒方程。
我感谢 Richard Gonsalves 他实现了开普勒方程的龙格-库塔数值积分，我将其从 C++ 转换为 C# 并放入一个名为 KeplerLib 的库中。我修改了它以迭代一个轨道，并以公里为单位显示距离，以公里/秒为单位显示速度。在 Visual Studio 解决方案资源管理器中右键单击“Kepler”，然后单击“设置为启动项目”。按 Ctrl-F5 启动而不调试，并显示上面图像中所示的 Windows 窗体。卫星绕地球的轨道具有以下轨道参数：

  • 偏心率 = 0.6
  • 半长轴 = 18125 公里
  • 倾角 = 0
  • Omega = 0

（在此图像中，我使用向下箭头键将摄像机的 Phi (F) 值设置为 65 度。摄像机的 Phi、Theta 和 R 值显示在控制台窗口中。）

单击“开始”按钮将启动一个演示卫星绕行星运行的动画。（按下“开始”后，按钮标签将变为“暂停”。我按下了“暂停”以获取屏幕截图。按下“继续”将继续动画。）构造函数 KeplerForm 调用库方法 KeplerCalc.Integrate()，该方法反过来求解运动方程并创建一个 *状态向量* 列表，List<StateVector> svList。每个状态向量代表轨道天体在某个瞬时时刻的位置和速度。这个状态向量列表及其卫星的位置和速度由 KeplerForm 的 DrawOrbit() 方法使用。当它迭代 svList 中的每个 StateVector 时，DrawOrbit() 会在 StateVector 的位置（一个 PointF 结构）处绘制卫星（一个小的红色方块）。一个名为 tmrPaint 的计时器用于驱动动画。计时器的间隔根据卫星的速度计算，因此在动画中，卫星将在其最接近地球（近地点）的点附近加速，并在到达其最远点（远地点）时减速，这是开普勒本人观察到的现象。

控制台窗口

与我的 NXT 蓝牙监视器项目 一样，我为 Windows Forms 项目附加了一个控制台用于调试目的，使用了 AllocConsole 和 .NET Interop。它默认启用；要禁用它，只需在 OrbitMechMain.cs 或 KeplerFormMain.cs 的 Main() 的前 7 行中注释掉即可。

IntPtr ptr = GetForegroundWindow();
int pid;
GetWindowThreadProcessId(ptr, out pid);
Process process = Process.GetProcessById(pid);
AllocConsole();
AttachConsole(process.Id);
Console.WriteLine("{0} OrbitMech started", DateTime.Now);

历史

  • 版本 2.0.0.0 Omega 现在显示为东经或西经（单击帮助以查看版本）。
© . All rights reserved.