用mathematica调用JPL的DE系列星历表 · DE405

【写在前面的话】

本人是师大天文系的学生，这学期课业要求学习DE系列星历表的使用，老师没有细讲具体方法，所以个人摸索了大概两个礼拜，也算是勉强实现了目标。写这篇文章的目的是为了给同样小白的朋友一点引导，可以少走些弯路提高效率。这里必须要说明，这些知识都是前辈老师们从jpl发布的信息、以及他们的实践中总结出来的，我只是知识的的搬运工而已。

这篇文章分为两大块，第一部分讲历表数据的内容和使用方法，第二部分主要介绍mathematica实现历表调用的具体操作，对于mma没有需求的朋友可以只关注前半部分。

一、DE系列历表的结构（以DE405为例）

虽然405比较老了，但是就其数据量和精度来说，还是能满足大多数需求的，而且所有DE历表的结构是一样的，这里就讲405这一种了。

1、DE星历表的获取

因为de历表是美国人的杰作，最官方的数据是在jpl的网站上，所以国内的小伙伴可能翻墙获取，网站在这里 http://ssd.jpl.nasa.gov/ ，好像社区里也有历表数据可以下载，这里就不多说了
（1）各位下载到的文件应当包括这几个文件夹：ascii、fortran、linux等等。jpl给了读历表的fortran、c++程序，您可以选用这些，相关的文章在社区里也能够找到，这里给一篇的链接吧
（2）真正的天体数据是ascii这个文件夹，里面应该包括：名为de***的历表文件夹、ascii_format.txt。这个txt文件就描述了历表文件里面的信息，您可以自己去看，不过这些内容我也会在后面的结构讲解里包括。
（3）对于de***文件夹，我们选用de405具体讲解。里面包括


可以看到，每一个数据文件包含了20年的信息，最后是header、和testpo。文件的扩展名是.405，不过直接用记事本打开也是能看到数据的。

2、header文件结构

历表文件的数据结构基本包含在header的描述里，可以用记事本打开看一下，对照我下面写的内容。##必须要先说明##，本人对历表的了解来源于这篇文章，我们在这里也算是重复其内容。

（1）“KSIZE= 2036 NCOEFF= 1018”

NCOEFF的含义是，文件中每个数据块中切比雪夫插值系数的个数。KSIZE是NCOEFF的2倍

（2）GROUP 1010

包含了历表名称及版本，以及分别以儒略日和公历表示的历表起始和终止历元信息。

（3）GROUP 1030

以儒略日的形式给出了历表的起止时间，即2305424.5~2525008.5
还有每个数据块的时间跨度：32天
（不过，header里的起止历元和数据文件里可能不一样，这个要以数据为准）

（4）GROUP 1040 & GROUP 1041

这两块说的是同一个东西。1040列出了一系列天文常数的名称，一共156个，1041给出了对应常数的数值，常数的单位可以在ascii_format.txt里面找到。我们后面再介绍

（5）GROUP 1050

这一块给出的是每个天体对应数据的位置。可以看到是有三行数字，每行13个。读的时候要按列读，第 i 列的三个数从上到下分别描述了：第 i 个天体数据开始的位置、第 i 个天体每个子块中每个坐标分量的数据量、数据块划分为子块的个数。

• 第 i 个天体数据开始的位置 ： 比如第一个天体（水星）从第3个数据开始读，读到第170个；而第二个天体（金星）从第171个数据开始读，读到第230。为什么不是从1开始？因为前两个数据是该数据块的起止时间，依然是儒略日。可以容易计算这个起止时间的跨度是32天
• 第 i 个天体每个坐标分量的数据量 ： 拿水星为例，水星有三个坐标xyz，每个坐标有14个数据，一共42个数据
• 子块个数 ： 同样以水星为例，4代表了32天要分为四块，相当于每个子块包含了8天的数据，而上面说到每个子块有42个数据，那么这个数据块里的水星数据就是168个，也就是从3到170.

（6）GROUP 1070

可以看到这个下面没东西，不是文件坏了。jpl希望你在用他的程序的时候，先把几个文件合并起来，把header放在第一个，就相当于合并完之后，GROUP 1070下面是真正的天体数据。

3、ascii_format.txt

这个文件的主要内容就是讲解header信息的含义。其他的内容还包括：

• 天体序号：
  1 Mercury 水星
  2 Venus 金星
  3 Earth-Moon barycenter 地月系质心
  4 Mars 火星
  5 Jupiter 木星
  6 Saturn 土星
  7 Uranus 天王星
  8 Neptune 海王星
  9 Pluto 冥王星
  10 Moon (geocentric) 地心坐标系下的月球
  11 Sun 太阳
  12 Earth Nutations in longitude and obliquity (IAU 1980 model) 地球经线章动和倾斜
  13 Lunar mantle libration 月面天平动
  14 Lunar mantle angular velocity 月面角速度
  15 TT-TDB (at geocenter)

可以看到有15个天体，但是旧版本历表往往只有13个，这个就看各位的个人需求了

• 坐标数
  如无特殊说明都是三个坐标分量；12号只有两个分量、15号只有1个分量。

• 单位
  位置量都是km，章动和天平动都是rad，月面角速度rad/day，TT-TDB单位是秒。

• 参考系说明
  除了月球运动数据是地心参考系，其他都是日心系，即标准的国际天球参考系ICRS；历表并不直接包含地球的运动，需要得到日心系下的地月系质心运动、以及地心系下的月球运动，做一下相对运动与坐标转换，才能得到地球的运动。

4、数据顺序

请注意，前面说过（拿水星为例）每个坐标有14个数据，一个子块里有42个数据，但是要注意数据的排列。它的顺序是：x1,x2,x3……y1,y2,y3……z1,z2,z3，四个子块各自都是这样排列，然后它们再接在一起。

二、插值方法

现在，给为已经知道数据是什么和数据在哪里了，但是想要用这些数据得到天体的轨道，还需要进行插值

• 切比雪夫多项式及其一阶导数：
  $$T_0(t)=1$$$$T_1(t)=t$$ T i ( t ) = 2 t T i − 1 ( t ) − T i − 2 ( t ) T_i(t)=2tT_{i-1}(t)-T_{i-2}(t) Ti​(t)=2tT