中國(guó)農(nóng)歷的朔望月是農(nóng)歷歷法的基礎(chǔ)，而朔望月又是嚴(yán)格以日月合朔發(fā)生的那一天作為月首，因此日月合朔時(shí)間的計(jì)算是制定農(nóng)歷歷法的關(guān)鍵。本文將介紹ELP-2000/82月球運(yùn)行理論，以及如何用ELP-2000/82月球運(yùn)行理論計(jì)算日月合朔時(shí)間。

        要計(jì)算日月合朔時(shí)間，首先要對(duì)日月合朔這一天文現(xiàn)象進(jìn)行數(shù)學(xué)定義。朔望月是在地球上觀察到的月相周期，平均長(zhǎng)度約等于29.53059日，而恒星月（天文月）是月亮繞地球公轉(zhuǎn)一周的時(shí)間，長(zhǎng)度約27.32166日。月相周期長(zhǎng)度比恒星月長(zhǎng)大約兩天，這是因?yàn)樵谠虑蚶@地球旋轉(zhuǎn)一周的同時(shí)，地球還帶著它繞太陽(yáng)旋轉(zhuǎn)了一定的角度的緣故，所以月相周期不僅與月球運(yùn)行有關(guān)，還和太陽(yáng)運(yùn)行有關(guān)。日月合朔的時(shí)候，太陽(yáng)、月亮和地球三者接近一條直線，月亮未被照亮的一面對(duì)著地球，因此地球上看不到月亮，此時(shí)又被稱為新月。圖（1）就是日月合朔天文現(xiàn)象的示意圖：

圖（1）日月天文現(xiàn)象示意圖

月亮繞太陽(yáng)公轉(zhuǎn)的白道面和地球繞太陽(yáng)公轉(zhuǎn)的黃道面存在一個(gè)最大約5°的夾角，因此大多數(shù)情況下，日月合朔時(shí)都不是嚴(yán)格在同一條直線上，不過(guò)也會(huì)發(fā)生在同一直線的情況，此時(shí)就會(huì)發(fā)生日食。圖（1－b）顯示了日月合朔時(shí)側(cè)切面上月亮的三種可能的位置情況，當(dāng)月亮處在位置2時(shí)就會(huì)發(fā)生日食。由圖（1）可知，日月合朔的數(shù)學(xué)定義就是太陽(yáng)和月亮的地心視黃經(jīng)差為0的時(shí)刻。

        要計(jì)算日月合朔，需要知道太陽(yáng)地心視黃經(jīng)和月亮地心視黃經(jīng)的計(jì)算方法?！叭諝v生成算法”系列的第三篇《用天文方法計(jì)算二十四節(jié)氣》一文已經(jīng)介紹了如何用VSOP82/87行星理論計(jì)算太陽(yáng)的地心視黃經(jīng)，本文將繼續(xù)介紹如何用ELP-2000/82月球理論計(jì)算月亮的地心視黃經(jīng)。ELP-2000/82月球理論是M. Chapront-Touze和J. Chapront在1983年提出的一個(gè)月球位置的半解析理論，和其它半解析理論一樣，ELP-2000/82理論也包含一套計(jì)算方法和相應(yīng)的迭代周期項(xiàng)。這套理論共包含37862個(gè)周期項(xiàng)，其中20560個(gè)用于計(jì)算月球經(jīng)度，7684個(gè)用于計(jì)算月球緯度，9618個(gè)用于計(jì)算地月距離。但是這些周期項(xiàng)中有很多都是非常小的值，例如一些計(jì)算經(jīng)緯度的項(xiàng)對(duì)結(jié)果的增益只有0.00001角秒，還有一些地月距離周期項(xiàng)對(duì)距離結(jié)果的增益只有0.02米，對(duì)于精度不高的歷法計(jì)算，完全可以忽略。

        有很多基于ELP-2000/82月球理論的改進(jìn)或簡(jiǎn)化理論，《天文算法》一書的第四十五章就介紹了一種改進(jìn)算法，其周期項(xiàng)參數(shù)都是從ELP-2000/82理論的周期項(xiàng)參數(shù)轉(zhuǎn)換來(lái)的，忽略了小的周期項(xiàng)。使用該方法計(jì)算的月球黃經(jīng)精度只有10”，月亮黃緯精度只有4”，但是只用計(jì)算60個(gè)周期項(xiàng)，速度很快，本文就采用這種修改過(guò)的ELP-2000/82理論計(jì)算月亮的地心視黃經(jīng)。這種計(jì)算方法的周期項(xiàng)分三部分，分別用來(lái)計(jì)算月球黃經(jīng)，月球黃緯和地月距離，三部分的周期項(xiàng)的內(nèi)容一樣，由四個(gè)計(jì)算輻角的系數(shù)和一個(gè)正弦（或余弦）振幅組成。計(jì)算月球黃經(jīng)和月球黃緯使用正弦表達(dá)式求和：A * sin(θ)，計(jì)算地月距離用余弦表達(dá)式求和：A * cos(θ)，其中輻角θ的計(jì)算公式是：

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

4.1式中的四個(gè)輻角系數(shù)a、b、c和d由每個(gè)迭代周期項(xiàng)給出，日月距角D、太陽(yáng)平近地角M、月亮平近地角M’以及月球生交點(diǎn)平角距F則分別有4.2式－4.5式進(jìn)行計(jì)算：

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式)

以上各式計(jì)算結(jié)果的單位是度，其中T是儒略世紀(jì)數(shù)，T計(jì)算由4.6式計(jì)算：

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

以計(jì)算月球黃經(jīng)的周期項(xiàng)第二項(xiàng)的計(jì)算為例，第二項(xiàng)數(shù)據(jù)如下，輻角系數(shù)a = 2，b = 0，c = -1，d = 0，振幅A = 1274027，黃經(jīng)計(jì)算用正弦表達(dá)式，則I₂的計(jì)算如下所示：

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

在套用4.7式計(jì)算出60個(gè)月球黃經(jīng)周期項(xiàng)值的時(shí)候，需要注意對(duì)包含了太陽(yáng)平近地角M的項(xiàng)進(jìn)行修正，因?yàn)镸的值與地球公轉(zhuǎn)軌道的離心率有關(guān)，因?yàn)殡x心率是個(gè)與時(shí)間有關(guān)的變量，導(dǎo)致振幅A實(shí)際上是個(gè)變量，需要根據(jù)時(shí)間進(jìn)行修正。月球黃經(jīng)周期項(xiàng)的修正方法是：如果輻角中包含了M或-M時(shí)，需要乘以系數(shù)E修正；如果輻角中包含了2M或-2M，則需要乘以系數(shù)E的平方進(jìn)行修正。系數(shù)E的計(jì)算表達(dá)式如下：

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

其中T值由4.6式計(jì)算。上面的計(jì)算月球黃經(jīng)的第二個(gè)周期項(xiàng)中M對(duì)應(yīng)的系數(shù)是0，因此I₂不需要修正，但是第五個(gè)周期項(xiàng)中M對(duì)應(yīng)的系數(shù)是1，因此I₅需要乘以E進(jìn)行修正。套用4.7式計(jì)算出60個(gè)月球黃經(jīng)周期項(xiàng)值I₁－I₆₀之和ΣI：

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

        月球黃緯的周期項(xiàng)和Σb的計(jì)算方法與月球黃經(jīng)周期項(xiàng)和ΣI的計(jì)算方法一樣，每個(gè)月球黃緯周期項(xiàng)也包含振幅A和四個(gè)輻角系數(shù)a、b、c和d，對(duì)于太陽(yáng)平近地角M的系數(shù)b不是0的情況也需要乘以E或E²進(jìn)行修正。地月距離的周期項(xiàng)和Σr也可以按照上面的方法計(jì)算，計(jì)算地月距離的目的是為了計(jì)算月亮光行差，因?yàn)榈卦戮嚯x較小，從地球觀察月亮產(chǎn)生的光行差也很小，相對(duì)于本文算法的精度（月球黃經(jīng)精度10”，月亮黃緯精度4”）來(lái)說(shuō)，可以忽略光行差修正，因此就不用計(jì)算地月距離。

        由于金星和木星對(duì)月球的攝動(dòng)影響，需要對(duì)計(jì)算出的月球黃經(jīng)周期項(xiàng)和ΣI和月球黃緯周期項(xiàng)和Σb金星攝動(dòng)修正，修正的方法如下：

Σ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式計(jì)算得到，L’是月球平黃經(jīng)，計(jì)算方法是：

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

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

A1、A2和A4是攝動(dòng)角修正量，計(jì)算方法如下：

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式最終得到月亮的地心視黃經(jīng)和地心視黃緯：

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

β = Σb / 1000000.0                                                  （4.17式）

ΣI和Σb最后要除以1000000.0是因?yàn)橹芷陧?xiàng)系數(shù)中振幅A的單位是0.000001度，因此λ和β的單位是度。下面給出計(jì)算月球地心視黃經(jīng)的代碼：

函數(shù)參數(shù)dbJD是力學(xué)時(shí)儒略日時(shí)間，返回以度為單位的月球視黃經(jīng)。其中GetMoonEclipticParameter()函數(shù)分別根據(jù)4.2式、4.3式、4.4式、4.5式、4.8式和4.12式計(jì)算日月距角D、太陽(yáng)平近地角M、月亮平近地角M’、月球生交點(diǎn)平角距F、修正系數(shù)E和月球平黃經(jīng)L’，不需多說(shuō)明，只要根據(jù)以上各式直接計(jì)算即可。CalcMoonECLongitudePeriodicTbl()函數(shù)計(jì)算60個(gè)月球黃經(jīng)周期項(xiàng)和，并根據(jù)M值系數(shù)的情況進(jìn)行修正，算法實(shí)現(xiàn)如下：

 CalcMoonLongitudePerturbation()函數(shù)計(jì)算月球黃經(jīng)攝動(dòng)修正量，使用了4.13式和4.14式給出的計(jì)算方法：

        至此，本文已經(jīng)介紹了使用ELP-2000/82月球理論計(jì)算任意時(shí)刻月亮地心視黃經(jīng)的方法，結(jié)合“日歷生成算法”系列的第三篇《用天文方法計(jì)算二十四節(jié)氣》一文介紹的計(jì)算太陽(yáng)地心視黃經(jīng)的方法，就可以計(jì)算日月合朔的準(zhǔn)確時(shí)間了。由于ELP-2000/82月球理論也沒(méi)有根據(jù)月球黃經(jīng)反算時(shí)間的方法，因此本文也采用和《用天文方法計(jì)算二十四節(jié)氣》一文中一樣的牛頓迭代法計(jì)算日月合朔時(shí)間。

        關(guān)于牛頓迭代法可以參考相關(guān)的數(shù)學(xué)資料，“日歷生成算法”系列的第三篇《用天文方法計(jì)算二十四節(jié)氣》一文對(duì)如何使用牛頓迭代法有簡(jiǎn)單的介紹，可以參考一下?？偟膩?lái)說(shuō)，就是要先定義需要求解的方程f(x)，根據(jù)上文的介紹，我們需要求解的是太陽(yáng)的地心黃經(jīng)和月亮的地心黃經(jīng)差值是0的時(shí)候的時(shí)間，《用天文方法計(jì)算二十四節(jié)氣》一文已經(jīng)介紹了求太陽(yáng)地心黃經(jīng)的函數(shù)GetSunEclipticLongitudeECDegree()，本文也給出了求月亮地心黃經(jīng)的函數(shù)GetMoonEclipticLongitudeECDegree()，因此可以定義方程為：

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

其中x是儒略日單位的，我們要用牛頓迭代法求方程f(x)=0時(shí)的解x，也就是時(shí)間值。牛頓迭代法求解的迭代式是：

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

這里也不多解釋了。導(dǎo)函數(shù)仍然使用近似公式，也不解釋了，直接上迭代求解的代碼了：

        至本文結(jié)束，我們已經(jīng)能夠使用半解析算法計(jì)算太陽(yáng)的黃經(jīng)和月亮的黃經(jīng)，并且能夠通過(guò)牛頓迭代法或者24節(jié)氣的準(zhǔn)確時(shí)間和日月合朔的準(zhǔn)確時(shí)間，在這基礎(chǔ)上就可以進(jìn)行中國(guó)農(nóng)歷的推算了，“日歷生成算法”系列的下一篇將介紹中國(guó)農(nóng)歷的歷法規(guī)則和推算方法。

        再次說(shuō)明一下，以上算法中討論的時(shí)間都是力學(xué)時(shí)時(shí)間（TD），與國(guó)際協(xié)調(diào)時(shí)（UTC）以及各個(gè)時(shí)區(qū)的本地時(shí)間都有不同，以上計(jì)算出來(lái)的時(shí)間都需要調(diào)整成本地時(shí)間，比如中國(guó)的中原地區(qū)就是東八區(qū)標(biāo)準(zhǔn)時(shí)（UTC + 8）。應(yīng)用本文的算法計(jì)算出2012年前后15個(gè)日月合朔時(shí)間如下（已經(jīng)轉(zhuǎn)換為東八區(qū)標(biāo)準(zhǔn)時(shí)）：

2011-11-25, 14:09:41.25

2011-12-25, 02:06:27.25

2012-01-23, 15:39:24.16

2012-02-22, 06:34:40.84

2012-03-22, 22:37:08.91

2012-04-21, 15:18:22.12

2012-05-21, 07:46:59.97

2012-06-19, 23:02:06.39

2012-07-19, 12:24:02.83

2012-08-17, 23:54:28.03

2012-09-16, 10:10:36.99

2012-10-15, 20:02:30.98

2012-11-14, 06:08:05.90

2012-12-13, 16:41:37.60

2013-01-12, 03:43:31.34