日歷在我們的生活中扮演著十分重要的角色，上班、上學、約會都離不開日歷。每年新年開始，人們都要更換新的日歷，你想知道未來一年的這么多天是怎么被確定下來的嗎？為什么去年的國慶節(jié)是星期五而今年的國慶節(jié)是星期三？那就來研究一下日歷算法吧。本文將介紹日歷的編排規(guī)則，確定某日是星期幾的計算方法，以及如何在計算機上打印某一年的年歷。

要研究日歷算法，首先要知道日歷的編排規(guī)則，也就是歷法。所謂歷法，指的就是推算年、月、日的時間長度和它們之間的關系，指定時間序列的法則。我國的官方歷法是中國公歷，也就是世界通用的格里歷（Gregorian Calendar），中國公歷的年分為平常年和閏年，平常年一年是365天，閏年一年是366天。判定一年是平常年還是閏年的規(guī)則如下：

1、如果年份是 4 的倍數(shù)，且不是 100 的倍數(shù)，則是閏年；

2、如果年份是 400 的倍數(shù)，則是閏年；

3、不滿足 1、2 條件的就是平常年。

總結成一句話就是：四年一閏，百年不閏，四百年再閏。

中國公歷關于月的規(guī)則是這樣的，一年分為十二個月，其中一月、三月、五月、七月、八月、十月和十二月是大月，一個月有 31 天。四月、六月、九月和十一月是小月，一個月有 30 天。二月天數(shù)要根據(jù)是否是閏年來定，如果是閏年，二月是 29 天，如果是平常年，二月是 28 天。

2、計算星期

除了年月日，人們日常生活中還對日期定義了另一個屬性，就是星期幾。星期并不是公歷范疇內的東西，但是人們已經(jīng)習慣用星期來管理和規(guī)劃時間，比如一個星期工作五天，休息兩天等等，星期的規(guī)則徹底改變了人們的生活習慣，因此星期已經(jīng)成為歷法中的一部分了。星期的命名最早起源于古巴比倫文化。公元前 7－6 世紀，巴比倫人就使用了星期制，一個星期中的每一天都有一個天神掌管。這一規(guī)則后來傳到古羅馬，并逐漸演變成現(xiàn)在的星期制度。

如何知道某一天到底是星期幾？除了查日歷之外，是否有辦法推算出來某一天是星期幾呢？答案是肯定的，星期不象年和月那樣有固定的歷法規(guī)則，但是星期的計算也有自己的規(guī)律。星期是固定的 7 天周期，其排列順序固定，不隨閏年、平常年以及大小月的天數(shù)變化影響。因此，只要確切地知道某一天是星期幾，就可以推算出其它日期是星期幾。推算的方法很簡單，就是計算兩個日期之間相差多少天，用相差的天數(shù)對 7 取余數(shù)，這個余數(shù)就是兩個日期的星期數(shù)的差值。舉個例子，假設已經(jīng)知道 1977 年 3 月 27 日是星期日，如何得知 1978 年 3 月 27 日是星期幾？按照前面的方法，計算出 1977 年 3 月 27 日到 1978 年 3 月 27 日之間相差 365 天，365 除以 7 余數(shù)是 1，所以 1978 年 3 月 27 日就是星期一。

上述方法計算星期幾的關鍵是求出兩個日期之間相隔的天數(shù)。有兩種常用的方法計算兩個日期之間相隔的天數(shù)，一種是利用公歷的月和年的規(guī)則直接計算，另一種是利用儒略日計算。利用公歷規(guī)則直接計算兩個日期之間相差的天數(shù)，簡單地講就是將兩個日期之間相隔的天數(shù)分成三個部分：前一個日期所在年份還剩下的天數(shù)、兩個日期之間相隔的整數(shù)年所包含的天數(shù)和后一個日期所在的年過去的天數(shù)。如果兩個日期是相鄰兩個年份的日期，則第二部分整年的天數(shù)就是 0。以 1977 年 3 月 27 日到 2005 年 5 月 31 日為例，1977 年還剩下的天數(shù)是 279 天，中間整數(shù)年是從 1978 年到 2005 年（不包括2005 年），共 26 年，包括 7 個閏年和 20 個平常年，總計 9862 天，最后是 2005 年從 1 月 1 日到 5 月 31 日經(jīng)過的天數(shù) 151 天。三者總和 10292 天。直接利用公歷規(guī)則計算日期相差天數(shù)的算法實現(xiàn)如下（為了簡化算法復雜度，這個實現(xiàn)假設用于定位星期的那個日期總是在需要計算星期幾的那個日期之前）：

int CalculateDays(int ys, int ms, int ds, int ye, int me, int de)

{

int days = CalcYearRestDays(ys, ms, ds);

if(ys != ye) /*不是同一年的日期*/

{

if((ye - ys) >= 2) /*間隔超過一年，要計算間隔的整年時間*/

{

days += CalcYearsDays(ys + 1, ye);

}

days += CalcYearPassedDays(ye, me, de);

}

else

{

days = days - CalcYearRestDays(ye, me, de);

}

return days;

}

/*計算一年中過去的天數(shù)，包括指定的這一天*/

int CalcYearPassedDays(int year, int month, int day)

{

int passedDays = 0;

int i;

for(i = 0; i < month - 1; i++)

{

passedDays += daysOfMonth[i];

}

passedDays += day;

if((month > 2) && IsLeapYear(year))

passedDays++;

return passedDays;

}

/*計算一年中還剩下的天數(shù)，不包括指定的這一天*/

int CalcYearRestDays(int year, int month, int day)

{

int leftDays = daysOfMonth[month - 1] - day;

int i;

for(i = month; i < MONTHES_FOR_YEAR; i++)

{

leftDays += daysOfMonth[i];

}

if((month <= 2) && IsLeapYear(year))

leftDays++;

return leftDays;

}

計算years年1月1日和yeare年1月1日之間的天數(shù)，

包括years年1月1日，但是不包括yeare年1月1日

int CalcYearsDays(int years, int yeare)

{

int days = 0;

int i;

for(i = years; i < yeare; i++)

{

if(IsLeapYear(i))

days += DAYS_OF_LEAP_YEAR;

else

days += DAYS_OF_NORMAL_YEAR;

}

return days;

}

另一種計算兩個日期相差天數(shù)的方法是利用儒略日（Julian Day，JD）進行計算。首先介紹一下儒略日，儒略日是一種不記年，不記月，只記日的歷法，是由法國學者 Joseph Justus Scaliger（1540－1609）在 1583 年提出來的一種以天數(shù)為計量單位的流水日歷。儒略日和儒略歷（Julian Calendar）沒有任何關系，命名為儒略日也僅僅他本人為了紀念他的父親--意大利學者 Julius Caesar Scaliger（1484－1558）。簡單來講，儒略日就是指從公元前 4713 年 1 月 1 日 UTC 12:00 開始所經(jīng)過的天數(shù)，JD0 就被指定為公元前 4713 年 1 月 1 日 12:00 到公元前 4713 年 1 月 2 日 12:00 之間的 24 小時，依次順推，每一天都被賦予一個唯一的數(shù)字。例如從 1996 年 1 月 1 日 12:00 開始的一天就是儒略日JD2450084。使用儒略日可以把不同歷法的年表統(tǒng)一起來，很方便地在各種歷法中追溯日期。如果計算兩個日期之間的天數(shù)，利用儒略日計算也很方便，先計算出兩個日期的儒略日數(shù)，然后直接相減就可以得到兩個日期相隔的天數(shù)。

由公歷的日期計算出儒略日數(shù)是一個很簡單的事情，有多個公式可以計算儒略日，本文選擇如下公式計算儒略日：

其中 y 是年份，m 是月份，d 是日期，如果 m 小于或等于 2，則 m 修正為 m+12，同時年份修正為 y-1。c 值由以下方法計算：

下面就是由公歷日期計算儒略日的算法實現(xiàn)：

119 int CalculateJulianDay(int year, int month, int day)

120 {

121 int B = 0;

122

123 if(month <= 2)

124 {

125 month += 12;

126 year -= 1;

127 }

128 if(IsGregorianDays(year, month, day))

129 {

130 B = year / 100;

131 B = 2 - B + year / 400;

132 }

133

134 double dd = day + 0.5000115740; /*本日12:00后才是儒略日的開始(過一秒鐘)*/

135 return int(365.25 * (year + 4716) + 0.01) + int(30.60001 * (month + 1)) + dd+ B - 1524.5;

136 }

儒略日的計算通常精確到秒，得到的 JD 數(shù)也是一個浮點數(shù)，本文僅僅是為了計算日期相隔的整數(shù)天數(shù)，因此都采用整數(shù)計算。由于儒略日的周期開始與每天中午12:00，而歷法中的天數(shù)通常是從0:00開始的，因此儒略日計算上對日期的天數(shù)進行了修正。1977年3月27日的儒略日是2443230，2005年5月31日的儒略日是2453522，差值是10292，和前一種方法計算的結果一致。

我們用兩種方法計算出兩個日期之間的天數(shù)都是10292，現(xiàn)在用10292除以7得到余數(shù)是2，也就是說2005年5月31日與1977年3月27日星期數(shù)差兩天，所以2005年5月31日就是是星期二。

上述計算星期的方法雖然步驟簡單，但是每次都要計算兩個日期的時間差，不是非常方便。如果能夠有一個公式可以直接根據(jù)日期計算出對應的星期豈不是更好？幸運的是，這樣的公式是存在的。此類公式的推導原理仍然是通過兩個日期的時間差來計算星期，只是通過選擇一個特殊的日期來簡化公式的推導。這個所謂的特殊日期指的是某一年的 12 月 31 日這天剛好是星期日這種情況。選擇這樣的日子有兩個好處，一個是計算上可以省去計算標準日期這一年的剩余天數(shù)，另一個是計算出來的日期差余數(shù)是幾就是星期幾，不需要再計算星期的差值。人們知道公元元年的 1 月 1 日是星期一，那么公元前 1 年的 12 月 31 日就是星期日，用這一天作為標準日期，就可以只計算整數(shù)年的時間和日期所在的年積累的天數(shù)，這個星期公式就是：

w = (L * 366 + N * 365 + D) % 7 (公式 2)

公式中的 L 是從公元元年到 y 年 m 月 d 日所在的年之間的閏年次數(shù)，N 是平常年次數(shù)，D 是 y 年內的積累天數(shù)。將整年數(shù) y - 1 = L + N 帶入上式，可得：

w = ( (y - 1) * 365 + L + D) % 7 (公式 3)

根據(jù)閏年規(guī)律，從公元元年到y(tǒng)年之間的閏年次數(shù)是可以計算出來的，即：

將L帶入公式2，得到星期w的最終計算公式：

還以2005年5月31日為例，利用公式5計算w的值為：

得到 2005 年 5 月 31 日是星期二，和前面的計算方法得到的結果一致。根據(jù)上述分析，可得寫出使用公式 5 計算星期的算法實現(xiàn)：

146 int TotalWeek(int year, int month, int day)

147 {

148 int d = CalcYearPassedDays(year, month, day);

149 int y = year - 1;

150 int w = y * DAYS_OF_NORMAL_YEAR + y / 4 - y / 100 + y / 400 + d;

151

152 return w % 7;

153 }

公式 5 的問題在于計算量大，不利于口算星期結果。于是人們就在公式 5 的基礎上繼續(xù)推導更簡單的公式。德國數(shù)學家克里斯蒂安·蔡勒（Christian Zeller, 1822- 1899）在 1886 年推導出了著名的蔡勒（Zeller）公式：

對計算出的 w 值除以7，得到的余數(shù)就是星期幾，如果余數(shù)是0，則為星期日。蔡勒公式中各符號的含義如下：

w ：星期；

c ：世紀數(shù) – 1的值，如21世紀，則 = 20；

m ：月數(shù)，的取值是大于等于3，小于等于14。在蔡勒公式中，某年的1月和2月看作上一年的13月和14月，比如2001年2月1日要當成2000年的14月1日計算；

y ：年份，取公元紀念的后兩位，如1998年， = 98，2001年， = 1；

d ：某月內的日數(shù)

為了方便口算，人們通常將公式6中的一項改成。

目前人們普遍認為蔡勒公式是計算某一天是星期幾的最好的公式。但是蔡勒公式有時候可能計算出的結果是負數(shù)，需要對結果+7進行修正。比如2006年7月1日，用蔡勒公式計算出的結果是 -1，實際上這天是星期六。根據(jù)前面分析的結果整理出的蔡勒公式算法實現(xiàn)如下：

155 int ZellerWeek(int year, int month, int day)

156 {

157 int m = month;

158 int d = day;

159

160 if(month <= 2) /*對小于2的月份進行修正*/

161 {

162 year--;

163 m = month + 12;

164 }

165

166 int y = year % 100;

167 int c = year / 100;

168

169 int w = (y + y / 4 + c / 4 - 2 * c + (13 * (m + 1) / 5) + d - 1) % 7;

170 if(w < 0) /*修正計算結果是負數(shù)的情況*/

171 w += 7;

172

173 return w;

174 }

蔡勒公式（公式6）和前面提到的公式 5 都只適用于格里歷法。羅馬教皇在 1582 年修改歷法，將 10 月 5 日指定為 10 月 15 日，從而正式廢止儒略歷法，開始啟用格里歷法。因此，上述求星期幾的公式只適用于 1582 年 10 月 15 日之后的日期，對于 1582 年將 10 月 4 日之前的日期，蔡勒也推導出了適用與儒略歷法的星期計算公式：

公式 7 適用于對 1582 年 10 月 4 日之前的日期計算星期，1582 年 10 月 5 日與 1582 年 10 月 15 日之間的日期是不存在的，因為它們都是同一天。

格里歷歷法簡單，除二月外每月天數(shù)固定，二月則根據(jù)是否是閏年確定是 28 天還是 29 天，每天的星期數(shù)可以通過蔡勒公式（公式 6）計算，有了這些信息，就可以按照一定的排版格式將某一年的日歷打印出來。排版打印的算法非常簡單，就是按照順序打印 12 個月的月歷，因此，打印月歷的函數(shù)就是輸出算法的重點。代碼沒什么特別之處，就是用一些小技巧確定每個月的第一天的開始位置，打印月歷的核心代碼如下：

229 void PrintMonthCalendar(int year, int month)

230 {

231 int days = GetDaysOfMonth(year, month); /*確定這個月的天數(shù)*/

232 if(days <= 0)

233 return;

234

235 PrintMonthBanner(nameOfMonth[month - 1]);

236 PrintWeekBanner();

237 int firstDayWeek = ZellerWeek(year, month, 1);

238 InsertRowSpace(firstDayWeek);

239 int week = firstDayWeek;

240 int i = 1;

241 while(i <= days)

242 {

243 printf("%-10d", i);

244 if(week == 6) /*到一周結束，切換到下一行輸出*/

245 {

246 SetNextRowStart();

247 }

248 i++;

249 week = (week + 1) % 7;

250 }

251 }

GetDaysOfMonth()函數(shù)其實就是從daysOfMonth表中查一下每月的天數(shù)，如果是閏年，則對二月的天數(shù)修正（＋1），daysOfMonth表定義如下：

int daysOfMonth[MONTHES_FOR_YEAR] = { 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};

計算星期不必對每一天都計算一次，只要對每個月的第一天計算一次就可以了，以后的日期可以用 week = (week + 1) % 7 直接推算出星期幾。下面就是我們的算法打印輸出的效果：

********************************************************************************

Calendar of 2012

********************************************************************************

----------January----------

Sunday Monday Tuesday Wednesday Thursday Friday Saturday

1 2 3 4 5 6 7

8 9 10 11 12 13 14

15 16 17 18 19 20 21

22 23 24 25 26 27 28

29 30 31

----------February----------

Sunday Monday Tuesday Wednesday Thursday Friday Saturday

1 2 3 4

5 6 7 8 9 10 11

12 13 14 15 16 17 18

19 20 21 22 23 24 25

26 27 28 29

----------March----------

Sunday Monday Tuesday Wednesday Thursday Friday Saturday

1 2 3

4 5 6 7 8 9 10

11 12 13 14 15 16 17

18 19 20 21 22 23 24

25 26 27 28 29 30 31

……

2、二十四節(jié)氣

二十四節(jié)氣在中國古代歷法中扮演著非常重要的角色，本文將介紹二十四節(jié)氣的基本知識，以及如何使用 VSOP82/87 行星運行理論計算二十四節(jié)氣發(fā)生的準確時間。

中國古代歷法都是以月亮運行規(guī)律為主，嚴格按照朔望月長度定義月，但是由于朔望月長度和地球回歸年長度無法協(xié)調，會導致農(nóng)歷季節(jié)和天氣的實際冷暖無法對應，因此聰明的古人將月亮運行規(guī)律和太陽運行規(guī)律相結合制定了中國農(nóng)歷的歷法規(guī)則。在這種特殊的陰陽結合的歷法規(guī)則中，二十四節(jié)氣就扮演著非常重要的作用，它是聯(lián)系月亮運行規(guī)律和太陽運行規(guī)律的紐帶。正是由于二十四節(jié)氣結合置閏規(guī)則，使得農(nóng)歷的春夏秋冬四季和地球繞太陽運動引起的天氣冷暖變化相一致，成為中國幾千年來生產(chǎn)、生活的依據(jù)。

二十四節(jié)氣起源于中國黃河流域。遠在春秋時代，古人就開始使用仲春、仲夏、仲秋和仲冬四個節(jié)氣指導農(nóng)耕種植。后來經(jīng)過不斷地改進與完善，到秦漢年間，二十四節(jié)氣已經(jīng)基本確立。公元前 104年，漢武帝頒布由鄧平等人制定的《太初歷》，正式把二十四節(jié)氣訂于歷法，明確了二十四節(jié)氣的天文位置。二十四節(jié)氣天文位置的定義，就是從太陽黃經(jīng)零度開始，沿黃經(jīng)每運行15度所經(jīng)歷的時日稱為“一個節(jié)氣”。太陽一個回歸年運行360度，共經(jīng)歷24個節(jié)氣，每個公歷月對應2個節(jié)氣。其中，每月第一個節(jié)氣為“節(jié)令”，即：立春、驚蟄、清明、立夏、芒種、小暑、立秋、白露、寒露、立冬、大雪和小寒等12個節(jié)令；每月的第二個節(jié)氣為“中氣”，即：雨水、春分、谷雨、小滿、夏至、大暑、處暑、秋分、霜降、小雪、冬至和大寒等12個中氣。“節(jié)令”和“中氣”交替出現(xiàn)，各歷時15天，人們習慣上把“節(jié)令”和“中氣”統(tǒng)稱為“節(jié)氣”。

為了更好地理解二十四節(jié)氣的天文位置，首先要解釋幾個天文學概念。“天球”是人們?yōu)榱搜芯刻祗w的位置和運動規(guī)律而引入的一個假象的球體，根據(jù)觀察點（也就是球心）的位置不同，可分為“日心天球”、“地心天球”等等。圖（1）就是天球概念的一個簡單示意圖：

圖（1）天球概念示意圖

天文學中常用的一個坐標體系就是“地心天球”，它與地球同心且有相同的自傳軸，理論上具有無限大的半徑。地球的赤道和南北極點延伸到天球上，對應著天赤道和南北天極點。和地球上用經(jīng)緯度定為位置一樣，天球也劃分了經(jīng)緯度，分別命名為“赤經(jīng)”和“赤緯”，地球上的經(jīng)度用的是度（分秒）為單位，赤經(jīng)以時（分秒）為單位。天空中的所有天體都可以投射到天球上，用赤經(jīng)和赤緯定為天體在天球上的位置。“黃道（Ecliptic）”是地球繞太陽公轉軌道的軌道平面與天球（地心天球）相交的大圓，由于地球公轉受月球和其它行星的攝動，地球的公轉軌道并不是嚴格的平面，因此黃道的嚴格定義是：地月系質心繞太陽公轉的瞬時平均軌道平面與天球相交的大圓。黃道和天赤道所在的兩個平面并不是重疊的，它們之間存在一個23度26分的交角，稱為“黃赤交角”。由于黃赤交角的存在，黃道和天赤道就在天球上有兩個交點，這兩個交點就是春分點和秋分點。在天球上以黃道為基圈可以形成黃道坐標系，在黃道坐標系中，也使用了經(jīng)緯度的概念，分別稱為“黃經(jīng)”和“黃緯”。天體的黃經(jīng)從春分點起沿黃道向東計量，春分點是黃經(jīng)0度，沿黃道一周是360度，使用的單位是度、分和秒。黃緯以黃道測量平面為準，向北記為0度到90度，向南記為0度到-90度。

黃道平面可以近似理解為地球繞太陽公轉的平面，以黃道為基圈的黃道坐標系根據(jù)觀測中心是太陽還是地球還可以區(qū)分為日心坐標系和地心坐標系，對應天體的黃道坐標分別被稱為“日心黃經(jīng)、日心黃緯”和“地心黃經(jīng)、地心黃緯”。日心黃經(jīng)和日心黃緯比較容易理解，因為太陽系的行星都是繞太陽公轉，以太陽為中心將這些行星向天球上投影是最簡單的確定行星位置關系的做法。但是人類自古觀察太陽的周年運動，都是以地球為參照，以太陽的周年視運動位置來計算太陽的運行軌跡，使用的其實都是地心黃經(jīng)和地心黃緯，要了解古代歷法，理解這一點非常重要。圖（2）就解釋了造成這種視覺錯覺的原因：

圖（2）太陽黃道視覺位置原理圖

古人定義二十四節(jié)氣的位置，是太陽沿著黃道運行時的視覺位置，每個節(jié)氣對應的黃道經(jīng)度其實是地心黃經(jīng)。從圖（2）可以看出日心黃經(jīng)和地心黃經(jīng)存在180度的轉換關系，同樣可以理解，日心黃緯和地心黃緯在方向上是反的，因此可以很方便地將兩類坐標相互轉換，轉換公式是：

太陽地心黃經(jīng) = 地球日心黃經(jīng) + 180° （3.1式）

太陽地心黃緯 = -地球日心黃緯（3.2式）

了解了以上的天文學基礎之后，就可以著手對二十四節(jié)氣的發(fā)生時間進行計算。我們常說的節(jié)氣發(fā)生時間，其實就是在太陽沿著黃道做視覺運動過程中，當太陽地心黃經(jīng)等于某個節(jié)氣黃經(jīng)度數(shù)時的那個瞬間的時間。所謂的用天文算法計算二十四節(jié)氣時間，就是根據(jù)牛頓力學原理或開普勒三大行星定律，計算出與歷法密切相關的地球、太陽和月亮三個天體的運行軌道和時間參數(shù)，以此得出當這些天體位于某個位置時的時間。這樣的天文計算需要計算者有扎實的微積分學、幾何學和球面三角學知識，令廣大天文愛好者望而卻步。但是隨著VSOP-82/87行星理論以及ELP-2000/82月球理論的出現(xiàn)，使得天文計算變得簡單易行，本文就是以VSOP-82/87行星理論為計算依據(jù)，計算二十四節(jié)氣的準確時間。

古代天文學家在對包括地球和月亮在內的行星運行軌道精確計算后發(fā)現(xiàn)，天體的運行因為受相近天體的影響，并不嚴格遵循理論方法計算出來的軌道，而是在理論軌道附近波動。這種影響在天文學上被稱為攝動，攝動很難被精確計算，只能根據(jù)經(jīng)驗估算。但是經(jīng)過長期的觀測和計算，天文學家發(fā)現(xiàn)行星軌道因為攝動影響而產(chǎn)生的波動其實也是有規(guī)律的，即在相當長的時間內呈現(xiàn)出周期變化的趨勢。于是天文學家開始研究這種周期變化，希望通過一種類似曲線擬合的方法，對一些周期計算項按照某種計算式迭代求和計算代替積分計算來模擬行星運行軌跡。這種計算式可以描述為：a + bt + ct² + … xcos(p + qt + rt² + …)，其中t是時間參數(shù)，這樣的理論通常被稱為半解析（semi-analytic）理論。其實早在十八世紀，歐洲學者Joseph Louis Lagrange就開始嘗試用這種周期項計算的方法修正行星軌道，但是他采用的周期項計算式是線性方程，精度不高。

1982年，P.Bretagnon公開發(fā)表了VSOP行星理論（這個理論的英文名稱是：Secular Variations of the Planetary Orbits，VSOP的縮寫其實是源于法文名稱：Variations Séculaires des Orbites Planétaires），VSOP理論是一個描述太陽系行星軌道在相當長時間范圍內周期變化的半分析（semi-analytic）理論。VSOP82理論是VSOP理論的第一個版本，提供了對太陽系幾大行星位置計算的周期序列，通過對周期序列進行正弦或余弦項累加求和，就可以得到這個行星在給定時間的軌道參數(shù)。不過VSOP82由于每次都會計算出全部超高精度的軌道參數(shù)，這些軌道參數(shù)對于歷法計算這樣的民用場合很不適用。1987年，Bretagnon 和 Francou 創(chuàng)建了VSOP87行星理論，VSOP87行星理論不僅能計算各種精密的軌道參數(shù)，還可以直接計算出行星的位置，行星位置可以是各種坐標系，包括黃道坐標系。VSOP87行星理論由6張周期項系數(shù)表組成，分別是VSOP87、VSOP87A、VSOP87B、VSOP87C、VSOP87D和VSOP87E，其中VSOP87D表可以直接計算行星日心黃經(jīng)（L）、日心黃緯（B）和到太陽的距離（R），此表計算出的結果適用于節(jié)氣位置判斷。

VSOP87D表包含了三部分數(shù)據(jù)，分別是計算行星日心黃經(jīng)的周期項系數(shù)表（L表）、計算行星日心黃緯的周期項系數(shù)表（B表）和計算行星和太陽距離的周期項系數(shù)表（R表）。VSOP87D表有太陽系8大行星的數(shù)據(jù)，本文的計算只關心與地球相關的數(shù)據(jù)。L表由L0－L5六部分組成，每一部分都包含若干個周期項系數(shù)條目，比如L0表有559個周期項系數(shù)條目，L1表有341個條目等等。L表的每個周期項系數(shù)條目包含若干個參數(shù)，用于計算各種軌道參數(shù)和位置參數(shù)，計算地球的日心黃經(jīng)只需要用到其中三個系數(shù)。計算所有的周期項系數(shù)并不是必須的，有時候減少一些系數(shù)比較小的周期項可以減少計算所花費的時間，當然，這會犧牲一點精度。假設計算地球日心黃經(jīng)的三個系數(shù)是A、B和C，則每個周期項的計算表達式是：

A * cos(B + Cτ) （3.3式）

其中τ是儒略千年數(shù)，τ的計算公式如下：

τ = (JDE - 2451545.0) / 365250 （3.4式）

JDE是計算軌道參數(shù)的時間，單位是儒略日，2451545.0是公元2000年1月1日 12時的儒略日數(shù)，關于儒略日的概念，請參考“日歷生成算法”的第一篇《中國公歷（格里歷）》中的說明以及計算方法。以L0表的第二個周期項為例，這個周期項數(shù)據(jù)中與日心黃經(jīng)計算有關的三個系數(shù)分別是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式計算出最終的地球日心黃經(jīng)，單位是弧度：

L = (L0 + L1 * τ+ L2 * τ² + L3 * τ³ + L4 * τ⁴ +L5 * τ⁵) / 10⁸ （3.5式）

用同樣的方法對地球日心黃緯的周期項系數(shù)表和計算行星和太陽距離的周期項系數(shù)表計算求和，可以得到地球日心黃緯B和日地距離R，B的單位是弧度，R的單位是天文單位（AU）^[1]。由于3.5式的計算方法需要多次計算τ的乘方，浮點數(shù)的乘方計算的速度比較慢，實際計算時，通常對3.5式進行變換，用乘法和加法代替直接的乘方計算，這是一種常用的轉換：

L = (((((L5 * τ + L4) * τ + L3) * τ + L2) * τ + L1) * τ + L0) / 10⁸ （3.6式）

本文就是使用3.6式代替3.5式進行計算。

VSOP82/87行星理論中的周期項系數(shù)對不同的行星具有不同的精度，對地球來說，在1900-2100年之間的200年跨度期間，計算精度是0.005"。前文曾說過，對于不需要這么高精度的計算應用時，可以適當減少一些系數(shù)比較小的周期項，減少計算量，提高計算速度。Jean Meeus在他的《天文算法》一書中就給出了一套精簡后的VSOP87D表的周期項，將計算地球黃經(jīng)的L0表由原來的559項精簡到64項，計算地球黃緯的B0表甚至被精簡到只有5項，從實際效果看，計算精度下降并不多，但是極大地減少了計算量。

使用VSOP87D周期項系數(shù)表計算得到的是J2000.0平黃道和平春分點（mean dynamic ecliptic and equinox）為基準的日心黃經(jīng)（L）和日心黃緯（B），其值與標準FK5系統(tǒng)略有差別，如果對精度要求很高可以采用下面的方法將計算得到的日心黃經(jīng)（L）和日心黃緯（B）轉到FK5系統(tǒng)^[2]：

首先然后 L'，單位是度：

L' = L - 1.397 * T - 0.00031*T² （3.7式）

3.7式中的T是儒略世紀數(shù)，它與儒略千年數(shù)τ的關系是：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()函數(shù)就是使用VSOP87行星理論計算行星日心黃經(jīng)的代碼實現(xiàn)，整個計算過程和前文描述一樣，首先根據(jù)VSOP87D表的數(shù)據(jù)計算出L0-L5，然后用3.6式計算出地球的日心黃經(jīng)，3.6式計算出來的單位是弧度，因此轉換成度分秒單位，最后使用3.1式將結果轉換成太陽的地心黃經(jīng)：

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;

/*地心黃經(jīng) = 日心黃經(jīng) + 180度*/

return (Mod360Degree(Mod360Degree(L / RADIAN_PER_ANGLE) + 180.0));

}

Mod360Degree()函數(shù)將大于360°或小于0°的值調整到0－360°之間，便于轉換顯示。CalcPeriodicTerm()函數(shù)使用3.3式對一個周期項系數(shù)表進行求和計算，可以指定需要計算的周期項數(shù)：

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()函數(shù)根據(jù)3.8式計算黃經(jīng)的修正量，longitude和latitude參數(shù)是由VSOP87理論計算出的太陽地心黃經(jīng)和地心黃緯，單位是度，dt是儒略千年數(shù)，返回值單位是度：

double AdjustSunEclipticLongitudeEC(double dt, double longitude, double latitude)

{

double T = dt * 10; //T是儒略世紀數(shù)

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;

}

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

【下篇將介紹修正理論和修正算法，請繼續(xù)關注】

小知識1：天文單位

天文單位（英文：Astronomical Unit，簡寫AU）是一個長度的單位，約等于地球跟太陽的平均距離。天文單位是天文常數(shù)之一，是天文學中測量距離，特別是測量太陽系內天體之間的距離的基本單位。地球到太陽的平均距離大約為一個天文單位，約等于1.496億千米。 1976年，國際天文學聯(lián)會把一天文單位定義為一顆質量可忽略、公轉軌道不受干擾而且公轉周期為365.2568983日（即一高斯年）的粒子與一個質量相等約一個太陽的物體的距離。當前普遍被接受并使用的天文單位的值是149,597,870,691±30米（約一億五千萬公里）。

小知識2：FK5系統(tǒng)

FK5常用的目視星表系統(tǒng)，又稱第五基本星表，是在FK4（第四基本星表）的基礎上發(fā)展出來的，對FK4星表進行了修正，于1984年正式啟用。它定義了一個以太陽質心為中心，J2000.0平赤道和春分點為基準的天球平赤道坐標系。近年來國際上又編制了FK6星表（第六基本星表），但是還沒有被正式啟用。

首先需要進行章動修正。章動是指地球沿自轉軸的指向繞黃道極緩慢旋轉過程中，由于地球上物質分布不均勻性和月球及其它行星的攝動力造成的輕微抖動。英國天文學家詹姆斯·布拉德利(1693—1762)最早發(fā)現(xiàn)了章動，章動可以沿著黃道分解為水平分量和垂直分量，黃道上的水平分量記為Δψ，稱為黃經(jīng)章動，它影響了天球上所有天體的經(jīng)度。黃道上的垂直分量記為Δε，稱為交角章動，它影響了黃赤交角。目前編制天文年歷所依據(jù)的章動理論是伍拉德在1953年建立的，它是以剛體地球模型為基礎的。1977年，國際天文聯(lián)合會的一個專家小組建議采用非剛體地球模型――莫洛堅斯基II模型代替剛體地球模型計算章動，1979年的國際天文學聯(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ù)（相位內項系數(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是儒略世紀數(shù)，計算出來的5個基本角距的單位都是度，在計算正弦或余弦時要轉換為弧度單位。計算每一個周期項的黃經(jīng)章動過程是這樣的，首先將3.10－3.14式計算出來的值與對應的5個基本角距系數(shù)組合，計算出輻角。以本文使用的章動周期項系數(shù)表中的第七項為例，5個基本角距對應的系數(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)代碼：

double CalcEarthLongitudeNutation(double dt)

{

double T = dt * 10;

double D,M,Mp,F,Omega;

GetEarthNutationParameter(dt, &D, &M, &Mp, &F, &Omega);

double resulte = 0.0 ;

for(int i = 0; i < sizeof(nutation) / sizeof(nutation[0]); i++)

{

double sita = nutation[i].D * D + nutation[i].M * M + nutation[i].Mp * Mp +nutation[i].F * F + nutation[i].omega * Omega;

resulte += (nutation[i].sine1 + nutation[i].sine2 * T ) * sin(sita);

}

/*先乘以章動表的系數(shù) 0.0001，然后換算成度的單位*/

return resulte * 0.0001 / 3600.0;

}

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

void GetEarthNutationParameter(double dt, double *D, double *M, double *Mp, double*F, double *Omega)

{

double T = dt * 10; /*T是從J2000起算的儒略世紀數(shù)*/

double T2 = T * T;

double T3 = T2 * T;

/*平距角（如月對地心的角距離）*/

*D = 297.85036 + 445267.111480 * T - 0.0019142 * T2 + T3 / 189474.0;

/*太陽（地球）平近點角*/

*M = 357.52772 + 35999.050340 * T - 0.0001603 * T2 - T3 / 300000.0;

/*月亮平近點角*/

*Mp = 134.96298 + 477198.867398 * T + 0.0086972 * T2 + T3 / 56250.0;

/*月亮緯度參數(shù)*/

*F = 93.27191 + 483202.017538 * T - 0.0036825 * T2 + T3 / 327270.0;

/*黃道與月亮平軌道升交點黃經(jīng)*/

*Omega = 125.04452 - 1934.136261 * T + 0.0020708 * T2 + T3 / 450000.0;

}

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

double CalcEarthObliquityNutation(double dt)

{

double T = dt * 10; /*T是從J2000起算的儒略世紀數(shù)*/

double D,M,Mp,F,Omega;

GetEarthNutationParameter(dt, &D, &M, &Mp, &F, &Omega);

double resulte = 0.0 ;

for(int i = 0; i < sizeof(nutation) / sizeof(nutation[0]); i++)

{

double sita = nutation[i].D * D + nutation[i].M * M + nutation[i].Mp * Mp +nutation[i].F * F + nutation[i].omega * Omega;

resulte += (nutation[i].cosine1 + nutation[i].cosine2 * T ) * cos(sita);

}

/*先乘以章動表的系數(shù) 0.001，然后換算成度的單位*/

return resulte * 0.0001 / 3600.0;

}

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

double CalcSunEarthRadius(double dt)

{

double R0 = CalcPeriodicTerm(Earth_R0, sizeof(Earth_R0) /sizeof(VSOP87_COEFFICIENT), dt);

double R1 = CalcPeriodicTerm(Earth_R1, sizeof(Earth_R1) /sizeof(VSOP87_COEFFICIENT), dt);

double R2 = CalcPeriodicTerm(Earth_R2, sizeof(Earth_R2) /sizeof(VSOP87_COEFFICIENT), dt);

double R3 = CalcPeriodicTerm(Earth_R3, sizeof(Earth_R3) /sizeof(VSOP87_COEFFICIENT), dt);

double R4 = CalcPeriodicTerm(Earth_R4, sizeof(Earth_R4) /sizeof(VSOP87_COEFFICIENT), dt);

double R = (((((R4 * dt) + R3) * dt + R2) * dt + R1) * dt + R0) / 100000000.0;

return R;

}

也可以不使用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都是儒略世紀數(shù)，M和C的單位都是度，帶入3.18式計算時需要轉換成弧度單位，計算出R以后，就可以這樣計算光行差修正量：

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

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

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

double GetSunEclipticLongitudeEC(double jde)

{

double dt = (jde - JD2000) / 365250.0; /*儒略千年數(shù)*/

// 計算太陽的地心黃經(jīng)

double longitude = CalcSunEclipticLongitudeEC(dt);

// 計算太陽的地心黃緯

double latitude = CalcSunEclipticLatitudeEC(dt) * 3600.0;

// 修正精度

longitude += AdjustSunEclipticLongitudeEC(dt, longitude, latitude);

// 修正天體章動

longitude += CalcEarthLongitudeNutation(dt);

// 修正光行差

/*太陽地心黃經(jīng)光行差修正項是: -20".4898/R*/

longitude -= (20.4898 / CalcSunEarthRadius(dt)) / (20 * PI);

return longitude;

}

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

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

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

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

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

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

34 double CalculateSolarTerms(int year, int angle)

35 {

36 double lJD, rJD;

37 EstimateSTtimeScope(year, angle, lJD, rJD); /*估算迭代起始時間區(qū)間*/

39 double solarTermsJD = 0.0;

40 double longitude = 0.0;

42 do

43 {

44 solarTermsJD = ((rJD - lJD) * 0.618) + lJD;

45 longitude = GetSunEclipticLongitudeECDegree(solarTermsJD);

46 /*

47 對黃經(jīng)0度迭代逼近時，由于角度360度圓周性，估算黃經(jīng)值可能在(345,360]和[0,15)兩個區(qū)間，

48 如果值落入前一個區(qū)間，需要進行修正

49 */

50 longitude = ((angle == 0) && (longitude > 345.0)) ? longitude - 360.0 :longitude;

52 (longitude > double(angle)) ? rJD = solarTermsJD : lJD = solarTermsJD;

53 }while((rJD - lJD) > 0.0000001);

55 return solarTermsJD;

56 }

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

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

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

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

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

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

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

74 double CalculateSolarTermsNewton(int year, int angle)

75 {

76 double JD0, JD1,stDegree,stDegreep;

78 JD1 = GetInitialEstimateSolarTerms(year, angle);

79 do

80 {

81 JD0 = JD1;

82 stDegree = GetSunEclipticLongitudeECDegree(JD0) - angle;

83 stDegreep = (GetSunEclipticLongitudeECDegree(JD0 + 0.000005)

84 - GetSunEclipticLongitudeECDegree(JD0 - 0.000005)) /0.00001;

85 JD1 = JD0 - stDegree / stDegreep;

86 }while((fabs(JD1 - JD0) > 0.0000001));

88 return JD1;

89 }

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

至此，我們就有了完整的計算節(jié)氣發(fā)生時間的方法，輸入年份和節(jié)氣對應的太陽黃經(jīng)度數(shù)，即可求的節(jié)氣發(fā)生的精確時間。最后說明一下，以上算法中討論的時間都是力學時時間（TD），與國際協(xié)調時（UTC）以及各個時區(qū)的本地時間都有不同，以上計算出來的時間都需要調整成本地時間，比如中國的中原地區(qū)就是東八區(qū)標準時（UTC + 8）。關于力學時、國際協(xié)調時（世界時）的定義，請參考文末的小知識3：力學時、原子時和國際協(xié)調時。應用本文的算法計算出2012年各個節(jié)氣的時間如下（已經(jīng)轉換為東八區(qū)標準時），與紫金山天文臺發(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：力學時、原子時和國際協(xié)調時

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

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

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

參考文章：

[1] 《Secular variations of the planetary orbits》http://www.worldlingo.com/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

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

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

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

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

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

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

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

4.1式中的四個輻角系數(shù)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是儒略世紀數(shù)，T計算由4.6式計算：

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

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

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

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

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

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

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

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

由于金星和木星對月球的攝動影響，需要對計算出的月球黃經(jīng)周期項和Σ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’是月球平黃經(jīng)，計算方法是：

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

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

β = Σb / 1000000.0 （4.17式）

ΣI和Σb最后要除以1000000.0是因為周期項系數(shù)中振幅A的單位是0.000001度，因此λ和β的單位是度。下面給出計算月球地心視黃經(jīng)的代碼：

123 double GetMoonEclipticLongitudeEC(double dbJD)

124 {

125 double Lp,D,M,Mp,F,E;

126 double dt = (dbJD - JD2000) / 36525.0; /*儒略世紀數(shù)*/

127

128 GetMoonEclipticParameter(dt, &Lp, &D, &M, &Mp, &F, &E);

129

130 /*計算月球地心黃經(jīng)周期項*/

131 double EI = CalcMoonECLongitudePeriodicTbl(D, M, Mp, F, E);

132

133 /*修正金星,木星以及地球扁率攝動*/

134 EI += CalcMoonLongitudePerturbation(dt, Lp, F);

135

136 /*計算月球地心黃經(jīng)*/

137 double longitude = Lp + EI / 1000000.0;

138

139 /*計算天體章動干擾*/

140 longitude += CalcEarthLongitudeNutation(dt / 10.0);

141

142 longitude = Mod360Degree(longitude);

143 return longitude;

144 }

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

42 double CalcMoonECLongitudePeriodicTbl(double D, double M, double Mp, double F,double E)

43 {

44 double EI = 0.0 ;

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 }

53 return EI;

54 }

CalcMoonLongitudePerturbation()函數(shù)計算月球黃經(jīng)攝動修正量，使用了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?數(shù)'ca?*/

90 double A1 = 119.75 + 131.849 * T;

91 double A2 = 53.09 + 479264.290 * T;

93 A1 = Mod360Degree(A1);

94 A2 = Mod360Degree(A2);

96 double result = 3958.0 * sin(DegreeToRadian(A1));

97 result += (1962.0 * sin(DegreeToRadian(Lp - F)));

98 result += (318.0 * sin(DegreeToRadian(A2)));

100 return result;

101 }

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

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

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

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

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

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

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 }

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

再次說明一下，以上算法中討論的時間都是力學時時間（TD），與國際協(xié)調時（UTC）以及各個時區(qū)的本地時間都有不同，以上計算出來的時間都需要調整成本地時間，比如中國的中原地區(qū)就是東八區(qū)標準時（UTC + 8）。應用本文的算法計算出2012年前后15個日月合朔時間如下（已經(jīng)轉換為東八區(qū)標準時）：

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

世界各國的日歷都是以天為最小單位，但是關于年和月的算法卻各不相同，大致可以分為三類：

陽歷--以天文年作為日歷的主要周期，例如：中國公歷（格里歷）

陰歷--以天文月作為日歷的主要周期，例如：伊斯蘭歷

陰陽歷--以天文年和天文月作為日歷的主要周期，例如：中國農(nóng)歷

我國古人很早就開始關注天象，定晝夜交替為“日”，月輪盈虧為“月”，寒暑交替為“年”，在總結日月變化規(guī)律的基礎上制定了兼有陰歷月和陽歷年性質的歷法，稱為中國農(nóng)歷。本文將介紹中國農(nóng)歷的歷法規(guī)則、天干地支（Heavenly Stems，Earthly Branches）的計算方法以、二十四節(jié)氣與中國農(nóng)歷的關系以及知道節(jié)氣和日月合朔的精確時間的情況下推算中國農(nóng)歷年歷的方法。

在介紹中國農(nóng)歷的歷法之前，必須要先介紹一下中國古代的紀年方法。中國古代用天干地支紀年，嚴格來講，天干地支紀年以及十二屬相并不是中國農(nóng)歷歷法的一部分，但是在中國歷史上直到今天，天干地支以及十二屬相一直都是做為中國農(nóng)歷紀年關系密切的一部分而存在，因此這里先介紹一下天干地支紀年法以及十二屬相。

中國古代紀年不用數(shù)字，而是采用天干地支組合。天干有十個，分別是：甲、乙、丙、丁、戊、己、庚、辛、壬、癸；地支有十二個，分別是：子、丑、寅、卯、辰、巳、午、未、申、酉、戌、亥。使用時天干地支各取一字，天干在前，地支在后，組合成干支，例如甲子、乙丑、丙寅等等，依次輪回可形成六十種組合，以這些天干地支組合紀年，每六十年一個輪回，稱為一個甲子。實際上中國古代紀月、紀日以及紀時辰都采用干支方法，這些干支組合起來就是我們熟悉的生辰八字。十二屬相又稱“十二生肖”，由十一種源自自然界的動物：鼠、牛、虎、兔、蛇、馬、羊、猴、雞、狗、豬以及傳說中的龍組成，用于紀年時，按順序和十二地支組合成子鼠、丑牛、寅虎、卯兔、辰龍、巳蛇、午馬、未羊、申猴、酉雞、戌狗和亥豬。天干地支以及十二生肖常組合起來描述農(nóng)歷年，比如公歷2011年就是農(nóng)歷辛卯兔年、2012年是壬辰龍年等等。

計算某一年的天干地支，有很多經(jīng)驗公式，如果知道某一年的天干地支，也可以直接推算其它年份的天干地支。舉個例子，如果知道2000年是庚辰龍年，則2012年的干支可以這樣推算：（2012-2000）% 10=2，2012年的天干就是從庚開始向后推2個天干，即壬；2012年的地支可以這樣推算：（2012 - 2000）% 12 = 0，2012年的地支仍然是辰，因此2012年的天干地支就是壬辰，十二生肖龍年。對于2000年以前的年份，計算出年份差后只要將天干和地支向前推算即可。例如1995年的干支可以這樣計算：（2000 – 1995）%10 = 5，（2000 – 1995）%12 = 5，庚向前推算5即是乙，辰向前推算5即是亥，因此1995年的干支就是乙亥，十二生肖豬年。這個干支推算算法的實現(xiàn)如下：

202 void CalculateYearGanZhi(int year, int *gan, int *zhi)

203 {

204 int sc = year - 2000;

205 *gan = (7 + sc) % 10;

206 *zhi = (5 + sc) % 12;

207

208 if(*gan < 0)

209 *gan += 10;

210 if(*zhi < 0)

211 *zhi += 12;

212 }

獲得2008年的干支紀年：

9 TCHAR *nameOfTianGan[COUNTS_FOR_TIANGAN] = {_T("甲"),_T("乙"),_T("丙"),_T("丁"),_T("戊"),_T("己"),_T("庚"),_T("辛"),_T("壬"),_T("癸") };

10 TCHAR *nameOfDiZhi[COUNTS_FOR_DIZHI] = {_T("子"),_T("丑"),_T("寅"),_T("卯"),_T("辰"),_T("巳"),_T("午"),_T("未"),_T("申"),_T("酉"),_T("戌"),_T("亥")};

146 int gan,zhi;

147

148 CalculateYearGanZhi(2008, &gan, &zhi);

149

150 text.Format(_T("農(nóng)歷【%s%s】%s年"),

151 year, m_curMonth, nameOfTianGan[gan - 1], nameOfDiZhi[zhi - 1], nameOfShuXiang[zhi- 1]);

結果是：農(nóng)歷戊子鼠年。

中國農(nóng)歷是以月亮運行周期為基礎，結合太陽運行規(guī)律（二十四節(jié)氣）制定的歷法，農(nóng)歷月的定義規(guī)則就是中國農(nóng)歷歷法的關鍵，因此要了解中國農(nóng)歷的歷法規(guī)則，就必須知道如何定義月，如何設置閏月？中國農(nóng)歷的一年有十二個月或十三個月，但是正統(tǒng)的叫法只有十二個月，分別是正月、二月、三月、四月、五月、六月、七月、八月、九月、十月、冬月和臘月（注意，正統(tǒng)的中國農(nóng)歷是沒有十一月和十二月的，如果你用的歷法軟件有顯示農(nóng)歷十一月和農(nóng)歷十二月，就說明非常不專業(yè)）。中國民間常用“十冬臘月天”來形容寒冷的天氣，其實指的就是十月，十一月和十二月這三個最冷的月份。一年有十三個月的情況是因為有閏月，多出來的這個閏月沒有月名，只是跟在某個月后面，稱為閏某月。比如公歷2009年對應的農(nóng)歷乙丑年，就是閏五月，于是這一年可以過兩個端午節(jié)。

中國農(nóng)歷為什么會有閏月？其實中國農(nóng)歷置閏月是為了協(xié)調回歸年和農(nóng)歷年的矛盾。前面提到過，中國農(nóng)歷是一種陰陽歷，農(nóng)歷的月分大月和小月，大月一個月是30天，小月一個月是29天。中國農(nóng)歷把日月合朔（太陽和月亮的黃經(jīng)相同，但是月亮不可見）的日期定位月首，也就是“初一”，把月圓的時候定為望日，也就是“十五”，月亮繞地球公轉一周稱為一個朔望月。天文學的朔望月長度是29.5306日，中國農(nóng)歷以朔望月為基礎，嚴格保證每個月的頭一天是朔日，這就使得每個月是大月還是小月的安排不能固定，通常需要通過天文學觀測和計算來確定。一個農(nóng)歷年由12個朔望月組成，這樣一個農(nóng)歷年的長度就是29.5306 12 = 354.3672日，而陽歷的一個天文學回歸年是365.2422日，這樣一個農(nóng)歷年就比一個回歸年少10.88天，這個誤差如果累計起來過16年就會出現(xiàn)“六月飛雪”的奇觀了。為了協(xié)調農(nóng)歷年和回歸年之間的矛盾，聰明的先人在天文觀測的基礎上，找到了“閏月”的方法，通過在適當?shù)脑路莶迦腴c月來保證每個農(nóng)歷年的正月到三月是春季，四月到六月是夏季，七月到九月是秋季，十月到十二月是冬季，也就是說，讓歷法和天文氣象能夠基本對上，不至于出現(xiàn)“六月飛雪”。

那么多長時間增加一個閏月比較合適呢？最早人們推算是“三年一閏”，后來是“五年兩潤”，隨著歷法計算的精確，最終定型為“十九年七閏”。這個“十九年七閏”又是怎么算出來的呢？其實就是求出回歸年日數(shù)和朔望月日數(shù)的最小公倍數(shù)，也就是m個回歸年的天數(shù)和n個朔望月的天數(shù)相等，即：

m 365.2422 = n 29.5306

這樣m和n的比例就是29.5306 : 365.2422 19 : 235，按照這個最接近的整數(shù)倍數(shù)關系，每19個回歸年需要添加的閏月就是：

235 – 12 19 = 7

也就是“十九年七閏”的由來。但是需要注意的是，“十九年七閏”也并不是精確的結果，每19年就會有0.0892天的誤差：

19 365.2422 - 235 29.5306 0.0892

這樣每213年就會積累約1天的誤差，因此，即使按照“十九年七閏”計算，中國農(nóng)歷每一兩百年就需要修正一次。正因為這樣，現(xiàn)行農(nóng)歷從唐代以后就已經(jīng)不再遵守“十九年七閏”法，而是采用更準確的“中氣置閏”法。“中氣置閏”法更準確的名稱應該是“定冬至”法，就是定兩個冬至節(jié)氣之間的時間為一個農(nóng)歷年，這樣農(nóng)歷年的長度就和太陽回歸年長度對應，不會產(chǎn)生誤差。

現(xiàn)在，我們知道農(nóng)歷通過置閏月的方式協(xié)調農(nóng)歷年和回歸年長度不相等的問題，也知道了置閏的方法是“中氣置閏”法，那么到底什么是“中氣”，又是如何定中氣置閏月呢？要回答這個問題，就需要介紹另一個天文現(xiàn)象――節(jié)氣。二十四節(jié)氣起源于黃河流域，遠在春秋時代，就定出仲春、仲夏、仲秋和仲冬等四個節(jié)氣。以后不斷地改進與完善，到秦漢年間，二十四節(jié)氣已完全確立，漢武帝太初元年（公元前104年）制定的《太初歷》，則第一次從歷法上明確了二十四節(jié)氣的天文位置。

地球沿著一個近似橢圓軌道繞太陽公轉，這個公轉軌道所在的平面就是“黃道面”，黃道面向外延伸與天球的交線就是“黃道”。古人由于觀測條件限制，只能根據(jù)視覺感覺認為是太陽沿著黃道繞地球運轉，因此設定太陽從黃經(jīng)（黃道經(jīng)度）零度起（以春分點為起點自西向東度量），將太陽沿黃經(jīng)每運行15度所經(jīng)歷的時日稱為“一個節(jié)氣”。太陽每年運行360度，共經(jīng)歷二十四個節(jié)氣，春季的節(jié)氣有立春（315度）、雨水（330度）、驚蟄（345度）、春分（0度、360度）、清明（15度）和谷雨（30度），夏季的節(jié)氣有立夏（45度）、小滿（60度）、芒種（75度）、夏至（90度）、小暑（105度）和大暑（120度），秋季的節(jié)氣有立秋（135度）、處暑（150度）、白露（165度）、秋分（180度）、寒露（195度）和霜降（210度）。冬季的節(jié)氣有立冬（225度）、小雪（240度）、大雪（255度）、冬至（270度）、小寒（285度）和大寒（300度）。二十四節(jié)氣又細分為十二節(jié)氣和十二中氣，二十四節(jié)氣按照順序排在奇數(shù)位置上的就是節(jié)氣，排在偶數(shù)位置上的就是中氣。也就是說，立春、驚蟄、清明、立夏、芒種、小暑、立秋、白露、寒露、立冬、大雪和小寒就是十二個節(jié)氣，而雨水、春分、谷雨、小滿、夏至、大暑、處暑、秋分、霜降、小雪、冬至和大寒就是十二個中氣。二十四個節(jié)氣平分在公歷的12個月中，每月一節(jié)氣一中氣。二十四節(jié)氣反映了太陽的周年運動（以地球為參照物的視運動），所以節(jié)氣在現(xiàn)行的公歷中日期基本固定，上半年在6日、21日，下半年在8日、23日，前后不差 1～2天。中國民間流傳的《二十四節(jié)氣歌》就是為了方便記憶這些節(jié)氣：

春雨驚春清谷天，

夏滿芒夏暑相連，

秋處露秋寒霜降，

冬雪雪冬小大寒，

每月兩節(jié)不變更，

最多相差一兩天。

傳統(tǒng)上一個農(nóng)歷年起于冬至，終于冬至，因此要確定在哪一年置閏，主要看那一年兩個冬至之間有幾個朔望月，如果是12個朔望月，則不置閏，如果是十三個朔望月，則置閏月，至于閏幾月，則要看節(jié)氣而定。對于有13個朔望月的農(nóng)歷年，置閏月的規(guī)則就是從農(nóng)歷二月開始到十月，第一個沒有中氣的月就是閏月，這個沒有中氣的朔望月跟在哪個月后面就是閏幾月。為什么會有沒有中氣的朔望月呢？黃道上兩個中氣之間相隔30度，一個回歸年的長度是365.2422日，則兩個中氣之間的平均間隔是365.2422 12 = 30.4368日，但是因為地球軌道是橢圓軌道，因此相鄰的兩個中氣的時間間隔是不均勻的，比如在遠地點附近的中氣間隔就會長一點，最長可能是31.45天。而農(nóng)歷的朔望月平均長度是29.5306日，這樣就會出現(xiàn)某個朔望月剛好落在兩個中氣之間的情況，比如，某個月的上一個月月末是一個中氣，但是下一個中氣落在這個月的下一個月的頭幾天里，這樣這個月就沒有中氣了。舉個例子，2001年農(nóng)歷辛已年的四月二十九（公歷5月21日）是小滿，農(nóng)歷四月之后的這個朔望月從公歷5月23日持續(xù)到公歷6月20日，而小滿后的下一個中氣夏至是在公歷的6月21日，也就是農(nóng)歷四月的下下個月的初一，這樣農(nóng)歷四月后的這個月就沒有中氣，跟在四月之后，就稱為閏四月。

由于節(jié)氣在回歸年中是均勻分布的，因此公歷中的節(jié)氣日期基本上是固定的，比如立春是在公歷的2月3－5日，不會超出這個日期范圍，這也就是《二十四節(jié)氣歌》所說的：每月兩節(jié)不變更，最多相差一兩天。但是在中國農(nóng)歷中哪個中氣屬于哪個月是有規(guī)定的，雨水是正月的中氣，春分是二月的中氣，谷雨是三月的中氣，小滿是四月的中氣，夏至是五月的中氣，大暑是六月的中氣，處暑是七月的中氣，秋分是八月的中氣，霜降是九月的中氣，小月是十月的中氣，冬至是十一月的中氣，大寒是十二月的中氣。

在了解了農(nóng)歷與節(jié)氣的關系以及農(nóng)歷如何置閏月的方法之后，還需要解決一個問題才能著手農(nóng)歷年歷的推算，那就是如何確定農(nóng)歷年的開始，或者說哪個月的初一是農(nóng)歷新年的開始？要回答這個問題，就需要了解中國農(nóng)歷特有的“月建”問題。

中國農(nóng)歷是陰陽合歷，需要同時考慮太陽和月亮的位置。所以在確定歲首（元旦）時，需要先確定它在某個季節(jié)，然后再選定與這個季節(jié)相近的朔望月作為歲首。由于一歲（一個回歸年）和12個陰歷月并不相等，相差約10.88天，因此每隔三年需要設置一個閏月調整季節(jié)。中國上古的天文學家想出了一個簡便的方法判斷月序與季節(jié)的關系，這就是以傍晚時北斗七星的斗柄的指向確定月序，稱為“十二月建”。從北方起向東轉，將地面劃分為十二個方位，傍晚時北斗所指的方位，就是該月的月建，其子月為冬至所在之月，對應十一月，丑月是冬至所在之月的次月，對應十二月，寅月在丑月之后，對應正月。中國在歷史上的不同時期，多次修改過歲首（元旦）的起始月份，上古時代就有“三正”之說，所謂“三正”，就是“夏正建寅、殷正建丑、周正建子”，意思是夏歷以寅月（正月）為歲首，殷歷以丑月（十二月）為歲首，周歷以子月（十一月）為歲首。從秦代到西漢前期又采用秦歷，秦歷建亥，也就是以亥月作為歲首之月，漢武帝太初元年（公元104年）改用太初歷，重新適用建寅的夏歷，以寅月（正月）為歲首。在這之后的兩千多年時間里，除王莽和魏明帝一度改用建丑的殷歷，唐武后和肅宗時改用建子的周歷外，各個朝代均使用建寅的夏歷直到清朝末年。辛亥革命勝利以后，南京國民政府將公歷1月1日改為元旦，但是人們仍習慣稱農(nóng)歷的正月初一為元旦。新中國成立初期召開的第一屆政治協(xié)商會議，正式將公歷的1月1日確定為元旦，將農(nóng)歷的正月初一定為“春節(jié)”，也就是說，農(nóng)歷的歲首仍然采用夏歷從寅月（正月）開始。

了解了“月建”問題，就解決了農(nóng)歷朔望月與公歷月的對應關系，那就是冬至節(jié)氣所在的朔望月就是農(nóng)歷的子月，對于目前適用的夏歷建寅的月建體系，就意味著冬至節(jié)氣所在的朔望月是農(nóng)歷的十一月，只要找到這個朔望月的起始日（日月合朔發(fā)生的時刻所在的那一日），就找到了公歷的日期月農(nóng)歷日期的對應關系。下面總結一下中國農(nóng)歷歷法的基本法則：

1、嚴格以日月合朔發(fā)生時刻為月首，這一天定為初一，通過計算兩次日月合朔的時間間隔確定每月是29天還是30天；

2、月以中氣得名，冬至節(jié)氣總是出現(xiàn)在農(nóng)歷十一月，包含雨水中氣的月為正月（即寅月），月無中氣者為閏月，與前一個月同名；

3、從某一年的冬至后第一天開始，到下一個冬至這段時間內，如果有十三個朔望月出現(xiàn)，則此期間要增加一個閏月，從二月到十月，第一個沒有中氣的月就是閏月，如果在此期間有超過兩個朔望月沒有中氣，則只有第一個沒有中氣的朔望月是閏月；

4、農(nóng)歷年以正月初一為歲首（關于農(nóng)歷歲首的說法，請參考文末附加的《小知識5：正月初一和立春節(jié)氣》），以臘月（十二月）廿九或三十為除夕；

5、如果節(jié)氣和日月合朔在同一天，則該節(jié)氣是這個新朔望月的節(jié)氣。（民間歷法）

規(guī)則5對節(jié)氣和朔日在同一天的處理，采用了民間歷法的處理原則，關于民間歷法和歷理歷法的區(qū)別，請參考文末附加的《小知識1：民間歷法和歷理歷法》。

了解了農(nóng)歷歷法的基本法則后，就可以根據(jù)歷法進行農(nóng)歷年歷的推算。農(nóng)歷年歷的推算是一件很復雜的事情，需要知道每年二十四個節(jié)氣和本年內每次日月合朔的精確時間，這些時間的獲取比較困難?，F(xiàn)在有很多可以顯示農(nóng)歷的日歷軟件，其實并不計算這些時間，而是事先從權威機構（如紫金山天文臺）獲取這些經(jīng)過推算的時間，然后用各種方法將這些信息存儲在設計好的數(shù)據(jù)結構中。當計算農(nóng)歷時采用查表的方法獲取每年的二十四節(jié)氣日期、大小月情況以及閏月情況，這樣的軟件受數(shù)據(jù)量的限制，往往只能顯示近一兩百年的年歷。

還有一種確定節(jié)氣時間和朔日時間的方法，就是在已知某個節(jié)氣或朔日的精確時間后，通過某些規(guī)律先前或向后推算其它節(jié)氣或朔日的時間。有一些經(jīng)驗公式可以用來計算節(jié)氣發(fā)生的日期，比如“通式壽星公式”，可以計算出某一年的某個節(jié)氣時間，但是只能精確到日。關于“通式壽星公式”的詳細內容，請參考文末附加的《小知識2：通式壽星公式》。至于精確的節(jié)氣或朔日時間，也只能從權威機構獲取。以節(jié)氣的時間推算為例，二十四個節(jié)氣就是黃道上的24各點，由于地球運動受其它天體的影響，導致這些節(jié)氣在每年的時間是不固定的，但是這些節(jié)氣之間的間隔時間基本上可以看作是固定的，下表就是二十四節(jié)氣的時間間隔表：

節(jié)氣名	與上一節(jié)氣之間的時間差	與小寒節(jié)氣的累積時間差
小寒	1271448.00	0.00
大寒	1272494.40	1272494.40
立春	1275526.20	2548020.60
雨水	1282123.20	3830143.80
驚蟄	1290082.80	5120226.60
春分	1300639.20	6420865.80
清明	1311153.00	7732018.80
谷雨	1323253.80	9055272.60
立夏	1333685.40	10388958.00
小滿	1344107.40	11733065.40
芒種	1351227.00	13084292.40
夏至	1357299.60	14441592.00
小暑	1358968.80	15800560.80
大暑	1358786.40	17159347.20
立秋	1354419.00	18513766.20
處暑	1348236.00	19862002.20
白露	1339003.20	21201005.40
秋分	1328654.40	22529659.80
寒露	1317185.40	23846845.20
霜降	1305760.80	25152606.00
立冬	1295081.40	26447687.40
小雪	1285764.00	27733451.40
大雪	1278469.80	29011921.20
冬至	1273556.40	30285477.60

表（1）二十四節(jié)氣時間間隔表（單位：秒鐘）

已知1900年小寒時刻為1月6日2:05:00，以這個節(jié)氣時刻為基準，推算其它年份節(jié)氣的算法實現(xiàn)如下：

8 static double s_stAccInfo[] =

9 {

10 0.00, 1272494.40, 2548020.60, 3830143.80, 5120226.60, 6420865.80,

11 7732018.80, 9055272.60, 10388958.00, 11733065.40, 13084292.40,14441592.00,

12 15800560.80, 17159347.20, 18513766.20, 19862002.20, 21201005.40,22529659.80,

13 23846845.20, 25152606.00, 26447687.40, 27733451.40, 29011921.20,30285477.60

14 };

16 //已知1900年小寒時刻為1月6日02:05:00

17 const double base1900_SlightColdJD = 2415025.5868055555;

19 double CalculateSolarTermsByExp(int year, int st)

20 {

21 if((st < 0) || (st > 24))

22 return 0.0;

24 double stJd = 365.24219878 * (year - 1900) + s_stAccInfo[st] / 86400.0;

26 return base1900_SlightColdJD + stJd;

28 }

base1900_SlightColdJD是北京時間1900年1月6日凌晨2:05:00的儒略日數(shù)，CalculateSolarTermsByExp()函數(shù)返回指定年份的節(jié)氣的儒略日數(shù)。已知某個朔日的精確時間推算其它朔日時間的方法也類似，以朔望月的長度為單位向前或向后累加即可。

這種推算的方法是建立在地球回歸年的長度是固定365.2422天、節(jié)氣的間隔是絕對固定的、朔望月長度是平均的29.5305天等假設之上的，由于天體運動的互相影響，這種假設不是絕對成立的，因此這種推算方法的誤差很大。以CalculateSolarTermsByExp()函數(shù)為例，計算1900年前后30年內的節(jié)氣時間的誤差還可以控制在30分鐘以內，但是到2000年的時候誤差已經(jīng)超過130分鐘了。人們還總結出了計算節(jié)氣和朔日時間的兩個經(jīng)驗公式，本文末尾附加的《小知識3：計算節(jié)氣和朔日的經(jīng)驗公式》一節(jié)會詳細介紹這兩個公式，不過這兩個公式的結果也只能精確到日，不能提供10秒以內精度的時間。要想精確地獲得幾千年乃至更長時間范圍內任意一年的節(jié)氣發(fā)生時間和日月合朔時間，就只能采用“天文算法”。

所謂的“天文算法”，就是利用經(jīng)典力學定律推導行星運轉軌道，對任意時刻的行星位置進行精確計算，從而獲得某種天文現(xiàn)象發(fā)生時的時間，比如日月合朔這一天文現(xiàn)象就是太陽和月亮的地心黃經(jīng)（視黃經(jīng)）差為0的那一瞬間。能夠計算任意時刻行星位置的一套理論就被稱為星歷表，比較著名的星歷表有美國國家航空航天局下屬的噴氣推進實驗室發(fā)布的DE系列星歷表，還有瑞士天文臺在DE406基礎上拓展的瑞士星歷表等等。根據(jù)行星運行軌道直接計算行星位置通常不是很方便，更何況大多數(shù)民用天文計算用不上那么多精確的軌道參數(shù)，于是天文學家在這些星歷表的基礎上推導出了很多可以做簡便計算，但是又能保證一定精度的行星運行理論，比較著名的有VSOP82/87太陽系行星運行理論和ELP-2000/82月球運行理論，這兩套理論在精度上已經(jīng)很接近DE系列星歷表了。關于如何應用這兩套倫理進行天文歷法計算，請參考“日歷生成算法”系列文章的第三篇《用天文方法計算二十四節(jié)氣》和第四篇《用天文方法計算日月合朔》，本文介紹的農(nóng)歷年歷推算是在已經(jīng)通過天文算法獲得了精確的節(jié)氣時間和日月合朔時間的基礎上進行的。

中國的官方紀時采用的是中國公歷（格里歷），因此農(nóng)歷年歷的推導應以公歷年的周期為主導，附上農(nóng)歷年的信息，也就是說，年歷以公歷的1月1日為起始，至12月31日結束，根據(jù)農(nóng)歷歷法推導出的農(nóng)歷日期信息，附加在公歷日期信息上形成雙歷。通常情況下，一個公歷年周期都不能完整地對應到一個農(nóng)歷年周期上，二者的偏差也不固定，因此不存在穩(wěn)定的對應關系，也就是說，不存在從公歷的日期到農(nóng)歷日期的轉換公式，只能根據(jù)農(nóng)歷的歷法規(guī)則推導出農(nóng)歷日期與公歷日期的對應關系。由農(nóng)歷歷法規(guī)則可知，上一個公歷年的冬至（）所在的朔望月是上一個農(nóng)歷年的十一月（冬月），所以在進行節(jié)氣計算時，需要計算包括上一年冬至節(jié)氣在內的二十五個節(jié)氣，才能對應上上一個農(nóng)歷年的十一月和當前農(nóng)歷年的十一月。在計算與之對應的朔日時，考慮到有閏月的情況，需要從上一年冬至節(jié)氣前的第一個朔日，連續(xù)計算15個朔日才能保證覆蓋兩個冬至之間的一整年時間，圖（1）顯示了2011年沒有閏月的情況下朔日和冬至的關系：

圖（1）沒有閏月情況下朔日與冬至節(jié)氣關系圖

圖中上排數(shù)字是公歷月的編號，黑色圓點代表朔日，黑色三角形代表冬至節(jié)氣。圖（2）顯示了2012年有閏月的情況下朔日和冬至的關系：

圖（2）有閏月情況下朔日與冬至節(jié)氣關系圖

通過計算得到能夠覆蓋兩個冬至節(jié)氣的所有朔日時間后，就可以著手建立公歷日期與農(nóng)歷日期的對應關系。以圖（1）所示的2011年為例，首先根據(jù)計算得到的15個朔日（2011年只會用到其中的前14個時間）時間，建立與2011年（公歷年）有關的朔望月關系表：

朔日編號	合朔時間	對應公歷日期	月長	月名
1	01:35:39.90	2010-12-06	29	冬月
2	17:02:34.26	2011-01-04	30	臘月
3	10:30:42.67	2011-02-03	30	正月
4	04:45:59.44	2011-03-05	29	二月
5	22:32:15.13	2011-04-03	30	三月
6	14:50:31.79	2011-05-03	30	四月
7	05:02:32.51	2011-06-02	29	五月
8	16:53:54.10	2011-07-01	30	六月
9	02:39:45.06	2011-07-31	29	七月
10	11:04:06.43	2011-08-29	29	八月
11	19:08:50.09	2011-09-27	30	九月
12	03:55:54.64	2011-10-27	29	十月
13	14:09:40.97	2011-11-25	30	冬月
14	02:06:27.05	2011-12-25	29	臘月
15	15:39:23.99	2012-01-23	30	正月

表（2）2011年朔望月與公歷日期關系表

編號為1和2的兩個朔日之間的朔望月是十一月，因為冬至節(jié)氣落在這個朔望月，其它月的月名依次類推，正月的朔日就是春節(jié)。輸出公歷和農(nóng)歷雙歷時，以月（公歷）為單位，從每月第一天開始，依次判斷每一天屬于哪個朔望月，確定這一天的農(nóng)歷月名，然后比較這一天和這個朔望月的朔日之間相差幾天，記為農(nóng)歷日期。以2011年1月1日為例，這一天在2010年12月6日（2010年農(nóng)歷十一月的朔日）和2011年1月4日之間（2010年農(nóng)歷十二月的朔日），查表（1）可知對應的農(nóng)歷月是十一月，這一天和2010年12月6日相差26天，因此這一天的農(nóng)歷日期就是“廿七”。再以2011年2月3日（春節(jié)）這一天為例，查朔望月表得知2月3日屬于從2月3日開始的朔望月，這個朔望月的月名是正月，而2月3日就是月首，農(nóng)歷日期是初一，正月初一就是春節(jié)。

先來介紹兩個函數(shù)，這兩個函數(shù)分別用于計算節(jié)氣和日月合朔發(fā)生的時間，函數(shù)算法的具體描述將在“日歷生成算法”系列文章的第三篇《用天文方法計算二十四節(jié)氣》和第四篇《用天文方法計算日月合朔》中介紹，此處只是簡單介紹一下用法。首先是計算節(jié)氣時間的函數(shù)：

5 double CalculateSolarTerms(int year, int angle);

這個函數(shù)用于計算指定的年份（year參數(shù)）中，太陽在黃道上運行（視運動）到指定角度時的時間，angle可以設定節(jié)氣發(fā)生時的角度，比如CalculateSolarTerms(2011, 270)就是計算2011年冬至的時間。這個函數(shù)返回的時間類型是儒略日，關于儒略日的說明請參考“日歷生成算法”系列文章的第一篇《中國公歷（格里歷）》。

接下來介紹計算日月合朔時間的函數(shù)：

8 double CalculateMoonShuoJD(double tdJD);

這個函數(shù)返回指定時間附近的朔日時間，搜索的范圍是tdJD參數(shù)指定時間的前一天到后29.5305天，tdJD參數(shù)和返回值的時間類型都是儒略日。

生成指定公歷年份的公歷和農(nóng)歷的雙歷年歷的流程如下：

圖（3）計算公農(nóng)歷雙歷年歷的流程

GetAllSolarTermsJD()函數(shù)從指定年份的指定節(jié)氣開始，連續(xù)計算25個節(jié)氣時間，時間可以跨年份，內部判斷過冬至節(jié)氣后自動轉到下一年的節(jié)氣繼續(xù)計算：

139 void CChineseCalendar::GetAllSolarTermsJD(int year, int start, double*SolarTerms)

140 {

141 int i = 0;

142 int st = start;

143 while(i < 25)

144 {

145 double jd = CalculateSolarTerms(year, st * 15);

147 if(st == WINTER_SOLSTICE)

148 {

149 year++;

150 }

151 st = (st + 1) % SOLAR_TERMS_COUNT;

152 }

153 }

start參數(shù)是節(jié)氣的索引，定義二十四節(jié)氣的索引如下：

38 const int VERNAL_EQUINOX = 0; // 春分

39 const int CLEAR_AND_BRIGHT = 1; // 清明

40 const int GRAIN_RAIN = 2; // 谷雨

41 const int SUMMER_BEGINS = 3; // 立夏

42 const int GRAIN_BUDS = 4; // 小滿

43 const int GRAIN_IN_EAR = 5; // 芒種

44 const int SUMMER_SOLSTICE = 6; // 夏至

45 const int SLIGHT_HEAT = 7; // 小暑

46 const int GREAT_HEAT = 8; // 大暑

47 const int AUTUMN_BEGINS = 9; // 立秋

48 const int STOPPING_THE_HEAT = 10; // 處暑

49 const int WHITE_DEWS = 11; // 白露

50 const int AUTUMN_EQUINOX = 12; // 秋分

51 const int COLD_DEWS = 13; // 寒露

52 const int HOAR_FROST_FALLS = 14; // 霜降

53 const int WINTER_BEGINS = 15; // 立冬

54 const int LIGHT_SNOW = 16; // 小雪

55 const int HEAVY_SNOW = 17; // 大雪

56 const int WINTER_SOLSTICE = 18; // 冬至

57 const int SLIGHT_COLD = 19; // 小寒

58 const int GREAT_COLD = 20; // 大寒

59 const int SPRING_BEGINS = 21; // 立春

60 const int THE_RAINS = 22; // 雨水

61 const int INSECTS_AWAKEN = 23; // 驚蟄

節(jié)氣索引乘以15就是節(jié)氣在黃道上對應的度數(shù)。GetNewMoonJDs()函數(shù)從指定時間開始連續(xù)計算15個朔日時間，從第一個冬至節(jié)氣前的第一個朔日開始。15個朔日可以形成14個完整的朔望月，保證在有閏月的情況下也能包含兩個冬至節(jié)氣：

137 void CChineseCalendar::GetNewMoonJDs(double jd, double *NewMoon)

138 {

139 for(int i = 0; i < NEW_MOON_CALC_COUNT; i++)

140 {

141 double shuoJD = CalculateMoonShuoJD(jd);

142 NewMoon[i] = shuoJD;

143

144 jd += 29.5; /*轉到下一個最接近朔日的時間*/

145 }

146 }

BuildAllChnMonthInfo()函數(shù)根據(jù)15個朔日時間組成14個朔望月，根據(jù)相鄰朔日的間隔計算出農(nóng)歷月天數(shù)用來判定大小月，并且從“十一月”開始依次為每個朔望月命名（月建名稱）：

170 bool CChineseCalendar::BuildAllChnMonthInfo()

171 {

172 CHN_MONTH_INFO info; //一年最多可13個農(nóng)歷月

173 int i;

174 int yuejian = 11; //采用夏歷建寅，冬至所在月份為農(nóng)歷11月

175 for(i = 0; i < (NEW_MOON_CALC_COUNT - 1); i++)

176 {

177 info.mmonth = i;

178 info.mname = (yuejian <= 12) ? yuejian : yuejian - 12;

179 info.shuoJD = m_NewMoonJD[i];

180 info.nextJD = m_NewMoonJD[i + 1];

181 info.mdays = int(info.nextJD + 0.5) - int(info.shuoJD + 0.5);

182 info.leap = 0;

183

184 CChnMonthInfo cm(&info);

185 m_ChnMonthInfo.push_back(cm);

186

187 yuejian++;

188 }

189

190 return (m_ChnMonthInfo.size() == (NEW_MOON_CALC_COUNT - 1));

191 }

CalcLeapChnMonth()函數(shù)根據(jù)節(jié)氣和朔日時間判斷在兩個冬至節(jié)氣之間的農(nóng)歷年是否有閏月，判斷的依據(jù)就是看第十四個朔日是否在第二個冬至節(jié)氣之前，如果第十四個朔日發(fā)生在第二個冬至節(jié)氣之前，就說明在兩個冬至節(jié)氣之間發(fā)生了十三次朔日，需要置閏月。因為農(nóng)歷中十二個中氣屬于哪個農(nóng)歷月是固定的，因此置閏月的過程就是依次判斷十二個中氣是否在對應的農(nóng)歷月中，如果本應該屬于某個農(nóng)歷月的中氣卻沒有落在這個農(nóng)歷月中，則這個農(nóng)歷月就是閏月，需要設置閏月標志，同時調整這個月之后的月名。調整農(nóng)歷月名的方法就是月名減一，比如原來是八月就要調整為七月，這樣就將十三個月對應上了十二個月名（其中多出來的一個農(nóng)歷月被命名為閏某月）。如果節(jié)氣和朔日發(fā)生在同一天，CalcLeapChnMonth()函數(shù)采用的是民間歷法的規(guī)則，與現(xiàn)行歷法一致：

194 void CChineseCalendar::CalcLeapChnMonth()

195 {

196 assert(m_ChnMonthInfo.size() > 0); /*陰歷月的初始化必須在這個之前*/

197

198 int i;

199

200 if(int(m_NewMoonJD[13] + 0.5) <= int(m_SolarTermsJD[24] + 0.5)) //第13月的月末沒有超過冬至，說明今年需要閏一個月

201 {

202 //找到第一個沒有中氣的月

203 i = 1;

204 while(i < (NEW_MOON_CALC_COUNT - 1))

205 {

206

207 /*m_NewMoonJD[i + 1]是第i農(nóng)歷月的下一個月的月首，本該屬于第i月的中氣如果比下一個月

208 的月首還晚，或者與下個月的月首是同一天（民間歷法），則說明第i月沒有中氣*/

209 if(int(m_NewMoonJD[i + 1] + 0.5) <= int(m_SolarTermsJD[2 * i] +0.5))

210 break;

211 i++;

212 }

213 if(i < (NEW_MOON_CALC_COUNT - 1)) /*找到閏月，對后面的農(nóng)歷月調整月名*/

214 {

215 m_ChnMonthInfo[i].SetLeapMonth(true);

216 while(i < (NEW_MOON_CALC_COUNT - 1))

217 {

218 m_ChnMonthInfo[i++].ReIndexMonthName();

219 }

220 }

221 }

222 }

從理論上講，本文介紹的算法在精度允許的范圍內可以計算前后幾千年的農(nóng)歷年歷，但是對古代的農(nóng)歷計算需要小心。首先是“平朔”和“定朔”的問題，唐代以前使用的是平朔方法定月首，本文介紹的計算方法采用的是“定朔”方法，因此計算出的年歷與唐代以前的歷史會不一致。另外，即是在唐代以后采用“定朔”的歷法，因為古代天文觀測和計算受條件限制，可能不夠精確，因此與現(xiàn)在用天文算法計算出的結果可能并不一致。所以對歷史農(nóng)歷的計算應該以歷史事實為主，天文計算為輔，當計算與歷史不一致時，要根據(jù)歷史數(shù)據(jù)進行校正。Calendar.exe是根據(jù)本文介紹的算法編寫的日歷小程序，沒有太多的功能，主要是為了驗證算法，因為沒有歷史數(shù)據(jù)用于修正結果，因此不支持1601年以前的農(nóng)歷計算（也就是說按照天文算法計算出來的結果可能和實際歷史上的歷法不符）。

圖（5）演示程序的界面

小知識1：民間歷法和歷理歷法

新中國成立以后沒有頒布新的“官方農(nóng)歷歷法”，將歷法和政治分離體現(xiàn)了時代的進步，但是由于沒有 “官方歷法”，也引起了一些問題。比如我國現(xiàn)在采用的農(nóng)歷歷法是《時憲歷》，它源于清朝順治年間（公元1645）頒布的《順治歷》，它有兩個不足之處：一個是日月合朔和節(jié)氣的時間以北京當?shù)貢r間為準，也就是東經(jīng)116度25分的當?shù)貢r間，其節(jié)氣和新月的觀察只適用于中原地區(qū)。其它經(jīng)度的地方，因為時間的關系，對導致日月合朔和節(jié)氣時間的差異導致置閏和月順序各不相同。另一個不足之處就是日月合朔時間和節(jié)氣時間判斷不精確，如果日月合朔時間和節(jié)氣時間在同一天，不管具體的時間是否有先后，一律將此節(jié)氣算做新月中的節(jié)氣，這樣一來，如果這個節(jié)氣是中氣，就會影響到閏月的設置。歷理歷法針對這兩點進行了改進，對節(jié)氣時間和日月合朔時間統(tǒng)一采用東經(jīng)120度即東八區(qū)標準時，這樣在任何時區(qū)的節(jié)氣和置閏結果都是一樣的，以東八區(qū)標準時為準。對于節(jié)氣時間和日月合朔時間在同一天的情況，精確計算到時、分、秒，只有日月合朔時間在節(jié)氣時間之前，這個節(jié)氣才包含在次月內。歷理歷法從理論上講更符合現(xiàn)代天文學的精確計算，但是需要注意的是，歷理歷法仍然只是存在于理論上的歷法，我國現(xiàn)行的農(nóng)歷歷法依然是民間歷法《時憲歷》或《順治歷》。

小知識2：通式壽星公式

“通式壽星公式”是前人整理出來的一個用于計算每年立春日期的經(jīng)驗公式:

Date = 向下取整(Y * D + C) - L

其中，Y是年份，D的值是0.2422，C是經(jīng)驗值，取決于節(jié)氣和年份，對于21世紀，立春節(jié)氣的C值是4.475，春分節(jié)氣的C值是20.646等等；

L是閏年數(shù)，其計算公式為：

L = 向下取整(Y/4) - 向下取整(Y/100) + 向下取整(Y/400)

用“通式壽星公式”確定2011年立春日期的過程如下：

L = int(2011/4) – int(2011/100) + int(2011/400) = 502 – 20 + 5 = 487

Date = int（2011×0.2422+4.475）- 487 = 491 – 487 = 4

所以，2011年的立春日期是2月4日。

小知識3：計算節(jié)氣和朔日的經(jīng)驗公式

以1900年1月0日（星期日）為基準日，之后的每一天與基準日的差值稱為“積日”， 1900年1月1日的積日是1，以后的時間依次類推，則計算第y年第x個節(jié)氣的積日公式是：

F = 365.242 * (y – 1900) + 6.2 + 15.22 *x - 1.9 * sin(0.262 * x)

其中x是節(jié)氣的索引，0代表小寒，1代表大寒，其它節(jié)氣按照順序類推。

計算從1900年開始第m個朔日的公式是：

M = 1.6 + 29.5306 * m + 0.4 * sin(1 - 0.45058 * m)

小知識4：平朔和定朔

中國農(nóng)歷的朔望月長度是平均29.5305天，所以農(nóng)歷月就有大月30天，小月29天之分，從先秦時期到唐代，農(nóng)歷歷法均是采用大小月輪流交替的方式設置每個農(nóng)歷月的天數(shù)，只有少數(shù)情況下才出現(xiàn)連續(xù)兩個大月的情況，采用這種方式的歷法就稱為“平朔”。“平朔”歷法簡單，但是不能保證日月合朔發(fā)生在初一這一天，有可能是上月的月末一天，也有可能是本月初二。南北朝時期，一種新的歷法被提出來，這種歷法嚴格按照日月合朔為月初制定農(nóng)歷月，采用這種方式的歷法就稱為“定朔”。“定朔”歷法嚴格將日月合朔時間確定月初，因為月球公轉是橢圓軌道，速度并不是均勻，所以會發(fā)生連續(xù)多個大月或連續(xù)多個小月的情況，導致“定朔”歷法推廣遇到很大的阻力，直到唐代，中國歷法才全面棄用“平朔”，改用“定朔”。

小知識5：正月初一和立春節(jié)氣

立春是二十四節(jié)氣之首，所以古代民間都是在“立春”這一天過節(jié)，相當于現(xiàn)代的春節(jié)（中國古代即是節(jié)氣也是節(jié)日的情況很多，比如清明、冬至等等）。1911年，孫中山領導的辛亥革命建立了中華民國，在從歷法上正式把農(nóng)歷正月初一定為“春節(jié)”，把公歷1月1日定為“元旦”，也就是“新年”。農(nóng)歷年從正月初一開始沒有爭議，但是農(nóng)歷生肖年從何時開始卻一直有爭議，目前多數(shù)人都認為“立春”節(jié)氣是農(nóng)歷生肖年的開始。因為在中國古代歷法中，十二生肖的計算與天干地支有很大關系，所以在“論天干地支、計算廿四節(jié)氣”的情況下，“立春”節(jié)氣應該是新生肖的開始。對于普通老百姓來說，習慣于認為正月初一是生肖年的開始，因此，正月初一和“立春”節(jié)氣之間出生的小孩，在確定屬相的時候就有點麻煩了。屬龍還是屬蛇？這是個問題。

小知識1：公歷的閏年

中國公歷（也就是格里歷）的置閏規(guī)則是四年一閏，百年不閏，四百年再閏，為什么會有這么奇怪的置閏規(guī)則呢？這實際上與天體運行周期與人類定義的歷法周期之間的誤差有關。地球繞太陽運轉的周期是365.2422天，即一個回歸年（Tropical Year），而公歷的一年是365天，這樣一年就比回歸年短了0.2422日，四年積累下來就多出0.9688天（約1天），于是設置一個閏年，這一年多一天。這樣一來，四個公歷年又比四個回歸年多了0.0312天，平均每年多0.0078天，這樣經(jīng)過四百年就會多出3.12天，也就是說每四百年要減少3個閏年才行，于是就設置了百年不閏，四百年再閏的置閏規(guī)則。

實際上公歷的置閏還有一條規(guī)則，就是對于數(shù)值很大的年份，如果能整除3200，同時能整除172800則是閏年。這是因為前面即使四百年一閏，仍然多了0.12天，平均就是每天多0.0003天，于是每3200年就又多出0.96天，也就是說每3200年還要減少一個閏年，于是能被3200整除的年就不是閏年了。然而誤差并沒有終結，每3200年減少一個閏年（減少一天）實際上多減了0.04天，這個誤差還要繼續(xù)累計計算，這已經(jīng)超出了本文的范圍，有興趣的讀者可以自己計算。

小知識2：儒略歷和格里歷

在公元1582年10月15日之前，人們使用的歷法是源自古羅馬的儒略歷，儒略歷的置閏規(guī)則就是四年一閏，但是沒有計算每年多出來的0.0078天，這樣從公元前46年到公元1582年一共累積多出了10天，為此，當時的教皇格里十三世將1582年10月5日人為指定為10月15日，并開始啟用新的置閏規(guī)則，這就是后來沿用至今的格里歷。

小知識3：約化儒略日

由于儒略日數(shù)字位數(shù)太多，國際天文聯(lián)合會于1973年8月決定對其修正，采用約化儒略日（MJD）進行天文計算，定義MJD = JD – 2400000.5，MJD相應的起始點是1858年11月17日 0:00。

小知識4：1752年9月到底是怎么回事兒

如果你用的操作系統(tǒng)是unix或linux，在控制臺輸入以下命令：

#cal 9 1752

你會看到這樣一個奇怪的月歷輸出：

September 1752

Su Mo Tu We Th Fr Sa

1 2 14 15 16

17 18 19 20 21 22 23

24 25 26 27 28 29 30

1752年的9月缺了11天，到底怎么回事兒？這其實還是因為從儒略歷到格里歷的轉換造成的。1582年10月5日，羅馬教皇格里十三世宣布啟用更為精確的格里歷，但是整個歐洲大陸并不是所有國家都立即采用格里歷，比如大英帝國就是直到1752年9月議會才批準采用格里歷，所以大英帝國及其所有殖民地的歷法一直到1752年9月才發(fā)生跳變，“跟上”了格里歷。德國和荷蘭到了1698年才采用格里歷，而俄羅斯則直到1918年革命才采用格里歷。Linux的cal指令起源與最初AT&T的UNIX，當然采用的是美國歷法，但是美國歷史太短，再往前就只能采用英國歷法，所以cal指令的結果就成了這樣。對于采用格里歷的國家來說，只要知道1582年10月發(fā)生了日期跳變就行了，可以不用關心1752年9月到底是怎么回事兒。但是對于研究歷史和考古的人來說，就必需要了解這個歷史，搞清楚每個歐洲國家改用格里歷的年份，否則就可能在一些問題上出錯。在歐洲研究歷史，你會發(fā)現(xiàn)很多事件都是有多個時間版本的，比如大科學家牛頓的生日就有兩個時間版本，一個是按照儒略歷歷法的1642年12月25日，另一個是格里歷歷法的1643年1月4日，對于英國人來說，1752年之前都是按照儒略歷計算的，所以英國的史書可能會記載牛頓出生在圣誕節(jié)，這也沒什么可奇怪的。

居天下之廣居，立天下之正位，行天下之大道，得志與民由之，不得志獨行其道，富貴不能淫，貧賤不能移，威武不能屈，此之謂大丈夫。

--noryes