算法系列之十八:用天文方法計算二十四節氣(上)
二十四節氣在中國古代曆法中扮演著非常重要的角色,本文將介紹二十四節氣的基本知識,以及如何使用VSOP82/87行星執行理論計算二十四節氣發生的準確時間。
中國古代曆法都是以月亮執行規律為主,嚴格按照朔望月長度定義月,但是由於朔望月長度和地球迴歸年長度無法協調,會導致農曆季節和天氣的實際冷暖無法對應,因此聰明的古人將月亮執行規律和太陽執行規律相結合制定了中國農曆的歷法規則。在這種特殊的陰陽結合的歷法規則中,二十四節氣就扮演著非常重要的作用,它是聯絡月亮執行規律和太陽執行規律的紐帶。正是由於二十四節氣結合置閏規則,使得農曆的春夏秋冬四季和地球繞太陽運動引起的天氣冷暖變化相一致,成為中國幾千年來生產、生活的依據。
二十四節氣起源於中國黃河流域。遠在春秋時代,古人就開始使用仲春、仲夏、仲秋和仲冬四個節氣指導農耕種植。後來經過不斷地改進與完善,到秦漢年間,二十四節氣已經基本確立。公元前 104年,漢武帝頒佈由鄧平等人制定的《太初曆》,正式把二十四節氣訂於曆法,明確了二十四節氣的天文位置。二十四節氣天文位置的定義,就是從太陽黃經零度開始,沿黃經每執行15度所經歷的時日稱為“一個節氣”。太陽一個迴歸年執行360度,共經歷24個節氣,每個公曆月對應2個節氣。其中,每月第一個節氣為“節令”,即:立春、驚蟄、清明、立夏、芒種、小暑、立秋、白露、寒露、立冬、大雪和小寒等12個節令;每月的第二個節氣為“中氣”,即:雨水、春分、穀雨、小滿、夏至、大暑、處暑、秋分、霜降、小雪、冬至和大寒等12
為了更好地理解二十四節氣的天文位置,首先要解釋幾個天文學概念。“天球”是人們為了研究天體的位置和運動規律而引入的一個假象的球體,根據觀察點(也就是球心)的位置不同,可分為“日心天球”、“地心天球”等等。圖(1)就是天球概念的一個簡單示意圖:
圖(1)天球概念示意圖
天文學中常用的一個座標體系就是“地心天球”,它與地球同心且有相同的自傳軸,理論上具有無限大的半徑。地球的赤道和南北極點延伸到天球上,對應著天赤道和南北天極點。和地球上用經緯度定為位置一樣,天球也劃分了經緯度,分別命名為“赤經”和“赤緯”,地球上的經度用的是度(分秒)為單位,赤經以時(分秒)為單位。天空中的所有天體都可以投射到天球上,用赤經和赤緯定為天體在天球上的位置。“黃道(Ecliptic)”是地球繞太陽公轉軌道的軌道平面與天球(地心天球)相交的大圓,由於地球公轉受月球和其它行星的攝動,地球的公轉軌道並不是嚴格的平面,因此黃道的嚴格定義是:地月系質心繞太陽公轉的瞬時平均軌道平面與天球相交的大圓。黃道和天赤道所在的兩個平面並不是重疊的,它們之間存在一個23
黃道平面可以近似理解為地球繞太陽公轉的平面,以黃道為基圈的黃道座標系根據觀測中心是太陽還是地球還可以區分為日心座標系和地心座標系,對應天體的黃道座標分別被稱為“日心黃經、日心黃緯”和“地心黃經、地心黃緯”。日心黃經和日心黃緯比較容易理解,因為太陽系的行星都是繞太陽公轉,以太陽為中心將這些行星向天球上投影是最簡單的確定行星位置關係的做法。但是人類自古觀察太陽的週年運動,都是以地球為參照,以太陽的週年視運動位置來計算太陽的執行軌跡,使用的其實都是地心黃經和地心黃緯,要了解古代曆法,理解這一點非常重要。圖(2)就解釋了造成這種視覺錯覺的原因:
圖(2)太陽黃道視覺位置原理圖
古人定義二十四節氣的位置,是太陽沿著黃道執行時的視覺位置,每個節氣對應的黃道經度其實是地心黃經。從圖(2)可以看出日心黃經和地心黃經存在180度的轉換關係,同樣可以理解,日心黃緯和地心黃緯在方向上是反的,因此可以很方便地將兩類座標相互轉換,轉換公式是:
太陽地心黃經 = 地球日心黃經 + 180° (3.1式)
太陽地心黃緯 = -地球日心黃緯 (3.2式)
瞭解了以上的天文學基礎之後,就可以著手對二十四節氣的發生時間進行計算。我們常說的節氣發生時間,其實就是在太陽沿著黃道做視覺運動過程中,當太陽地心黃經等於某個節氣黃經度數時的那個瞬間的時間。所謂的用天文演算法計算二十四節氣時間,就是根據牛頓力學原理或開普勒三大行星定律,計算出與曆法密切相關的地球、太陽和月亮三個天體的執行軌道和時間引數,以此得出當這些天體位於某個位置時的時間。這樣的天文計算需要計算者有紮實的微積分學、幾何學和球面三角學知識,令廣大天文愛好者望而卻步。但是隨著VSOP-82/87行星理論以及ELP-2000/82月球理論的出現,使得天文計算變得簡單易行,本文就是以VSOP-82/87行星理論為計算依據,計算二十四節氣的準確時間。
古代天文學家在對包括地球和月亮在內的行星執行軌道精確計算後發現,天體的執行因為受相近天體的影響,並不嚴格遵循理論方法計算出來的軌道,而是在理論軌道附近波動。這種影響在天文學上被稱為攝動,攝動很難被精確計算,只能根據經驗估算。但是經過長期的觀測和計算,天文學家發現行星軌道因為攝動影響而產生的波動其實也是有規律的,即在相當長的時間內呈現出週期變化的趨勢。於是天文學家開始研究這種週期變化,希望通過一種類似曲線擬合的方法,對一些週期計算項按照某種計算式迭代求和計算代替積分計算來模擬行星執行軌跡。這種計算式可以描述為:a + bt + ct2 + … xcos(p + qt + rt2 + …),其中t是時間引數,這樣的理論通常被稱為半解析(semi-analytic)理論。其實早在十八世紀,歐洲學者Joseph Louis Lagrange就開始嘗試用這種週期項計算的方法修正行星軌道,但是他採用的週期項計算式是線性方程,精度不高。
1982年,P.Bretagnon公開發表了VSOP行星理論(這個理論的英文名稱是:Secular Variations of the Planetary Orbits,VSOP的縮寫其實是源於法文名稱:Variations Séculaires des Orbites Planétaires),VSOP理論是一個描述太陽系行星軌道在相當長時間範圍內週期變化的半分析(semi-analytic)理論。VSOP82理論是VSOP理論的第一個版本,提供了對太陽系幾大行星位置計算的週期序列,通過對週期序列進行正弦或餘弦項累加求和,就可以得到這個行星在給定時間的軌道引數。不過VSOP82由於每次都會計算出全部超高精度的軌道引數,這些軌道引數對於曆法計算這樣的民用場合很不適用。1987年,Bretagnon 和 Francou 建立了VSOP87行星理論,VSOP87行星理論不僅能計算各種精密的軌道引數,還可以直接計算出行星的位置,行星位置可以是各種座標系,包括黃道座標系。VSOP87行星理論由6張週期項係數表組成,分別是VSOP87、VSOP87A、VSOP87B、VSOP87C、VSOP87D和VSOP87E,其中VSOP87D表可以直接計算行星日心黃經(L)、日心黃緯(B)和到太陽的距離(R),此表計算出的結果適用於節氣位置判斷。
VSOP87D表包含了三部分資料,分別是計算行星日心黃經的週期項係數表(L表)、計算行星日心黃緯的週期項係數表(B表)和計算行星和太陽距離的週期項係數表(R表)。VSOP87D表有太陽系8大行星的資料,本文的計算只關心與地球相關的資料。L表由L0-L5六部分組成,每一部分都包含若干個週期項係數條目,比如L0表有559個週期項係數條目,L1表有341個條目等等。L表的每個週期項係數條目包含若干個引數,用於計算各種軌道引數和位置引數,計算地球的日心黃經只需要用到其中三個係數。計算所有的週期項係數並不是必須的,有時候減少一些係數比較小的週期項可以減少計算所花費的時間,當然,這會犧牲一點精度。假設計算地球日心黃經的三個係數是A、B和C,則每個週期項的計算表示式是:
A * cos(B + Cτ) (3.3式)
其中τ是儒略千年數,τ的計算公式如下:
τ = (JDE - 2451545.0) / 365250 (3.4式)
JDE是計算軌道引數的時間,單位是儒略日,2451545.0是公元2000年1月1日 12時的儒略日數,關於儒略日的概念,請參考“日曆生成演算法”的第一篇《中國公曆(格里曆)》中的說明以及計算方法。以L0表的第二個週期項為例,這個週期項資料中與日心黃經計算有關的三個係數分別是A= 3341656.456,B=4.66925680417,C=6283.07584999140,則第二個週期項的計算方法是:3341656.456 * cos(4.66925680417 + 6283.0758499914 * τ)。對L0表的各項分別計算後求和可得到L0表週期項總和L0,對L表的其它幾個部分使用相同的方法計算週期項和,可以得到L1、L2、L3、L4和L5,然後用用3.5式計算出最終的地球日心黃經,單位是弧度:
L = (L0 + L1 * τ+ L2 * τ2 + L3 * τ3 + L4 * τ4 +L5 * τ5) / 108 (3.5式)
用同樣的方法對地球日心黃緯的週期項係數表和計算行星和太陽距離的週期項係數表計算求和,可以得到地球日心黃緯B和日地距離R,B的單位是弧度,R的單位是天文單位(AU)[1]。由於3.5式的計算方法需要多次計算τ的乘方,浮點數的乘方計算的速度比較慢,實際計算時,通常對3.5式進行變換,用乘法和加法代替直接的乘方計算,這是一種常用的轉換:
L = (((((L5 * τ + L4) * τ + L3) * τ + L2) * τ + L1) * τ + L0) / 108 (3.6式)
本文就是使用3.6式代替3.5式進行計算。
VSOP82/87行星理論中的週期項係數對不同的行星具有不同的精度,對地球來說,在1900-2100年之間的200年跨度期間,計算精度是0.005"。前文曾說過,對於不需要這麼高精度的計算應用時,可以適當減少一些係數比較小的週期項,減少計算量,提高計算速度。Jean Meeus在他的《天文演算法》一書中就給出了一套精簡後的VSOP87D表的週期項,將計算地球黃經的L0表由原來的559項精簡到64項,計算地球黃緯的B0表甚至被精簡到只有5項,從實際效果看,計算精度下降並不多,但是極大地減少了計算量。
使用VSOP87D週期項係數表計算得到的是J2000.0平黃道和平春分點(mean dynamic ecliptic and equinox)為基準的日心黃經(L)和日心黃緯(B),其值與標準FK5系統略有差別,如果對精度要求很高可以採用下面的方法將計算得到的日心黃經(L)和日心黃緯(B)轉到FK5系統[2]:
首先然後 L',單位是度:
L' = L - 1.397 * T - 0.00031*T2 (3.7式)
3.7式中的T是儒略世紀數,它與儒略千年數τ的關係是:T = 10 *τ。然後使用L'計算L和B的修正值ΔL和ΔB:
ΔL = -0.09033 + 0.03916 * ( cos(L') + sin(L') ) * tan(B) (3.8式)
ΔB = +0.03916 * ( cos(L') - sin(L') ) (3.9式)
ΔL和ΔB的單位都是",是度分秒角度單位體系,需要將3.6式計算出得L和B轉換成以度(°)為單位的值後再進行修正。
CalcSunEclipticLongitudeEC()函式就是使用VSOP87行星理論計算行星日心黃經的程式碼實現,整個計算過程和前文描述一樣,首先根據VSOP87D表的資料計算出L0-L5,然後用3.6式計算出地球的日心黃經,3.6式計算出來的單位是弧度,因此轉換成度分秒單位,最後使用3.1式將結果轉換成太陽的地心黃經:
double CalcSunEclipticLongitudeEC(double dt) { double L0 = CalcPeriodicTerm(Earth_L0, sizeof(Earth_L0) / sizeof(VSOP87_COEFFICIENT), dt); double L1 = CalcPeriodicTerm(Earth_L1, sizeof(Earth_L1) / sizeof(VSOP87_COEFFICIENT), dt); double L2 = CalcPeriodicTerm(Earth_L2, sizeof(Earth_L2) / sizeof(VSOP87_COEFFICIENT), dt); double L3 = CalcPeriodicTerm(Earth_L3, sizeof(Earth_L3) / sizeof(VSOP87_COEFFICIENT), dt); double L4 = CalcPeriodicTerm(Earth_L4, sizeof(Earth_L4) / sizeof(VSOP87_COEFFICIENT), dt); double L5 = CalcPeriodicTerm(Earth_L5, sizeof(Earth_L5) / sizeof(VSOP87_COEFFICIENT), dt); double L = (((((L5 * dt + L4) * dt + L3) * dt + L2) * dt + L1) * dt + L0) / 100000000.0; /*地心黃經 = 日心黃經 + 180度*/ return (Mod360Degree(Mod360Degree(L / RADIAN_PER_ANGLE) + 180.0)); } |
Mod360Degree()函式將大於360°或小於0°的值調整到0-360°之間,便於轉換顯示。CalcPeriodicTerm()函式使用3.3式對一個週期項係數表進行求和計算,可以指定需要計算的週期項數:
double CalcPeriodicTerm(const VSOP87_COEFFICIENT *coff, int count, double dt) { double val = 0.0; for(int i = 0; i < count; i++) val += (coff[i].A * cos((coff[i].B + coff[i].C * dt))); return val; } |
同樣的方法計算太陽的地心黃緯:
double CalcSunEclipticLatitudeEC(double dt) { double B0 = CalcPeriodicTerm(Earth_B0, sizeof(Earth_B0) / sizeof(VSOP87_COEFFICIENT), dt); double B1 = CalcPeriodicTerm(Earth_B1, sizeof(Earth_B1) / sizeof(VSOP87_COEFFICIENT), dt); double B2 = CalcPeriodicTerm(Earth_B2, sizeof(Earth_B2) / sizeof(VSOP87_COEFFICIENT), dt); double B3 = CalcPeriodicTerm(Earth_B3, sizeof(Earth_B3) / sizeof(VSOP87_COEFFICIENT), dt); double B4 = CalcPeriodicTerm(Earth_B4, sizeof(Earth_B4) / sizeof(VSOP87_COEFFICIENT), dt); double B = (((((B4 * dt) + B3) * dt + B2) * dt + B1) * dt + B0) / 100000000.0; /*地心黃緯 = -日心黃緯*/ return -(B / RADIAN_PER_ANGLE); } |
AdjustSunEclipticLongitudeEC()函式根據3.8式計算黃經的修正量,longitude和latitude引數是由VSOP87理論計算出的太陽地心黃經和地心黃緯,單位是度,dt是儒略千年數,返回值單位是度:
double AdjustSunEclipticLongitudeEC(double dt, double longitude, double latitude) { double T = dt * 10; //T是儒略世紀數 double dbLdash = longitude - 1.397 * T - 0.00031 * T * T; // 轉換為弧度 dbLdash *= RADIAN_PER_ANGLE; return (-0.09033 + 0.03916 * (cos(dbLdash) + sin(dbLdash)) * tan(latitude * RADIAN_PER_ANGLE)) / 3600.0; } |
經過上述計算轉換得到座標值是理論值,或者說是天體的幾何位置,但是FK5系統是一個目視系統,也就是說體現的是人眼睛觀察效果(光學位置),這就需要根據地球的物理環境、大氣環境等資訊做進一步的修正,使其和人類從地球上觀察星體的觀測結果一致。
【下篇將介紹修正理論和修正演算法,請繼續關注】
小知識1:天文單位
天文單位(英文:Astronomical Unit,簡寫AU)是一個長度的單位,約等於地球跟太陽的平均距離。天文單位是天文常數之一,是天文學中測量距離,特別是測量太陽系內天體之間的距離的基本單位。地球到太陽的平均距離大約為一個天文單位,約等於1.496億千米。 1976年,國際天文學聯會把一天文單位定義為一顆質量可忽略、公轉軌道不受干擾而且公轉週期為365.2568983日(即一高斯年)的粒子與一個質量相等約一個太陽的物體的距離。當前普遍被接受並使用的天文單位的值是149,597,870,691±30米(約一億五千萬公里)。
小知識2:FK5系統
FK5常用的目視星表系統,又稱第五基本星表,是在FK4(第四基本星表)的基礎上發展出來的,對FK4星表進行了修正,於1984年正式啟用。它定義了一個以太陽質心為中心,J2000.0平赤道和春分點為基準的天球平赤道座標系。近年來國際上又編制了FK6星表(第六基本星表),但是還沒有被正式啟用。