jplephem 2.17

安装

jplephem依赖的唯一第三方包是NumPy，当您运行时， pip将自动尝试与pyephem一起安装：

$ pip install jplephem

如果您看到 NumPy 编译错误，请尝试直接从其网站下载和安装 NumPy，或者直接使用已安装科学工具的 Python 发行版，例如 Anaconda。

请注意，jplephem仅提供生成纯 3 维向量所需的逻辑。大多数对天文学感兴趣的程序员会想看看Skyfield ，它使用jplephem但将数字转换为更传统的测量值，如赤经和赤纬。

大多数用户将jplephem与 NASA JPL 的 NAIF 设施提供的卫星行星内核 (SPK) 文件一起使用，以便与他们自己的 SPICE 工具包一起使用。他们在该目录下收集了最有用的内核：

http://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/

要了解有关 SPK 文件的更多信息， 可从 NAIF 设施网站上的 NASA JPL 域下获得官方SPK必读文档。

命令行工具

如果你已经下载了一个.bsp文件，你可以从命令行运行jplephem来显示其中的数据：

python -m jplephem comment de421.bsp
python -m jplephem daf de421.bsp
python -m jplephem spk de421.bsp

您还可以通过限制其涵盖的日期范围来获取较大的星历并生成较小的摘录：

python -m jplephem excerpt 2018/1/1 2018/4/1 de421.bsp excerpt421.bsp

如果您的起始年份为负数，您将收到错误消息，因为 Unix 命令在看到破折号时需要一个选项列表。解决方法是提供一个特殊的参数——它说“我已经完成了传递选项，即使下一个参数带有破折号”：

python -m jplephem excerpt -- -800/1/1 800/1/1 de422.bsp excerpt422.bsp

您还可以按所需目标的整数代码进行过滤。无法识别的目标不会引发错误，让您可以将目标主列表应用于可能具有或不具有所有目标的整个系列的 SPK 文件：

python -m jplephem excerpt --targets 1,2,3 2018/1/1 2018/4/1 de421.bsp excerpt421.bsp

如果输入的星历是一个 URL，那么jplephem将尝试通过仅获取覆盖您指定的日期所需的远程文件块来节省带宽。例如，木星卫星星历jup310.bsp以大而闻名，重达近 1 GB。但是，如果您只需要几个月的木星卫星，您可以下载的数据要少得多：

$ python -m jplephem excerpt 2018/1/1 2018/4/1 \
    https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/satellites/jup310.bsp \
    excerpt.bsp
$ ls -lh excerpt.bsp
-rw-r----- 1 brandon brandon 1.2M Feb 11 13:36 excerpt.bsp

在这种情况下，只需下载大约千分之一的星历数据。

DE421 入门

DE421 星历表是一个有用的起点。它的大小为 17 MB，但提供了 1900-2050 年间的预测：

https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/a_old_versions/de421.bsp

内核下载后，您可以使用jplephem加载此 SPK 文件并了解它提供的段：

>>> from jplephem.spk import SPK
>>> kernel = SPK.open('de421.bsp')
>>> print(kernel)
File type DAF/SPK and format LTL-IEEE with 15 segments:
2414864.50..2471184.50  Type 2  Solar System Barycenter (0) -> Mercury Barycenter (1)
2414864.50..2471184.50  Type 2  Solar System Barycenter (0) -> Venus Barycenter (2)
2414864.50..2471184.50  Type 2  Solar System Barycenter (0) -> Earth Barycenter (3)
2414864.50..2471184.50  Type 2  Solar System Barycenter (0) -> Mars Barycenter (4)
2414864.50..2471184.50  Type 2  Solar System Barycenter (0) -> Jupiter Barycenter (5)
2414864.50..2471184.50  Type 2  Solar System Barycenter (0) -> Saturn Barycenter (6)
2414864.50..2471184.50  Type 2  Solar System Barycenter (0) -> Uranus Barycenter (7)
2414864.50..2471184.50  Type 2  Solar System Barycenter (0) -> Neptune Barycenter (8)
2414864.50..2471184.50  Type 2  Solar System Barycenter (0) -> Pluto Barycenter (9)
2414864.50..2471184.50  Type 2  Solar System Barycenter (0) -> Sun (10)
2414864.50..2471184.50  Type 2  Earth Barycenter (3) -> Moon (301)
2414864.50..2471184.50  Type 2  Earth Barycenter (3) -> Earth (399)
2414864.50..2471184.50  Type 2  Mercury Barycenter (1) -> Mercury (199)
2414864.50..2471184.50  Type 2  Venus Barycenter (2) -> Venus (299)
2414864.50..2471184.50  Type 2  Mars Barycenter (4) -> Mars (499)

由于接下来的几个示例涉及向量输出，让我们告诉 NumPy 使向量输出具有吸引力。

>>> import numpy as np
>>> np.set_printoptions(precision=3)

文件的每个部分都可以让您预测对象相对于其他参考点的位置。如果你想要火星在 2457061.5（2015 年 2 月 8 日）相对于太阳系中心的坐标，这个星历只需要你一步：

>>> position = kernel[0,4].compute(2457061.5)
>>> print(position)
[2.057e+08 4.251e+07 1.394e+07]

但是了解火星相对于地球的位置需要三个步骤，从火星到太阳系质心到地月质心，最后到地球本身：

>>> position = kernel[0,4].compute(2457061.5)
>>> position -= kernel[0,3].compute(2457061.5)
>>> position -= kernel[3,399].compute(2457061.5)
>>> print(position)
[ 3.161e+08 -4.679e+07 -2.476e+07]

你可以看到这个星历的输出是以公里为单位的。如果您使用其他星历，请检查其文档以确保其使用的单位。

如果您将日期作为 NumPy 数组提供，则返回的每个组件本身就是一个向量，只要您的日期：

>>> jd = np.array([2457061.5, 2457062.5, 2457063.5, 2457064.5])
>>> position = kernel[0,4].compute(jd)
>>> print(position)
[[2.057e+08 2.053e+08 2.049e+08 2.045e+08]
 [4.251e+07 4.453e+07 4.654e+07 4.855e+07]
 [1.394e+07 1.487e+07 1.581e+07 1.674e+07]]

一些星历表通过返回 6 向量而不是 3 向量来包含内联速度。对于没有的星历，您可以要求对 Chebyshev 多项式进行微分以产生速度，该速度作为第二个返回值传递：

>>> position, velocity = kernel[0,4].compute_and_differentiate(2457061.5)
>>> print(position)
[2.057e+08 4.251e+07 1.394e+07]
>>> print(velocity)
[-363896.059 2019662.996  936169.773]

默认情况下，速度将是每天行进的距离，以星历表恰好使用的距离单位。要获得每秒的速度，只需除以一天中的秒数：

>>> velocity_per_second = velocity / 86400.0
>>> print(velocity_per_second)
[-4.212 23.376 10.835]

API 的详细信息

这里有一些细节可供准备超越上面提供的高级 API 并通读代码以了解更多信息的人们。

• jplephem不是将整个星历读入内存，而是 对底层文件进行内存映射，以便操作系统可以仅将代码正在使用的数据有效地分页到 RAM。

• 一旦元数据从二进制 SPK 文件中解析出来，多项式系数本身就会通过构建一个可以访问原始二进制文件内容的 NumPy 数组对象来加载。令人高兴的是，NumPy 已经知道如何解释双精度浮点数的压缩数组。您可以通过模块jplephem.daf中的DAF类了解底层 DAF “双精度数组文件”格式，以防您需要在 Python 中打开其他此类数组文件。

• SPK 文件由段组成。当你第一次创建一个SPK 内核对象k时，它会检查文件并创建一个 Segment对象列表，它保存在一个名为 k.segments的属性下的列表中，你可以通过循环在你自己的代码中自由地检查它。

• 除了打印出 SPK 文件时获得的单行摘要之外，还有关于每个段的更多信息，您可以通过要求段详细打印自身来查看：

  >>> segment = kernel[3,399]
  >>> print(segment.describe())
  2414864.50..2471184.50  Type 2  Earth Barycenter (3) -> Earth (399)
    frame=1 source=DE-0421LE-0421
  

• 从内核加载的每个Segment都有许多从 SPK 文件加载的属性：

  >>> from jplephem.spk import BaseSegment
  >>> help(BaseSegment)
  Help on class BaseSegment in module jplephem.spk:
  ...
   |  segment.source - official ephemeris name, like 'DE-0430LE-0430'
   |  segment.start_second - initial epoch, as seconds from J2000
   |  segment.end_second - final epoch, as seconds from J2000
   |  segment.start_jd - start_second, converted to a Julian Date
   |  segment.end_jd - end_second, converted to a Julian Date
   |  segment.center - integer center identifier
   |  segment.target - integer target identifier
   |  segment.frame - integer frame identifier
   |  segment.data_type - integer data type identifier
   |  segment.start_i - index where segment starts
   |  segment.end_i - index where segment ends
  ...
  

• 如果要访问原始系数，请使用段 load_array()方法。它返回两个浮点数和一个 NumPy 数组：

  >>> initial_epoch, interval_length, coefficients = segment.load_array()
  >>> print(coefficients.shape)
  (3, 14080, 13)
  

• 方括号查找机制kernel[3,399]是一种非标准的便利，它只返回文件中的最后一个匹配段。虽然 SPK 标准确实规定最后一段优先，但它也表示对于最后一段未涵盖的日期，应回退特定中心目标对的较早段。因此，如果您曾经处理过复杂的内核，您将需要实施后备规则，将一些日期发送到给定中心和目标的最终段，但将其他日期发送到有资格覆盖它们的早期段。

• 如果您正在考虑轻旅行时间并需要重复计算位置，但最后需要速度，并且希望避免重复昂贵的位置计算，那么试试segment.generate()方法 - 它会让你询问位置，然后只有在确定光时误差现在足够小时才继续速度。

高精度日期

由于所有现代儒略日期都是大于 240 万的数字，因此标准的 64 位 Python 或 NumPy 浮点数必然只留下有限数量的位可用于小数部分。 美国海军天文台天文应用部门的技术说明 2011-02表明， 64 位浮点儒略日期的可能精度约为 20.1 µs。

如果您需要以高于 20.1 µs 的精度提供时间并接收回行星位置，那么您有两种选择。

首先，您可以使用特殊的float96 NumPy 类型提供时间，该类型也别名为longfloat。如果您将float96标量或float96数组作为您的tdb参数提供给任何jplephem例程，您应该得到一个高精度的结果。

其次，您可以将一个或多个日期拆分为两部分，并将它们作为一对参数提供两个tdb和tdb2。如何拆分日期的一种流行方法是使用tdb float 表示整数 Julian 日期，使用tdb2表示指定一天中的时间的分数。如果您希望提供它，几乎所有jplephem例程都接受这个可选 的 tdb2参数，这要感谢 Marten van Kerkwijk 的工作！

支持二进制 PCK

您还可以从二进制 PCK 文件加载和生成旋转矩阵。它的段可通过返回对象的段属性获得。

>>> from jplephem.pck import PCK
>>> p = PCK.open('moon_pa_de421_1900-2050.bpc')
>>> p.segments[0].body
31006
>>> p.segments[0].frame
1
>>> p.segments[0].data_type
2

给定一个太阳系重心儒略日期，该片段将返回构建旋转矩阵所需的三个角度：极点的赤经、极点的赤纬和体轴的累积旋转。通常这些都以弧度表示。

>>> tdb = 2454540.34103
>>> print(p.segments[0].compute(tdb, 0.0, False))
[3.928e-02 3.878e-01 3.253e+03]

您也可以要求速度。

>>> r, v = p.segments[0].compute(tdb, 0.0, True)
>>> print(r)
[3.928e-02 3.878e-01 3.253e+03]
>>> print(v)
[6.707e-09 4.838e-10 2.655e-06]

报告问题

您可以在 GitHub 存储库中报告任何问题、错误或问题：

https://github.com/brandon-rhodes/python-jplephem/

变更日志

2021 年 12 月 31 日 — 版本 2.17

• 修复了摘录命令中的AttributeError 。

2021 年 7 月 3 日 — 2.16 版

• 修复了当星历段需要完全跳过时在摘录命令中引发的ValueError ，因为它与用户指定的日期范围没有重叠。

• 在包的顶层添加了一个__version__常量。

2020 年 9 月 2 日 — 2.15 版

• excerpt子命令现在接受--targets选项以通过仅将匹配的段复制到输出 SPK 文件中来节省空间。

• 儒略日分数tdb2的处理比以前更加小心，当连续时间之间的差异下降约 0.1 µs 时，在连续位置之间提供更平滑的增量。

2020 年 3 月 26 日 — 2.14 版

• 在支持fileno()但不支持mmap()的平台上回退到普通文件 I/O ，例如Pyodide 平台。

2020 年 2 月 22 日 — 2.13 版

• 当一个段被赋予一个超出段日期范围的儒略日期时引发的异常现在是ValueError子类 OutOfRangeError的一个实例，它提醒调用者 SPK 段支持的日期范围，并带有一个数组属性，指示哪些输入日期是有过错。

2019 年 12 月 13 日 — 2.12 版

• 在发现该功能是最近添加的一些用户安装的 NumPy 不支持后，用反向切片[::-1]替换了 NumPy flip()的使用。

2019 年 12 月 13 日 — 2.11 版

• 颠倒切比雪夫多项式的计算顺序以略微提高速度，简化代码，并在一种情况下（将 PCK 输出与 NASA 进行比较）获得部分数字的额外精度。

2019 年 12 月 11 日 — 2.10 版

• 通过新的jplephem.pck模块记录和发布对.bcp二进制 PCK 内核文件的支持。

2019 年 1 月 3 日 — 2.9 版

• 向段类添加了load_array()方法。

2018 年 7 月 22 日 — 2.8 版

• 切换到制作整个文件的单个内存映射，以避免在用户加载具有数百个段的星历时耗尽文件描述符。

2018 年 2 月 11 日 — 2.7 版

• 扩展了命令行工具，最显着的是能够通过 HTTP 仅获取覆盖特定日期范围的大型星历的那些部分，从而生成更小的.bsp文件。

2016 年 12 月 19 日 — 2.6 版

• 修复了使用python -m jplephem从命令行调用模块的能力 ，并添加了一个测试来保持它的修复。

2015 年 11 月 9 日 — 2.5 版

• 将fileno()调用移出DAF构造函数，以支持从StringIO对象中获取至少摘要信息。

2015 年 11 月 1 日 — 2.4 版

• 通过将mmap()从使用 PAGESIZE切换到ALLOCATIONGRANULARITY来添加 Windows 兼容性。

• 通过小心在 NumPy形状元组中仅使用整数来避免新的 NumPy 弃用警告。

• 将名称“TDB”和“TT”添加到 DE430 的名称数据库中。

2015 年 8 月 16 日 — 2.3 版

• 向主DAF类本身添加了对旧 NAIF/DAF 内核（如 de405.bsp ）的自动检测和支持，而不需要笨拙地使用完全不同的替代类。

2015 年 8 月 5 日 — 2.2 版

• 您现在可以从命令行调用jplephem。

• 修复了针对系数计数仅为 2 的 SPK 段引发的异常，例如提供 Mercury 与 Mercury 重心偏移的 DE421 和 DE430 段。

• 支持旧的 NAIF/DAF 内核，例如de405.bsp。

• SPK()构造函数现在更简单，采用DAF对象而不是打开的文件。这被认为是内部 API 更改——公共 API 是构造函数SPK.open()。

2015 年 2 月 24 日 — 2.1 版

• 从一次将整个 SPK 文件映射到内存中切换到按需分别对每个段进行内存映射。

2015 年 2 月 8 日 — 2.0 版

• 添加了对直接从 NASA 下载的 SPICE SPK 内核文件的支持，并将旧的 Python 打包星历指定为“遗留”。

2013 年 11 月 26 日 — 1.2 版

• Helge Eichhorn 修复了position_and_velocity() 参数tdb2的默认值，因此它默认为零天而不是 2.0 天。添加了测试以防止任何未来的回归。

2013 年 7 月 10 日 — 1.1 版

• 弃用旧的compute()方法，支持单独的 position()和position_and_velocity()方法。

• 通过保存由compute_bundle()返回的系数“束”，支持在两个单独的阶段计算位置和速度。

• 来自 Marten van Kerkwijk：第二个tdb2时间参数，适用于希望从两个 64 位浮点数构建更高精度日期的用户。

2013 年 1 月 18 日 — 1.0 版

• 初始发行

参考

喷气推进实验室的“太阳系动力学”页面介绍了进行太阳系位置计算的各种选项： http ://ssd.jpl.nasa.gov/?ephemerides

构建jplephem Python 星历包的纯 ASCII 格式元素集以及文档可在以下位置找到：ftp: //ssd.jpl.nasa.gov/pub/eph/planets/ascii/

使用星历表的等效 FORTRAN 代码可在同一 FTP 站点找到：ftp: //ssd.jpl.nasa.gov/pub/eph/planets/fortran/