【接上篇】

        經(jīng)過上述計算轉(zhuǎn)換得到坐標(biāo)值是理論值，或者說是天體的幾何位置，但是FK5系統(tǒng)是一個目視系統(tǒng)，也就是說體現(xiàn)的是人眼睛觀察效果（光學(xué)位置），這就需要根據(jù)地球的物理環(huán)境、大氣環(huán)境等信息做進(jìn)一步的修正，使其和人類從地球上觀察星體的觀測結(jié)果一致。

        首先需要進(jìn)行章動修正。章動是指地球沿自轉(zhuǎn)軸的指向繞黃道極緩慢旋轉(zhuǎn)過程中，由于地球上物質(zhì)分布不均勻性和月球及其它行星的攝動力造成的輕微抖動。英國天文學(xué)家詹姆斯·布拉德利(1693—1762)最早發(fā)現(xiàn)了章動，章動可以沿著黃道分解為水平分量和垂直分量，黃道上的水平分量記為Δψ，稱為黃經(jīng)章動，它影響了天球上所有天體的經(jīng)度。黃道上的垂直分量記為Δε，稱為交角章動，它影響了黃赤交角。目前編制天文年歷所依據(jù)的章動理論是伍拉德在1953年建立的，它是以剛體地球模型為基礎(chǔ)的。1977年，國際天文聯(lián)合會的一個專家小組建議采用非剛體地球模型――莫洛堅斯基II模型代替剛體地球模型計算章動，1979年的國際天文學(xué)聯(lián)合會第十七屆大會正式通過了這一建議，并決定于1984年正式實施。

        地球章動主要是月球運動引起的，也具有一定的周期性，可以描述為一些周期項的和，主要項的周期是6798.4日(18.6年)，但其它項是一些短周期項(小于10天)。本文采用的計算方法取自國際天文聯(lián)合會的IAU1980章動理論，周期項系數(shù)數(shù)據(jù)來源于《天文算法》一書第21章的表21－A，該表忽略了IAU1980章動理論中系數(shù)小于0.0003'的周期項，因此只有63項。每個周期項包括計算黃經(jīng)章動（Δψ）的正弦系數(shù)（相位內(nèi)項系數(shù)）、計算交角章動的（Δε）余弦系數(shù)（相位外項系數(shù)）以及計算輻角的5個基本角距（M、M'、D、F、Ω）的線性組合系數(shù)。5個基本角距的計算公式是：

平距角(日月對地心的角距離)：
D = 297.85036 455267.111480 * T - 0.0019142 * T² T³ / 189474        （3.10式）
太陽（地球）平近點角：
M = 357.52772 35999.050340 * T - 0.0001603 * T² - T³ / 300000         （3.11式）
月球平近點角
M'= 134.96298 477198.867398 * T 0.0086972 * T² T³ / 56250        （3.12式）

月球緯度參數(shù)：
F = 93.27191 483202.017538 * T - 0.0036825 * T² T³ / 327270          （3.13式）
黃道與月球平軌道升交點黃經(jīng)：
Ω= 125.04452 - 1934.136261 * T 0.0020708 * T² T³ / 450000            （3.14式）

以上各式中的T是儒略世紀(jì)數(shù)，計算出來的5個基本角距的單位都是度，在計算正弦或余弦時要轉(zhuǎn)換為弧度單位。計算每一個周期項的黃經(jīng)章動過程是這樣的，首先將3.10－3.14式計算出來的值與對應(yīng)的5個基本角距系數(shù)組合，計算出輻角。以本文使用的章動周期項系數(shù)表中的第七項為例，5個基本角距對應(yīng)的系數(shù)分別是1、0、-2、2和2，輻角θ的值就是：-2D M 2F 2Ω。計算出輻角后就可以計算周期項的值：

S = (S1 S2 * T) * sin(θ)                          （3.15式）

仍以第七項為例，S的值就是（-517 1.2 * T）* sin(θ)。對每一項的值S累加就可得到黃經(jīng)章動，單位是0.0001'。交角章動的計算方法與黃經(jīng)章動的計算類似，輻角θ的值是一樣的，只是計算章動使用的是余弦系數(shù)：

C = (C1 C2 * T) * cos(θ)                          （3.16式）

CalcEarthLongitudeNutation()函數(shù)就是計算黃經(jīng)章動的實現(xiàn)代碼：

 GetEarthNutationParameter()輔助函數(shù)用于計算5個基本角距：

同樣，計算交角章動的實現(xiàn)代碼是：

         除了章動修正，對于目測系統(tǒng)來說，還要進(jìn)行光行差修正。光行差是指在同一瞬間，運動中的觀察者所觀測到的天體視方向與靜止的觀測者所觀測到天體的真方向之差。造成光行差的原因有兩個，一個是光的有限速度，另一個是觀察者的運動。在地球上的天文觀測者因和地球一起運動（自傳＋公轉(zhuǎn)），他所看到的星光方向與假設(shè)地球不動時看到的方向不一樣。以太陽為例，光線從太陽傳到地球需要約8分鐘的時間，在這8分鐘多的時間中，地球沿著公轉(zhuǎn)軌道移動了一段距離人們根據(jù)現(xiàn)在的觀察認(rèn)定太陽在那個視位置，事實上那是8分鐘前太陽的位置。在精確的天文計算中，需要考慮這種光行差引起的視位置差異，在計算太陽的地心視黃經(jīng)時，要對其進(jìn)行光行差修正。地球上的觀測者可能會遇到幾種光行差，分別是因地球公轉(zhuǎn)引起的周年光行差，因地球自傳引起的周日光行差，還有因太陽系或銀河系運動形成的長期光行差等等，對于從地球上觀察太陽這種情況，只需要考慮周年光行差和周日光行差。因太陽公轉(zhuǎn)速度比較快，周年光行差最大可達(dá)到20.5角秒，在計算太陽視黃經(jīng)時需要考慮修正。地球自傳速度比較慢，周日光行差最大約為零點幾個角秒，因此計算太陽視黃經(jīng)時忽略周日光行差。

        下面是一個粗略計算太陽地心黃經(jīng)光行差修正量的公式，其中R是地球和太陽的距離：

AC = -20'.4898 / R                                    （3.17式）

分子20.4898并不是一個常數(shù)，但是其只的變化非常緩慢，在0年是20'.4893，在4000年是20'.4904。前文提到過，太陽到地球的距離R可以用VSOP87D表的R0－R5周期項計算出來，R的單位是“天文單位（AU）”，和計算太陽地心黃經(jīng)和地心黃緯類似，太陽到地球的距離可以這樣算出來：

也可以不使用VSOP，而用下面的公式直接計算日地距離R：

R = 1.000001018 (1 - e²) / (1 e * cos(v))                 （3.18式）

其中e是地球軌道的離心率：

e = 0.016708617 - 0.000042037 * T - 0.0000001236 * T²      （3.19式）

v的計算公式是v = M C，其中M是太陽平近地角：

M = 357.52910 35999.05030 * T - 0.0001559 * T² - 0.00000048 * T³     （3.20式）

中心C的太陽方程：

C = (1.914600 - 0.004817 * T - 0.000014 * T²) * sin(M)

(0.019993 - 0.000101 * T) * sin(2M)

0.000290 * sin(3M)                                              （3.21式）

以上各式中的T都是儒略世紀(jì)數(shù)，M和C的單位都是度，帶入3.18式計算時需要轉(zhuǎn)換成弧度單位，計算出R以后，就可以這樣計算光行差修正量：

AC = K / R (K是光行差常數(shù)，K = 20'.49552)          （3.22式）

無論是使用3.17式還是使用3.22式，最終計算出來的太陽光行差修正單位都是角秒。

        由VSOP87理論計算出來的幾何位置黃經(jīng)，經(jīng)過坐標(biāo)轉(zhuǎn)換，章動修正和光行差修正后，就可以得到比較準(zhǔn)確的太陽地心視黃經(jīng)，GetSunEclipticLongitudeEC()函數(shù)就是整個過程的代碼：

參數(shù)jde是力學(xué)時時間，單位是儒略日，返回太陽地心視黃經(jīng)，單位是度。

        到現(xiàn)在為止，我們已經(jīng)知道如何使用VSOP82/87理論計算以儒略日為單位的任意時刻的太陽地心視黃經(jīng)，但是這和實際歷法計算需求還不一致，歷法計算需要根據(jù)太陽地心視黃經(jīng)反求出此時的時間。VSOP82/87理論沒有提供反向計算的方法，但是可以采用根據(jù)時間正向計算太陽視黃經(jīng)，配合誤差修正進(jìn)行迭代計算的方法，使正向計算出來的結(jié)果向已知結(jié)果收斂，當(dāng)達(dá)到一定的迭代次數(shù)或計算結(jié)果與已知結(jié)果誤差滿足精度要求時，停止迭代，此時的正向輸入時間就是所求的時間。地球公轉(zhuǎn)軌道是近似橢圓軌道，軌道方程不具備單調(diào)性，但是在某個節(jié)氣附件的一小段時間區(qū)間中，軌道方程具有單調(diào)性，這個是本文迭代算法的基礎(chǔ)。

        實際上，我們要做的事情就是求解方程的根，但是我們面臨的這個方程沒有求根公式。對此類問題，數(shù)學(xué)上通常采用的迭代求解方法有二分逼近法和牛頓迭代法，事實上二分逼近法可以用更好的策略，比如用黃金分割代替二分法進(jìn)行逼近區(qū)間的選擇。接下來我們將分別介紹這兩種方法在計算二十四節(jié)氣中的應(yīng)用，首先介紹黃金分割逼近法。

        已知太陽視黃經(jīng)的值，反求對應(yīng)的時間的過程是這樣的，首先根據(jù)節(jié)氣對應(yīng)的視黃經(jīng)角度值W，估算出節(jié)氣可能的時間區(qū)間[A, B]，然后找到這個時間區(qū)間內(nèi)黃金分割點對應(yīng)的時間值C，C的計算采用3.23式：

C = ((B - A) * 0.618) A                               (3.23式)

用C值估算出太陽視黃經(jīng)W’，如果W’ > W，則調(diào)整調(diào)迭代時間區(qū)間為[A, C]，如果W’ < W，則調(diào)整迭代時間區(qū)間為[C, B]，然后重復(fù)上述過程，直到W’ 與W的差值滿足精度要求為止（區(qū)間上下限A和B的差值小于門限制也可以作為迭代退出條件）。采用黃金分割法進(jìn)行逼近求值的算法實現(xiàn)如下：

這里要特別說明一下，由于角度的360度圓周性，當(dāng)在太陽黃經(jīng)0度附近逼近時，區(qū)間的上下界可能分別位于（345, 360]和[0, 15)兩個區(qū)間上，此時需要將（345, 360]區(qū)間修正為（-15, 0]，使得逼近區(qū)間邊界的選取能夠正常進(jìn)行。EstimateSTtimeScope()函數(shù)估算節(jié)氣的時間區(qū)間，估算的依據(jù)是每個月的節(jié)氣時間比較固定，最多相差一兩天,考慮的幾千年后歲差的影響，這個估算范圍還可以再放寬一點，比如，對于月內(nèi)的第一個節(jié)氣，可以將時間范圍估算為4日到9日，對于月內(nèi)的第二個節(jié)氣，可以將時間范圍估算為16日到24日，保證迭代范圍內(nèi)有解。EstimateSTtimeScope()函數(shù)算法簡單，這里就不列出代碼了。

        二分逼近或黃金分割逼近算法實現(xiàn)簡單，很容易控制，但是也存在效率低，收斂速度慢的問題，現(xiàn)在我們介紹牛頓迭代法，牛頓迭代法是一種在實數(shù)域和復(fù)數(shù)域上近似求解方程的方法。假設(shè)我們要求解的方程是f(x) = 0，如果f(x)的導(dǎo)函數(shù)f’(x)是連續(xù)的，則在真實解x附近的區(qū)域內(nèi)任意一點x₀開始迭代，則牛頓迭代法必收斂，特別當(dāng)f’(x)不等于0的時候，牛頓迭代法是平方收斂的，也就是說，每迭代一次，結(jié)果的有效數(shù)字將增加一倍。

        簡單的說，對于方程f(x) = 0，f(x)的導(dǎo)函數(shù)是f’(x)，則牛頓迭代法的迭代公式是：

X_{n 1} = x_n – f(x_n)/f’(x_n)                              (3.24式)

現(xiàn)在問題就是如何確定方程f(x)。對于我們面臨的問題，可以理解為已知angle，通過GetSunEclipticLongitudeEC(solarTermsJD)函數(shù)反向求解solarTermsJD的值，因此我們的方程可以理解為：

f(x) = GetSunEclipticLongitudeEC(x) – angle = 0

確定了方程f(x)，剩下的問題就是求導(dǎo)函數(shù)f’(x)。嚴(yán)格的求解，應(yīng)該根據(jù)GetSunEclipticLongitudeEC()函數(shù)，以儒略千年數(shù)dt為自變量，按照函數(shù)求導(dǎo)的規(guī)則求出導(dǎo)函數(shù)。因為GetSunEclipticLongitudeEC()函數(shù)內(nèi)部是調(diào)用其他函數(shù)，因此可以理解為是一個多個函數(shù)組合的復(fù)合函數(shù)，類似f(x) = g(x) h(x, k(x)) p(x)這樣的形式，可以按照求導(dǎo)規(guī)則逐步對其求導(dǎo)得到導(dǎo)函數(shù)。但是我不打算這么做，因為我有更簡單的方法，那就是使用計算導(dǎo)數(shù)的近似公式。其實求導(dǎo)函數(shù)的目的就是為了得到某一點的導(dǎo)數(shù)，如果有近似公式可以直接得到這一點的導(dǎo)數(shù)，就不用費勁求導(dǎo)函數(shù)了。

        如果函數(shù)f(x)是單調(diào)函數(shù)，或者是在某個區(qū)間上是單調(diào)函數(shù)，則在此函數(shù)的其單調(diào)區(qū)間上某一點的導(dǎo)數(shù)值可以用近似公式計算，這個近似公式是：

f’(x₀) = (f(x₀ 0.000005) – f(x₀ – 0.000005)) / 0.00001            (3.25式)

這是一個精度很高的近似公式，完全可以滿足民用歷法計算的精度要求。

        根據(jù)以上分析結(jié)果，使用牛頓迭代法求解節(jié)氣的算法就很容易實現(xiàn)了，以下就是牛頓迭代法求解節(jié)氣的代碼：

經(jīng)過驗證，牛頓迭代法具有非常好的收斂效果，一般只需3次迭代就可以得到滿足精度的結(jié)果。

        至此，我們就有了完整的計算節(jié)氣發(fā)生時間的方法，輸入年份和節(jié)氣對應(yīng)的太陽黃經(jīng)度數(shù)，即可求的節(jié)氣發(fā)生的精確時間。最后說明一下，以上算法中討論的時間都是力學(xué)時時間（TD），與國際協(xié)調(diào)時（UTC）以及各個時區(qū)的本地時間都有不同，以上計算出來的時間都需要調(diào)整成本地時間，比如中國的中原地區(qū)就是東八區(qū)標(biāo)準(zhǔn)時（UTC 8）。關(guān)于力學(xué)時、國際協(xié)調(diào)時（世界時）的定義，請參考文末的小知識3：力學(xué)時、原子時和國際協(xié)調(diào)時。應(yīng)用本文的算法計算出2012年各個節(jié)氣的時間如下（已經(jīng)轉(zhuǎn)換為東八區(qū)標(biāo)準(zhǔn)時），與紫金山天文臺發(fā)布的《2012中國天文年歷》中發(fā)布的時間在分鐘級別上完全吻合（此年歷只精確到分鐘）：

2012-01-06, 06:43:54.28   小寒

2012-01-21, 00:09:49.08   大寒

2012-02-04, 18:22:22.53   立春

2012-02-19, 14:17:35.37   雨水

2012-03-05, 12:21:01.56   驚蟄

2012-03-20, 13:14:24.17   春分

2012-04-04, 17:05:34.65   清明

2012-04-20, 00:12:03.28   谷雨

2012-05-05, 10:19:39.54   立夏

2012-05-20, 23:15:30.28   小滿

2012-06-05, 14:25:52.96   芒種

2012-06-21, 07:08:46.98   夏至

2012-07-07, 00:40:42.66   小暑

2012-07-22, 18:00:50.72   大暑

2012-08-07, 10:30:31.88   立秋

2012-08-23, 01:06:48.41   處暑

2012-09-07, 13:28:59.41   白露

2012-09-22, 22:48:57.14   秋分

2012-10-08, 05:11:41.45   寒露

2012-10-23, 08:13:32.83   霜降

2012-11-07, 08:25:56.47   立冬

2012-11-22, 05:50:08.09   小雪

2012-12-07, 01:18:55.23   大雪

2012-12-21, 19:11:35.61   冬至

小知識3：力學(xué)時、原子時和國際協(xié)調(diào)時

        力學(xué)時全稱是“牛頓力學(xué)時”，也被稱作是“歷書時”。它描述天體運動的動力學(xué)方程中作為時間自變量所體現(xiàn)的時間，或天體歷表中應(yīng)用的時間，是由天體力學(xué)的定律確定的均勻時間。力學(xué)時的初始?xì)v元取為1900年初附近，太陽幾何平黃經(jīng)為279°41′48″.04的瞬間，秒長定義為1900.0年回歸年長度的1／31556925.9747。1958年國際天文學(xué)聯(lián)合會決議決定：自1960年開始用力學(xué)時代替世界時作為基本的時間計量系統(tǒng)，規(guī)定天文年歷中太陽系天體的位置都按力學(xué)時推算。力學(xué)時與世界時之差由觀測太陽系天體（主要是月球）定出，因此力學(xué)時的測定精度較低，1967年起被原子時代替作為基本時間計量系統(tǒng)。

        國際協(xié)調(diào)時又稱世界時，是以本初子午線的平子夜起算的平太陽時，又稱格林威治時間。世界各地地方時與世界時之差等于該地的地理經(jīng)度。世界時1960年以前曾作為基本時間計量系統(tǒng)被廣泛應(yīng)用。由于地球自轉(zhuǎn)速度變化的影響，它不是一種均勻的時間系統(tǒng)。后來世界時先后被歷書時和原子時所取代。

        原子時是以物質(zhì)的原子內(nèi)部發(fā)射的電磁振蕩頻率為基準(zhǔn)的時間計量系統(tǒng)。原子時的初始?xì)v元規(guī)定為1958年1月1日世界時0時，秒長定義為銫-133原子基態(tài)的兩個超精細(xì)能級間在零磁場下躍遷輻射9192631770周所持續(xù)的時間。這是一種均勻的時間計量系統(tǒng)。1967年起，原子時已取代力學(xué)時作為基本時間計量系統(tǒng)。

參考文章：

[1] 《Secular variations of the planetary orbits》http://www./ma/enwiki/en/Secular_variations_of_the_planetary_orbits

[2] Jean.Meeus.Astronomical.Algorithms(天文算法)

[3] M.Chapront-Touze and J.Chapront.ELP 2000-85 - A semi-analytical lunar ephemeris adequate for historical times.Astronomy And Astrophysics,1998

[4] P.Bretagnon and G.Francou.Planetray theories in rectangular and spherical variables VSOP87 solutions. Astronomy And Astrophysics,1998