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

任万廷 2015-01-04

展开全文

中国农历的朔望月是农历历法的基础，而朔望月又是严格以日月合朔发生的那一天作为月首，因此日月合朔时间的计算是制定农历历法的关键。本文将介绍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 * T²

+ T³ / 545868 - T⁴ / 113065000 (4.2式)
M = 357.5291092 + 35999.0502909 * T - 0.0001536 * T²

+ T³ / 24490000 (4.3式)
M' = 134.9634114 + 477198.8676313 * T + 0.0089970 * T²

+ T³ / 69699 - T⁴ / 14712000 (4.4式)
F = 93.2720993 + 483202.0175273 * T - 0.0034029 * T²

- T³ / 3526000 + T⁴ / 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，黄经计算用正弦表达式，则I₂的计算如下所示：

I₂ = 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 * T² （4.8式）

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

ΣI = I₁+ I₂ + … + I₆₀ （4.9式）

月球黄纬的周期项和Σb的计算方法与月球黄经周期项和ΣI的计算方法一样，每个月球黄纬周期项也包含振幅A和四个辐角系数a、b、c和d，对于太阳平近地角M的系数b不是0的情况也需要乘以E或E²进行修正。地月距离的周期项和Σ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 * T²

+ T³ / 538841 - T⁴ / 65194000 （4.12式）

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

A1 = 119.75 + 131.849 * T                                             （4.13式）
A2 = 53.09 + 479264.290 * T                                           （4.14式）
A3 = 313.45 + 481266.484 * T                                          （4.15式）

完成所有修正后，就可以用4.16式和4.17式最终得到月亮的地心视黄经和地心视黄纬：

λ = L'+ ΣI / 1000000.0 （4.16式）

β = Σb / 1000000.0 （4.17式）

ΣI和Σb最后要除以1000000.0是因为周期项系数中振幅A的单位是0.000001度，因此λ和β的单位是度。下面给出计算月球地心视黄经的代码：

123 double GetMoonEclipticLongitudeEC(double dbJD)
124 {
125     double Lp,D,M,Mp,F,E;
126     double dt = (dbJD - JD2000) / 36525.0; /*儒略世纪数*/
127 
128     GetMoonEclipticParameter(dt, &Lp, &D, &M, &Mp, &F, &E);
129 
130     /*计算月球地心黄经周期项*/
131     double EI = CalcMoonECLongitudePeriodicTbl(D, M, Mp, F, E);
132 
133     /*修正金星,木星以及地球扁率摄动*/
134     EI += CalcMoonLongitudePerturbation(dt, Lp, F);
135 
136     /*计算月球地心黄经*/
137     double longitude = Lp + EI / 1000000.0;
138 
139     /*计算天体章动干扰*/
140     longitude += CalcEarthLongitudeNutation(dt / 10.0);
141 
142     longitude = Mod360Degree(longitude);
143     return longitude;
144 }

函数参数dbJD是力学时儒略日时间，返回以度为单位的月球视黄经。其中GetMoonEclipticParameter()函数分别根据4.2式、4.3式、4.4式、4.5式、4.8式和4.12式计算日月距角D、太阳平近地角M、月亮平近地角M’、月球生交点平角距F、修正系数E和月球平黄经L’，不需多说明，只要根据以上各式直接计算即可。CalcMoonECLongitudePeriodicTbl()函数计算60个月球黄经周期项和，并根据M值系数的情况进行修正，算法实现如下：

42 double CalcMoonECLongitudePeriodicTbl(double D, double M, double Mp, double F, double E)
43 {
44     double EI = 0.0 ;
45 
46     for(int i = 0; i < COUNT_ITEM(Moon_longitude); i++)
47     {
48         double sita = Moon_longitude[i].D * D + Moon_longitude[i].M * M + Moon_longitude[i].Mp * Mp + Moon_longitude[i].F * F;
49         sita = DegreeToRadian(sita);
50         EI += (Moon_longitude[i].eiA * sin(sita) * pow(E, fabs(Moon_longitude[i].M)));
51     }
52 
53     return EI;
54 }

CalcMoonLongitudePerturbation()函数计算月球黄经摄动修正量，使用了4.13式和4.14式给出的计算方法：

87 double CalcMoonLongitudePerturbation(double dt, double Lp, double F)
88 {
89     double T = dt; /*T是'ca?从'b4?J2000起'c6?算'cb?的'b5?儒'c8?略'c2?世'ca?纪'bc?数'ca?*/
90     double A1 = 119.75 + 131.849 * T;
91     double A2 = 53.09 + 479264.290 * T;
92 
93     A1 = Mod360Degree(A1);
94     A2 = Mod360Degree(A2);
95 
96     double result = 3958.0 * sin(DegreeToRadian(A1));
97     result += (1962.0 * sin(DegreeToRadian(Lp - F)));
98     result += (318.0 * sin(DegreeToRadian(A2)));
99 
100     return result;
101 }

至此，本文已经介绍了使用ELP-2000/82月球理论计算任意时刻月亮地心视黄经的方法，结合“日历生成算法”系列的第三篇《用天文方法计算二十四节气》一文介绍的计算太阳地心视黄经的方法，就可以计算日月合朔的准确时间了。由于ELP-2000/82月球理论也没有根据月球黄经反算时间的方法，因此本文也采用和《用天文方法计算二十四节气》一文中一样的牛顿迭代法计算日月合朔时间。

关于牛顿迭代法可以参考相关的数学资料，“日历生成算法”系列的第三篇《用天文方法计算二十四节气》一文对如何使用牛顿迭代法有简单的介绍，可以参考一下。总的来说，就是要先定义需要求解的方程f(x)，根据上文的介绍，我们需要求解的是太阳的地心黄经和月亮的地心黄经差值是0的时候的时间，《用天文方法计算二十四节气》一文已经介绍了求太阳地心黄经的函数GetSunEclipticLongitudeECDegree()，本文也给出了求月亮地心黄经的函数GetMoonEclipticLongitudeECDegree()，因此可以定义方程为：

f(x) = GetMoonEclipticLongitudeECDegree(x) – GetSunEclipticLongitudeECDegree(x) = 0

其中x是儒略日单位的，我们要用牛顿迭代法求方程f(x)=0时的解x，也就是时间值。牛顿迭代法求解的迭代式是：

X_n+1 = X_n – f(X_n)/f’(X_n)

这里也不多解释了。导函数仍然使用近似公式，也不解释了，直接上迭代求解的代码了：

102 double CalculateMoonShuoJD(double tdJD)
103 {
104     double JD0, JD1,stDegree,stDegreep;
105 
106     JD1 = tdJD;
107     do
108     {
109         JD0 = JD1;
110         double moonLongitude = GetMoonEclipticLongitudeECDegree(JD0);
111         double sunLongitude = GetSunEclipticLongitudeECDegree(JD0);
112 
113         stDegree = moonLongitude - sunLongitude;
114 
115 
116         stDegreep = (GetMoonEclipticLongitudeECDegree(JD0 + 0.000005) - GetSunEclipticLongitudeECDegree(JD0 + 0.000005) - GetMoonEclipticLongitudeECDegree(JD0 - 0.000005) + GetSunEclipticLongitudeECDegree(JD0 - 0.000005)) / 0.00001;
117         JD1 = JD0 - stDegree / stDegreep;
118     }while((fabs(JD1 - JD0) > 0.00000001));
119 
120     return JD1;
121 }