jplephem 2.22

安装

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

$ pip install jplephem

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

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

大多数用户会使用jplephem与NASA喷气推进实验室（JPL）的NAIF设施提供的卫星行星核（SPK）文件一起使用，这些文件与自己的SPICE工具包一起使用。他们在以下目录下收集了他们最有用的内核：

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

要了解更多关于SPK文件的信息，可以在NASA JPL域下的NAIF设施网站上找到官方的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非常大，近1GB。但如果您只需要几个月内木星的卫星数据，您可以下载大量的数据

$ 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星历是一个有用的起点。它的重量为17MB，但提供了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:
1899-07-29..2053-10-09  Type 2  Solar System Barycenter (0) -> Mercury Barycenter (1)
1899-07-29..2053-10-09  Type 2  Solar System Barycenter (0) -> Venus Barycenter (2)
1899-07-29..2053-10-09  Type 2  Solar System Barycenter (0) -> Earth Barycenter (3)
1899-07-29..2053-10-09  Type 2  Solar System Barycenter (0) -> Mars Barycenter (4)
1899-07-29..2053-10-09  Type 2  Solar System Barycenter (0) -> Jupiter Barycenter (5)
1899-07-29..2053-10-09  Type 2  Solar System Barycenter (0) -> Saturn Barycenter (6)
1899-07-29..2053-10-09  Type 2  Solar System Barycenter (0) -> Uranus Barycenter (7)
1899-07-29..2053-10-09  Type 2  Solar System Barycenter (0) -> Neptune Barycenter (8)
1899-07-29..2053-10-09  Type 2  Solar System Barycenter (0) -> Pluto Barycenter (9)
1899-07-29..2053-10-09  Type 2  Solar System Barycenter (0) -> Sun (10)
1899-07-29..2053-10-09  Type 2  Earth Barycenter (3) -> Moon (301)
1899-07-29..2053-10-09  Type 2  Earth Barycenter (3) -> Earth (399)
1899-07-29..2053-10-09  Type 2  Mercury Barycenter (1) -> Mercury (199)
1899-07-29..2053-10-09  Type 2  Venus Barycenter (2) -> Venus (299)
1899-07-29..2053-10-09  Type 2  Mars Barycenter (4) -> Mars (499)

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

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

文件的每个段都允许您根据给定的儒略日预测一个天体相对于另一个天体的位置。提供了一个小的例程，将公历日期转换为儒略日

>>> from jplephem.calendar import compute_julian_date
>>> compute_julian_date(2015, 2, 8)
2457061.5

以下是计算火星（目标4）相对于太阳系质心（目标0）在2015年2月8日午夜TDB（质心动力时间）坐标的方法，使用我们刚刚计算的儒略日

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

相比之下，学习火星相对于地球的位置需要三个步骤：从火星到太阳系质心，到地球-月球质心（3），最后到地球本身（399）。

>>> 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]

您可以看到这个星历DE421的输出是以千米为单位的。如果您使用另一个星历，请检查其文档以确保它使用的单位。

如果您提供日期作为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时，它会检查文件并创建一个段对象列表，该列表被保存在名为k.segments的属性下，您可以在自己的代码中通过遍历它来查看。

• 关于每个段的信息比您打印SPK文件时获得的单行摘要更为详细，您可以通过要求段详细打印自身来查看这些信息。

  >>> segment = kernel[3,399]
  >>> print(segment.describe())
  1899-07-29..2053-10-09  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微秒。

如果您需要提供时间并以超过20.1微秒的精度返回行星位置，那么您有两个选择。

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

其次，您可以先将您的日期或日期拆分为两部分，然后将它们作为一对参数 tdb 和 tdb2 提供。一种拆分日期的流行方法是将 tdb float 用于整数儒略日，将 tdb2 用于指定一天中时间的分数。几乎所有 jplephem 程序都接受这个可选的 tdb2 参数，感谢 Marten van Kerkwijk 的工作！

二进制 PCK 支持

您还可以从二进制 PCK 文件加载并生成旋转矩阵。这些段可以通过返回对象的 segments 属性访问。

>>> 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]

关闭历书

要释放与历书相关的所有打开的文件和内存映射，请调用其 close() 方法。

>>> kernel.close()
>>> p.close()

报告问题

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

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

变更日志

2024 年 4 月 24 日 — 版本 2.22

• 当打印时，段现在使用公历而不是打印原始儒略日来打印其开始和结束日期。

• 现在提供了一个小的 compute_julian_date 程序，用于将日历日期转换为儒略日。

• 修复了当 PCK 段 compute() 方法接收到超出范围的日期时引发的 ValueError 文本；它错误地报告了不正确的儒略日期范围，因为 PCK 使用 J2000 之前或之后的秒来计算时间，而不是年份。

2023 年 12 月 1 日 — 版本 2.21

• 调整了一个导入，以避免在 Python 2 中发生致命异常，以防有人仍在使用它。

2023 年 11 月 13 日 — 版本 2.20

• 现在每个段都受到一个锁的保护，以防两个线程同时触发执行段数据的初始加载的代码；症状是罕见异常 ValueError: cannot reshape array。

2023 年 9 月 6 日 — 版本 2.19

• 修复了 excerpt 命令中的错误，该错误导致当输入历书有二十多个段时截断其输出。现在该命令的输出应包括来自非常大的历书的所有匹配段。

• 修复了 excerpt 命令，以便在命令行上指定的日历日期产生以分数 .5 结尾的儒略日，这使得摘录端点更精确。

2022 年 9 月 28 日 — 版本 2.18

• 添加了对大端处理器的支持，并创建了一个 GitHub Actions CI 构建，其中包含大端和小端架构。

2021 年 12 月 31 日 — 版本 2.17

• 修复了 excerpt 命令中的 AttributeError。

2021 年 7 月 3 日 — 版本 2.16

• 修复了在历书段需要完全跳过因为它与用户指定的日期范围没有重叠时在 excerpt 命令中引发的 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）获得额外的一位精度，反转了 Chebyshev 多项式的计算顺序。

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 shape 元组中使用，避免了新的 NumPy 弃用警告。

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

2015 年 8 月 16 日 — 版本 2.3

• 将自动检测和对旧 NAIF/DAF 内核（如 de405.bsp）的支持添加到主 DAF 类本身，而不是要求使用一个完全不同的替代类。

2015 年 8 月 5 日 — 版本 2.2

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

• 修复了对于像 DE421 和 DE430 段（提供水星相对于水星质心的偏移）这样的只有 2 个系数的 SPK 段抛出异常的问题。

• 支持旧 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 的默认值，使其默认为0天而不是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](http://ssd.jpl.nasa.gov/?ephemerides)

在相同的FTP站点上可以找到使用星历的等效FORTRAN代码：[ftp://ssd.jpl.nasa.gov/pub/eph/planets/fortran/](ftp://ssd.jpl.nasa.gov/pub/eph/planets/fortran/)