用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)=2t{T}_{i-1}(t)-T_{i-2}(t)$$$$t\in [-1,1],i>=2$$
  $$T_0(t)=0$$$$T_1(t)=1$$$$T_i(t)=2T_{i-1}(t)+2t{T}_{i-1}^{\prime }(t)-T_{i-2}^{\prime }(t)$$$$t\in [-1,1],i>=2$$

具体方法:
之前我们提到的GROUP 1050的三行数据，分别记为data(1,#)、data(2,#)、data(3,#)（“#”位置是天体的序号）
设定数据的起止时间分别为 t1 和 t2。
对于第 i 个天体，所计算的时刻为 tjd。每个子块的跨度为 $$dt(i)=\frac{t2-t1}{data(3,i)}$$那么 tjd 所在的子块序数 j 为$$j=取整[\frac{tjd-t1}{dt(i)}]+1$$在子块 j 内，tjd 相应的规格化时间为$$t=\frac{2[tjd-t1-(j-1)dt(i)]}{dt(i)}-1$$最后，第 i 个天体（含章动、天平动）在历元的第k个坐标量、速度量分别为$${\Sigma }_{m=m1}^{m2}{T}_{m^{\prime }}(t)C_m$$$${\Sigma }_{m=m1}^{m2}{T}_{m^{\prime }}^{\prime }(t)C_m$$其中，m1、m2分别是你使用的数据开始和截止的位置，（可以是一个子块内该天体的该坐标的数据的起止位置）；m’=(m-m1)；Cm是m位置上的数据。

讲到这里，您应该可以得到所有天体的数据了，用它来做什么样的工作就看您的个人需求了。

################################################################

【MMA内容】

接下来的内容是使用历表的实例，基于mathematica，会有很多mma特有的计算方法，并不通用于fortran、c++等语言，各位各取所需吧。
还要指出，这里的算法介绍不包括基础语法，初学的小伙伴可以把代码拿走自己看帮助文档，写的都很详细的哦。

三、对数据的导入和使用（以金星为例）

1、数据导入

本人在这里使用了ReadList这个函数来导入数据，有关导入数据的函数，可以在帮助文档里搜索“列表”或者“文件”。以数字形式导入数据，并形成列表，赋值在de405上

de405=ReadList[“C:\\ascp1600.405”,Number]；

对数据进行切分，把每个数据块单独赋值在变量上。注意，我在划分数据块的时候，每个数据块有1022个元素，而不是1018，这是因为每一个数据块前面都有一行标题，是序号+1018，比如第15个块就是“15 1018”，而结尾处则有两个数据是零，可能是为了把数据补齐，这样的话一个数据块就有1022个数据

data[i_] := Take[de405, {1022 i - 1021, 1022 i}];

从每一个数据块中提取金星的数据。注意，开头的两个数据是不用的，所以虽然从header里读到的金星数据的位置是171~ 230，实际上应该是173~ 232。

VenusData[i_] := Take[data[i], {171,232}];

为了操作简便，我们把所有的金星数据合在一起，形成一个列表。这样就得到了金星在20年内的全部数据。229的意思是，这个文件里一共有229个数据块。（您自己用20年除以32天就得到了）

VenusDataAll = Flatten[Table[VenusData[i], {i, 1, 229}]];

然后为了计算需要，我们要把每一个子块的数据分别列出来

VenusData8day[i_] := Take[VenusDataAll,{39 i - 38, 39 i}];

再从每一个子块中提取三个坐标各自的数据

VenusData8dayx[i_] := Take[VenusData8day[i], {1, 14}];
VenusData8dayy[i_] := Take[VenusData8day[i], {15, 28}];
VenusData8dayz[i_] := Take[VenusData8day[i], {29, 42}];

然后把每个坐标的的数据分别合成一个大列表。916的意思是一共分出来916个子块，您可以自己用4乘以229.

VenusDatax = Flatten[Table[VenusData8dayx[i], {i, 1, 916}]];
VenusDatay = Flatten[Table[VenusData8dayy[i], {i, 1, 916}]];
VenusDataz = Flatten[Table[VenusData8dayz[i], {i, 1, 916}]];

这样，您就拥有了金星每个坐标在20年内的所有数据

2、插值

根据上面说过的插值方法，直接插值就能够得到。mma内置了切比雪夫多项式，可以直接调用。

t1= 2305424.5; t2 = 2312752.5; (开始和结束的儒略日)
dt = 8; (子块间的跨度)
j[t_] := IntegerPart[(t - t1)/dt] + 1; (子块序数)
tt[t_] := 2 (t - t1 - dt (j[t] - 1))/dt - 1; (规格化时间)
xm1[t_] := 14 (j[t]) - 13; xm2[t_] := 14 j[t]);
ym1[t_] := 14 (j[t]) - 13; ym2[t_] := 14 j[t]);
zm1[t_] := 14 (j[t]) - 13; zm2[t_] := 14 j[t]); (切比雪夫系数数据的位置)
Venusx[t_] := Sum[ChebyshevT[xm - xm1[t], tt[t]] Take[ VenusDatax, {xm, xm}], {xm, xm1[t], xm2[t]}];
Venusy[t_] := Sum[ChebyshevT[ym - ym1[t], tt[t]] Take[ VenusDatay, {ym, ym}], {ym, ym1[t], ym2[t]}];
Venusz[t_] := Sum[ChebyshevT[zm - zm1[t], tt[t]] Take[ VenusDataz, {zm, zm}], {zm, zm1[t],
zm2[t]}]; (地E月M系质心的三个坐标，对儒略日t的函数)

最终得到三个坐标对时间的函数，Venusx[t],Venusy[t],Venusz[t]，这三个插值的来的函数是连续的，后面可以直接使用。

简单应用与测试，可以用如下代码画出三维图像，时间就选取其中的三百天。这里就不贴图了，感兴趣的小伙伴自己去画吧嘿嘿嘿。

ParametricPlot3D[ Flatten@{EMx[t], EMy[t], EMz[t]}, {t, 2305425, 2305725}]

四、天体信息的应用实例

这一部分是我本人自己的作业，分享出来大家一起学习

（一）水星西大距时刻

1、数据采集（前面都讲过了，举一反三）

水星的三个坐标分别是Mex(t)、Mey(t)、Mez(t)
地球的三个坐标分别是Eax(t)、Eay(t)、Eaz(t)

这里有一点值得说明，做这个项目需要先把地球的运动从地月系质心的运动中分离出来。由于这种分离不涉及坐标系的相对转动，所以非常容易，就不多赘述了。

2、大距的判定

本人采用的方法是，写出水星和太阳对地心的夹角对时间的函数，求极大值得到西大距时刻。

先写出地球Earth到水星Mercury的距离EaMe[t]

EaMe[t_] := ((Mex[t] - Eax[t])^ 2 + (Mey[t] - Eay[t]) ^ 2 + (Mez[t] - Eaz[t])^ 2)^(1/2)

再分别写出水星和地球到太阳的距离

rEa[t_] := (Eax[t]^2 + Eay[t]^2 + Eaz[t]^ 2)^(1/2);
rMe[t_] := (Mex[t]^2 + Mey[t]^2 + Mez[t]^ 2)^(1/2);

根据余弦定理表示水星和太阳对地心的夹角

θ[t_] := ArcCos[(EaMe[t]^2 + rEa[t]^2 - rMe[t]^2)/(2 EaMe[t] rEa[t])];

然后就是一阶导为零、二阶导为负，求极大值的点。不过，很可能因为数据点太多而导致程序慢，这样的话我们就采用离散的判定方法，比较某一时刻 θ（t）与前一个和后一个，它最大他就是极大值。

基于这样的思路，我们可以设定测试点的间隔为1/10天

deltat=1/10;

可以生成一系列测试点（时间跨度您自己来修改）

TestPoints = Table[t, {t, 2305425, 2305725, deltat}];

可以自定义判定函数

greatestQ[t_] := If[First@θ[t - deltat] <= First@θ[t] && First@θ[t + deltat] <= First@θ[t], t, Nothing];

把判定函数作用在测试点上就能够得到水星西大距发生的时刻，精度为1/10天

greatestQ /@ TestPoints;

可能运行时间比较长，需要耐心等一下。

（二）金星凌日

1、数据采集 不再赘述

2、判定条件

我在这里采用的判定方法是区域交集法。有天文基础的朋友可以知道，金星凌日的区域在宇宙空间中是一个锥体，如果地球的一部分在这个锥体内，人们就能够观测到金星凌日。所以只要我们能够找到这个锥体，再利用mma神奇的判定区域交集的内置函数，我们就能够得到某一时刻是否发生金星凌日。一步一步来

（1）凌日区域

凌日区域是空间中的圆锥形。圆锥的顶点由太阳半径和金星半径的比值决定。我们利用mma内置的恒星和行星数据获得半径信息。

SunRadius = QuantityMagnitude@StarData[“Sun”, “Radius”];
VenusRadius = QuantityMagnitude@PlanetData[“Venus”, “Radius”];
EarthRadius = QuantityMagnitude@PlanetData[“Earth”, “Radius”];

锥体的顶点离太阳的距离，和金星离太阳的距离，的比值设定为λ

λ= SunRadius/(SunRadius + VenusRadius);

于是可以得到锥体顶点的坐标

ConePoint[t_] := [Lambda] Flatten@{Vex[t], Vey[t], Vez[t]};

在mma里面绘制锥形区域需要顶点、底面中心、底面半径。显然底面中心在金星矢径所在直线上，我们把它的距离设置为 h 倍的金星-日心距，也就是说

Point[t_]:=h {Vex[t], Vey[t], Vez[t]}

同时根据中学几何知识相似性可以算出这个距离上，底面的半径应该是

r=VenusRadius (h - λ)/(1 - λ)

这样就能够构建凌日区域

SelectRegion1[t_] := Region[Cone[{Flatten@Point[t], ConePoint[t]}, r];

还要构建地球区域

EarthRegion[t_] := Region[Sphere[Flatten@{Eax[t], Eay[t], Eaz[t]}, EarthRadius]];

2、分级判定算法

众所周知，金星凌日的发生周期特别长，如果直接用标准的凌日区域和地球区域求交集，那就需要一次性遍历上万个时间点，这样的计算量很大，而且计算时间很长。对此，我选择了不是很巧的办法，分级筛选

先建立较大的凌日范围，比如把上面的 r 改成一百倍，这样即允许我们选取点的时间跨度小一点、测试点就能少一点。但是一定要注意，如果时间跨度过大，导致两个点之间地球走过的范围，比凌日区域大，就有可能存在漏点的情况，所以要适当选取凌日区域和时间间隔。

我们构建多级的凌日区域

SelectRegion1[t_] := Region[Cone[{Flatten@Point[t], ConePoint[t]}, r];
SelectRegion3[t_] := Region[Cone[{Flatten@Point[t], ConePoint[t]},3 r];
SelectRegion10[t_] := Region[Cone[{Flatten@Point[t], ConePoint[t]},10 r];
SelectRegion50[t_] := Region[Cone[{Flatten@Point[t], ConePoint[t]}, 50r];
SelectRegion200[t_] := Region[Cone[{Flatten@Point[t], ConePoint[t]},200 r];

对于每一级筛选，都要构建适当的测试点列表。我们这里只展示第一级判定的情形，后续的判定由您各位自行举一反三。这里的起止时间是我随手写的，如果要使用请自己修改一下。精度定为50天，这是适用于200倍区域的。

TestPointsLevel1= Table[t,{t,2305425, 2307425,50}];

定义判定函数。这里多说一句，RegionDisjoint是用来判定两个区域没有交集的，如果没有交集就返回True，有交集则返回False

transitQ1[t_] := If[! RegionDisjoint[EarthRegion[t], SelectRegion200[t]], t, Nothing];

把判定函数作用在测试点上，就能够得到可能会发生凌日的时间。

TestResultLevel1= transitQ1 /@ TestPointsLevel1;

对于这个结果，必须要注意这不是凌日的日期，真正的凌日会发生在其前后50天的范围内，所以下一级筛选一定要包含这次结果前后50天的范围。

3、结果

按照上述方法，经过多级筛选后就能够得到准确地凌日时间，理论上精度可以很小很小，仅仅受制于切比雪夫插值法的精度。

五、最后的话

因为本人也是第一次学着使用历表，学习过程中踩了很多坑，所以希望下写这篇文章能够帮助大家。自我感觉这篇文章能解决一些问题，不过还是很稚嫩，算法很笨拙、叙述也很繁琐，非常欢迎各位读者能够给我提些意见和建议。如果有问题请尽管提出，我看到后会尽力解答的。