算法系列之十九：用天文方法计算日月合朔（新月）

        中国农历的朔望月是农历历法的基础，而朔望月又是严格以日月合朔发生的那一天作为月首，因此日月合朔时间的计算是制定农历历法的关键。本文将介绍ELP-2000/82月球运行理论，以及如何用ELP-2000/82月球运行理论计算日月合朔时间。

        要计算日月合朔时间，首先要对日月合朔这一天文现象进行数学定义。朔望月是在地球上观察到的月相周期，平均长度约等于29.53059日，而恒星月（天文月）是月亮绕地球公转一周的时间，长度约27.32166日。月相周期长度比恒星月长大约两天，这是因为在月球绕地球旋转一周的同时，地球还带着它绕太阳旋转了一定的角度的缘故，所以月相周期不仅与月球运行有关，还和太阳运行有关。日月合朔的时候，太阳、月亮和地球三者接近一条直线，月亮未被照亮的一面对着地球，因此地球上看不到月亮，此时又被称为新月。图（1）就是日月合朔天文现象的示意图：

图（1）日月天文现象示意图

月亮绕太阳公转的白道面和地球绕太阳公转的黄道面存在一个最大约5°的夹角，因此大多数情况下，日月合朔时都不是严格在同一条直线上，不过也会发生在同一直线的情况，此时就会发生日食。图（1－b）显示了日月合朔时侧切面上月亮的三种可能的位置情况，当月亮处在位置2时就会发生日食。由图（1）可知，日月合朔的数学定义就是太阳和月亮的地心视黄经差为0的时刻。

        要计算日月合朔，需要知道太阳地心视黄经和月亮地心视黄经的计算方法。“日历生成算法”系列的第三篇《用天文方法计算二十四节气》一文已经介绍了如何用VSOP82/87行星理论计算太阳的地心视黄经，本文将继续介绍如何用ELP-2000/82月球理论计算月亮的地心视黄经。ELP-2000/82月球理论是M. Chapront-Touze和J. Chapront在1983年提出的一个月球位置的半解析理论，和其它半解析理论一样，ELP-2000/82理论也包含一套计算方法和相应的迭代周期项。这套理论共包含37862个周期项，其中20560个用于计算月球经度，7684个用于计算月球纬度，9618个用于计算地月距离。但是这些周期项中有很多都是非常小的值，例如一些计算经纬度的项对结果的增益只有0.00001角秒，还有一些地月距离周期项对距离结果的增益只有0.02米，对于精度不高的历法计算，完全可以忽略。

        有很多基于ELP-2000/82月球理论的改进或简化理论，《天文算法》一书的第四十五章就介绍了一种改进算法，其周期项参数都是从ELP-2000/82理论的周期项参数转换来的，忽略了小的周期项。使用该方法计算的月球黄经精度只有10”，月亮黄纬精度只有4”，但是只用计算60个周期项，速度很快，本文就采用这种修改过的ELP-2000/82理论计算月亮的地心视黄经。这种计算方法的周期项分三部分，分别用来计算月球黄经，月球黄纬和地月距离，三部分的周期项的内容一样，由四个计算辐角的系数和一个正弦（或余弦）振幅组成。计算月球黄经和月球黄纬使用正弦表达式求和：A * sin(θ)，计算地月距离用余弦表达式求和：A * cos(θ)，其中辐角θ的计算公式是：

θ = a * D + b * M + c * M’ + d * F                           (4.1式)

4.1式中的四个辐角系数a、b、c和d由每个迭代周期项给出，日月距角D、太阳平近地角M、月亮平近地角M’以及月球生交点平角距F则分别有4.2式－4.5式进行计算：

D = 297.8502042 + 445267.1115168 * T - 0.0016300 * T2

+ T3 / 545868 - T4 / 113065000                   (4.2式)
M = 357.5291092 + 35999.0502909 * T - 0.0001536 * T2

+ T3 / 24490000                                (4.3式)
M' = 134.9634114 + 477198.8676313 * T + 0.0089970 * T2

+ T3 / 69699 - T4 / 14712000                     (4.4式)
F = 93.2720993 + 483202.0175273 * T - 0.0034029 * T2

- T3 / 3526000 + T4 / 863310000                  (4.5式)

以上各式计算结果的单位是度，其中T是儒略世纪数，T计算由4.6式计算：

T = (JDE - 2451545.0) / 36525.0                          (4.6式)

以计算月球黄经的周期项第二项的计算为例，第二项数据如下，辐角系数a = 2，b = 0，c = -1，d = 0，振幅A = 1274027，黄经计算用正弦表达式，则I2的计算如下所示：

I2 = 1274027 * sin(2D – M’)                            （4.7式）

在套用4.7式计算出60个月球黄经周期项值的时候，需要注意对包含了太阳平近地角M的项进行修正，因为M的值与地球公转轨道的离心率有关，因为离心率是个与时间有关的变量，导致振幅A实际上是个变量，需要根据时间进行修正。月球黄经周期项的修正方法是：如果辐角中包含了M或-M时，需要乘以系数E修正；如果辐角中包含了2M或-2M，则需要乘以系数E的平方进行修正。系数E的计算表达式如下：

E = 1 - 0.002516 * T - 0.0000074 * T2                      （4.8式）

其中T值由4.6式计算。上面的计算月球黄经的第二个周期项中M对应的系数是0，因此I2不需要修正，但是第五个周期项中M对应的系数是1，因此I5需要乘以E进行修正。套用4.7式计算出60个月球黄经周期项值I1－I60之和ΣI：

ΣI = I1 + I2 + … + I60                                    （4.9式）

        月球黄纬的周期项和Σb的计算方法与月球黄经周期项和ΣI的计算方法一样，每个月球黄纬周期项也包含振幅A和四个辐角系数a、b、c和d，对于太阳平近地角M的系数b不是0的情况也需要乘以E或E2进行修正。地月距离的周期项和Σr也可以按照上面的方法计算，计算地月距离的目的是为了计算月亮光行差，因为地月距离较小，从地球观察月亮产生的光行差也很小，相对于本文算法的精度（月球黄经精度10”，月亮黄纬精度4”）来说，可以忽略光行差修正，因此就不用计算地月距离。

        由于金星和木星对月球的摄动影响，需要对计算出的月球黄经周期项和ΣI和月球黄纬周期项和Σb金星摄动修正，修正的方法如下：

ΣI += +3958 * sin( A1 ) + 1962 * sin( L' - F ) + 318 * sin( A2 )             （4.10式）

Σb += -2235 * sin( L' ) + 382 * sin( A3) + 175 * sin( A1 - F ) + 175 * sin( A1 + F )
       + 127 * sin( L' - M') - 115 * sin( L' + M')                           （4.11式）

其中M’和F分别由4.4式和4.5式计算得到，L’是月球平黄经，计算方法是：

L'=218.3164591 + 481267.88134236 * T - 0.0013268 * T2

+ T3 / 538841 - T4 / 65194000                         （4.12式）

A1、A2和A4是摄动角修正量，计算方法如下：

A1 = 119.75 + 131.849 * T                                             （4.13式）
A2 = 53.09 + 479264.290 * T