全文引用自：

https://www.cnblogs.com/noryes/p/5716673.html

朔方

算法系列：日歷算法

Copyright © 1900-2016, NORYES, All Rights Reserved.

http://www.cnblogs.com/noryes/

歡迎轉(zhuǎn)載，請保留此版權(quán)聲明。

---------------------------------------------------------------------------------------

    轉(zhuǎn)載自http://blog.csdn.net/orbit/article/details/7749723

1、概述

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

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

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

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

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

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

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

 

2、計(jì)算星期

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

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

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

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) /*間隔超過一年，要計(jì)算間隔的整年時(shí)間*/

                            {

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

                               }

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

                    }

            else

                {

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

                     }

                return days;

        }

 /*計(jì)算一年中過去的天數(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;

        }

    /*計(jì)算一年中還剩下的天數(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;

        }

     /*

      計(jì)算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;

}

    另一種計(jì)算兩個(gè)日期相差天數(shù)的方法是利用儒略日（Julian Day，JD）進(jìn)行計(jì)算。首先介紹一下儒略日，儒略日是一種不記年，不記月，只記日的歷法，是由法國學(xué)者 Joseph Justus Scaliger（1540－1609）在 1583 年提出來的一種以天數(shù)為計(jì)量單位的流水日歷。儒略日和儒略歷（Julian Calendar）沒有任何關(guān)系，命名為儒略日也僅僅他本人為了紀(jì)念他的父親--意大利學(xué)者 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í)，依次順推，每一天都被賦予一個(gè)唯一的數(shù)字。例如從 1996 年 1 月 1 日 12:00 開始的一天就是儒略日J(rèn)D2450084。使用儒略日可以把不同歷法的年表統(tǒng)一起來，很方便地在各種歷法中追溯日期。如果計(jì)算兩個(gè)日期之間的天數(shù)，利用儒略日計(jì)算也很方便，先計(jì)算出兩個(gè)日期的儒略日數(shù)，然后直接相減就可以得到兩個(gè)日期相隔的天數(shù)。

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

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

    下面就是由公歷日期計(jì)算儒略日的算法實(shí)現(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 }

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

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

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

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

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

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

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

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

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

得到 2005 年 5 月 31 日是星期二，和前面的計(jì)算方法得到的結(jié)果一致。根據(jù)上述分析，可得寫出使用公式 5 計(jì)算星期的算法實(shí)現(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 的問題在于計(jì)算量大，不利于口算星期結(jié)果。于是人們就在公式 5 的基礎(chǔ)上繼續(xù)推導(dǎo)更簡單的公式。德國數(shù)學(xué)家克里斯蒂安·蔡勒（Christian Zeller, 1822- 1899）在 1886 年推導(dǎo)出了著名的蔡勒（Zeller）公式：

 

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

w ：星期；

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

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

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

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

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

    目前人們普遍認(rèn)為蔡勒公式是計(jì)算某一天是星期幾的最好的公式。但是蔡勒公式有時(shí)候可能計(jì)算出的結(jié)果是負(fù)數(shù)，需要對結(jié)果+7進(jìn)行修正。比如2006年7月1日，用蔡勒公式計(jì)算出的結(jié)果是 -1，實(shí)際上這天是星期六。根據(jù)前面分析的結(jié)果整理出的蔡勒公式算法實(shí)現(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的月份進(jìn)行修正*/ 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) /*修正計(jì)算結(jié)果是負(fù)數(shù)的情況*/ 171         w += 7; 172  173     return w; 174 }

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

 

公式 7 適用于對 1582 年 10 月 4 日之前的日期計(jì)算星期，1582 年 10 月 5 日與 1582 年 10 月 15 日之間的日期是不存在的，因?yàn)樗鼈兌际峭惶臁?/p>

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

 

229 void PrintMonthCalendar(int year, int month) 230 { 231     int days = GetDaysOfMonth(year, month); /*確定這個(gè)月的天數(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) /*到一周結(jié)束，切換到下一行輸出*/ 245         { 246             SetNextRowStart(); 247         } 248         i++; 249         week = (week + 1) % 7; 250     } 251 }

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

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

計(jì)算星期不必對每一天都計(jì)算一次，只要對每個(gè)月的第一天計(jì)算一次就可以了，以后的日期可以用 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é)氣的基本知識(shí)，以及如何使用 VSOP82/87 行星運(yùn)行理論計(jì)算二十四節(jié)氣發(fā)生的準(zhǔn)確時(shí)間。

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

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

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

圖（1）天球概念示意圖

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

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

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

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

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

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

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

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

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

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

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

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

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

JDE是計(jì)算軌道參數(shù)的時(shí)間，單位是儒略日，2451545.0是公元2000年1月1日 12時(shí)的儒略日數(shù)，關(guān)于儒略日的概念，請參考“日歷生成算法”的第一篇《中國公歷（格里歷）》中的說明以及計(jì)算方法。以L0表的第二個(gè)周期項(xiàng)為例，這個(gè)周期項(xiàng)數(shù)據(jù)中與日心黃經(jīng)計(jì)算有關(guān)的三個(gè)系數(shù)分別是A= 3341656.456，B=4.66925680417，C=6283.07584999140，則第二個(gè)周期項(xiàng)的計(jì)算方法是：3341656.456 * cos(4.66925680417 + 6283.0758499914 * τ)。對L0表的各項(xiàng)分別計(jì)算后求和可得到L0表周期項(xiàng)總和L0，對L表的其它幾個(gè)部分使用相同的方法計(jì)算周期項(xiàng)和，可以得到L1、L2、L3、L4和L5，然后用用3.5式計(jì)算出最終的地球日心黃經(jīng)，單位是弧度：

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

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

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

本文就是使用3.6式代替3.5式進(jìn)行計(jì)算。

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

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

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

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

3.7式中的T是儒略世紀(jì)數(shù)，它與儒略千年數(shù)τ的關(guān)系是：T = 10 *τ。然后使用L'計(jì)算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式計(jì)算出得L和B轉(zhuǎn)換成以度（°）為單位的值后再進(jìn)行修正。

        CalcSunEclipticLongitudeEC()函數(shù)就是使用VSOP87行星理論計(jì)算行星日心黃經(jīng)的代碼實(shí)現(xiàn)，整個(gè)計(jì)算過程和前文描述一樣，首先根據(jù)VSOP87D表的數(shù)據(jù)計(jì)算出L0-L5，然后用3.6式計(jì)算出地球的日心黃經(jīng)，3.6式計(jì)算出來的單位是弧度，因此轉(zhuǎn)換成度分秒單位，最后使用3.1式將結(jié)果轉(zhuǎn)換成太陽的地心黃經(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°的值調(diào)整到0－360°之間，便于轉(zhuǎn)換顯示。CalcPeriodicTerm()函數(shù)使用3.3式對一個(gè)周期項(xiàng)系數(shù)表進(jìn)行求和計(jì)算，可以指定需要計(jì)算的周期項(xiàng)數(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;  }

同樣的方法計(jì)算太陽的地心黃緯：

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

double AdjustSunEclipticLongitudeEC(double dt, double longitude, double latitude)  {      double T = dt * 10; //T是儒略世紀(jì)數(shù)      double dbLdash = longitude - 1.397 * T - 0.00031 * T * T;      // 轉(zhuǎn)換為弧度      dbLdash *= RADIAN_PER_ANGLE;      return (-0.09033 + 0.03916 * (cos(dbLdash) + sin(dbLdash)) * tan(latitude *RADIAN_PER_ANGLE)) / 3600.0;  }

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

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

小知識(shí)1：天文單位

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

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

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

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

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

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

平距角(日月對地心的角距離)：
D = 297.85036 + 455267.111480 * T - 0.0019142 * T² + T³ / 189474        （3.10式）
太陽（地球）平近點(diǎn)角：
M = 357.52772 + 35999.050340 * T - 0.0001603 * T² - T³ / 300000         （3.11式）
月球平近點(diǎn)角
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式）
黃道與月球平軌道升交點(diǎn)黃經(jīng)：
Ω= 125.04452 - 1934.136261 * T + 0.0020708 * T² + T³ / 450000            （3.14式）

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

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

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

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

CalcEarthLongitudeNutation()函數(shù)就是計(jì)算黃經(jīng)章動(dòng)的實(shí)現(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);      }      /*先乘以章動(dòng)表的系數(shù) 0.0001，然后換算成度的單位*/      return resulte * 0.0001 / 3600.0;  }

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

void GetEarthNutationParameter(double dt, double *D, double *M, double *Mp, double*F, double *Omega)  {      double T = dt * 10; /*T是從J2000起算的儒略世紀(jì)數(shù)*/      double T2 = T * T;      double T3 = T2 * T;      /*平距角（如月對地心的角距離）*/      *D = 297.85036 + 445267.111480 * T - 0.0019142 * T2 + T3 / 189474.0;      /*太陽（地球）平近點(diǎn)角*/      *M = 357.52772 + 35999.050340 * T - 0.0001603 * T2 - T3 / 300000.0;      /*月亮平近點(diǎn)角*/      *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;      /*黃道與月亮平軌道升交點(diǎn)黃經(jīng)*/      *Omega = 125.04452 - 1934.136261 * T + 0.0020708 * T2 + T3 / 450000.0;  }

 

同樣，計(jì)算交角章動(dòng)的實(shí)現(xiàn)代碼是：

double CalcEarthObliquityNutation(double dt)  {      double T = dt * 10; /*T是從J2000起算的儒略世紀(jì)數(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);      }      /*先乘以章動(dòng)表的系數(shù) 0.001，然后換算成度的單位*/      return resulte * 0.0001 / 3600.0;  }

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

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

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

分子20.4898并不是一個(gè)常數(shù)，但是其只的變化非常緩慢，在0年是20".4893，在4000年是20".4904。前文提到過，太陽到地球的距離R可以用VSOP87D表的R0－R5周期項(xiàng)計(jì)算出來，R的單位是“天文單位（AU）”，和計(jì)算太陽地心黃經(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，而用下面的公式直接計(jì)算日地距離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的計(jì)算公式是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式計(jì)算時(shí)需要轉(zhuǎn)換成弧度單位，計(jì)算出R以后，就可以這樣計(jì)算光行差修正量：

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

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

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

 

double GetSunEclipticLongitudeEC(double jde)  {      double dt = (jde - JD2000) / 365250.0; /*儒略千年數(shù)*/      // 計(jì)算太陽的地心黃經(jīng)      double longitude = CalcSunEclipticLongitudeEC(dt);      // 計(jì)算太陽的地心黃緯      double latitude = CalcSunEclipticLatitudeEC(dt) * 3600.0;      // 修正精度      longitude += AdjustSunEclipticLongitudeEC(dt, longitude, latitude);      // 修正天體章動(dòng)      longitude += CalcEarthLongitudeNutation(dt);      // 修正光行差      /*太陽地心黃經(jīng)光行差修正項(xiàng)是: -20".4898/R*/      longitude -= (20.4898 / CalcSunEarthRadius(dt)) / (20 * PI);      return longitude;  }

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

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

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

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

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

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

34 double CalculateSolarTerms(int year, int angle) 35 { 36     double lJD, rJD; 37     EstimateSTtimeScope(year, angle, lJD, rJD); /*估算迭代起始時(shí)間區(qū)間*/ 38  39     double solarTermsJD = 0.0; 40     double longitude = 0.0; 41  42     do 43     { 44         solarTermsJD = ((rJD - lJD) * 0.618) + lJD; 45         longitude = GetSunEclipticLongitudeECDegree(solarTermsJD); 46         /* 47             對黃經(jīng)0度迭代逼近時(shí)，由于角度360度圓周性，估算黃經(jīng)值可能在(345,360]和[0,15)兩個(gè)區(qū)間， 48             如果值落入前一個(gè)區(qū)間，需要進(jìn)行修正 49         */ 50         longitude = ((angle == 0) && (longitude > 345.0)) ? longitude - 360.0 :longitude; 51  52         (longitude > double(angle)) ? rJD = solarTermsJD : lJD = solarTermsJD; 53     }while((rJD - lJD) > 0.0000001); 54  55     return solarTermsJD; 56 }

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

        二分逼近或黃金分割逼近算法實(shí)現(xiàn)簡單，很容易控制，但是也存在效率低，收斂速度慢的問題，現(xiàn)在我們介紹牛頓迭代法，牛頓迭代法是一種在實(shí)數(shù)域和復(fù)數(shù)域上近似求解方程的方法。假設(shè)我們要求解的方程是f(x) = 0，如果f(x)的導(dǎo)函數(shù)f’(x)是連續(xù)的，則在真實(shí)解x附近的區(qū)域內(nèi)任意一點(diǎn)x₀開始迭代，則牛頓迭代法必收斂，特別當(dāng)f’(x)不等于0的時(shí)候，牛頓迭代法是平方收斂的，也就是說，每迭代一次，結(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ù)。因?yàn)镚etSunEclipticLongitudeEC()函數(shù)內(nèi)部是調(diào)用其他函數(shù)，因此可以理解為是一個(gè)多個(gè)函數(shù)組合的復(fù)合函數(shù)，類似f(x) = g(x) + h(x, k(x)) + p(x)這樣的形式，可以按照求導(dǎo)規(guī)則逐步對其求導(dǎo)得到導(dǎo)函數(shù)。但是我不打算這么做，因?yàn)槲矣懈唵蔚姆椒ǎ蔷褪鞘褂糜?jì)算導(dǎo)數(shù)的近似公式。其實(shí)求導(dǎo)函數(shù)的目的就是為了得到某一點(diǎn)的導(dǎo)數(shù)，如果有近似公式可以直接得到這一點(diǎn)的導(dǎo)數(shù)，就不用費(fèi)勁求導(dǎo)函數(shù)了。

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

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

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

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

74 double CalculateSolarTermsNewton(int year, int angle) 75 { 76     double JD0, JD1,stDegree,stDegreep; 77  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)); 87  88     return JD1; 89 }

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

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

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   冬至

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

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

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

        原子時(shí)是以物質(zhì)的原子內(nèi)部發(fā)射的電磁振蕩頻率為基準(zhǔn)的時(shí)間計(jì)量系統(tǒng)。原子時(shí)的初始?xì)v元規(guī)定為1958年1月1日世界時(shí)0時(shí)，秒長定義為銫-133原子基態(tài)的兩個(gè)超精細(xì)能級(jí)間在零磁場下躍遷輻射9192631770周所持續(xù)的時(shí)間。這是一種均勻的時(shí)間計(jì)量系統(tǒng)。1967年起，原子時(shí)已取代力學(xué)時(shí)作為基本時(shí)間計(jì)量系統(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)歷歷法的基礎(chǔ)，而朔望月又是嚴(yán)格以日月合朔發(fā)生的那一天作為月首，因此日月合朔時(shí)間的計(jì)算是制定農(nóng)歷歷法的關(guān)鍵。本文將介紹ELP-2000/82月球運(yùn)行理論，以及如何用ELP-2000/82月球運(yùn)行理論計(jì)算日月合朔時(shí)間。

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

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

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

        要計(jì)算日月合朔，需要知道太陽地心視黃經(jīng)和月亮地心視黃經(jīng)的計(jì)算方法。“日歷生成算法”系列的第三篇《用天文方法計(jì)算二十四節(jié)氣》一文已經(jīng)介紹了如何用VSOP82/87行星理論計(jì)算太陽的地心視黃經(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)對結(jié)果的增益只有0.00001角秒，還有一些地月距離周期項(xiàng)對距離結(jié)果的增益只有0.02米，對于精度不高的歷法計(jì)算，完全可以忽略。

        有很多基于ELP-2000/82月球理論的改進(jìn)或簡化理論，《天文算法》一書的第四十五章就介紹了一種改進(jìn)算法，其周期項(xiàng)參數(shù)都是從ELP-2000/82理論的周期項(xiàng)參數(shù)轉(zhuǎn)換來的，忽略了小的周期項(xiàng)。使用該方法計(jì)算的月球黃經(jīng)精度只有10”，月亮黃緯精度只有4”，但是只用計(jì)算60個(gè)周期項(xiàng)，速度很快，本文就采用這種修改過的ELP-2000/82理論計(jì)算月亮的地心視黃經(jīng)。這種計(jì)算方法的周期項(xiàng)分三部分，分別用來計(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、太陽平近地角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í)候，需要注意對包含了太陽平近地角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對應(yīng)的系數(shù)是0，因此I₂不需要修正，但是第五個(gè)周期項(xiàng)中M對應(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，對于太陽平近地角M的系數(shù)b不是0的情況也需要乘以E或E²進(jìn)行修正。地月距離的周期項(xiàng)和Σr也可以按照上面的方法計(jì)算，計(jì)算地月距離的目的是為了計(jì)算月亮光行差，因?yàn)榈卦戮嚯x較小，從地球觀察月亮產(chǎn)生的光行差也很小，相對于本文算法的精度（月球黃經(jīng)精度10”，月亮黃緯精度4”）來說，可以忽略光行差修正，因此就不用計(jì)算地月距離。

        由于金星和木星對月球的攝動(dòng)影響，需要對計(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)的代碼：

123 double GetMoonEclipticLongitudeEC(double dbJD) 124 { 125     double Lp,D,M,Mp,F,E; 126     double dt = (dbJD - JD2000) / 36525.0; /*儒略世紀(jì)數(shù)*/ 127  128     GetMoonEclipticParameter(dt, &Lp, &D, &M, &Mp, &F, &E); 129  130     /*計(jì)算月球地心黃經(jīng)周期項(xiàng)*/ 131     double EI = CalcMoonECLongitudePeriodicTbl(D, M, Mp, F, E); 132  133     /*修正金星,木星以及地球扁率攝動(dòng)*/ 134     EI += CalcMoonLongitudePerturbation(dt, Lp, F); 135  136     /*計(jì)算月球地心黃經(jīng)*/ 137     double longitude = Lp + EI / 1000000.0; 138  139     /*計(jì)算天體章動(dòng)干擾*/ 140     longitude += CalcEarthLongitudeNutation(dt / 10.0); 141  142     longitude = Mod360Degree(longitude); 143     return longitude; 144 }

函數(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、太陽平近地角M、月亮平近地角M’、月球生交點(diǎn)平角距F、修正系數(shù)E和月球平黃經(jīng)L’，不需多說明，只要根據(jù)以上各式直接計(jì)算即可。CalcMoonECLongitudePeriodicTbl()函數(shù)計(jì)算60個(gè)月球黃經(jīng)周期項(xiàng)和，并根據(jù)M值系數(shù)的情況進(jìn)行修正，算法實(shí)現(xiàn)如下：

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

87 double CalcMoonLongitudePerturbation(double dt, double Lp, double F) 88 { 89     double T = dt; /*T是'ca?從'b4?J2000起'c6?算'cb?的'b5?儒'c8?略'c2?世'ca?紀(jì)'bc?數(shù)'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 }

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

        關(guān)于牛頓迭代法可以參考相關(guān)的數(shù)學(xué)資料，“日歷生成算法”系列的第三篇《用天文方法計(jì)算二十四節(jié)氣》一文對如何使用牛頓迭代法有簡單的介紹，可以參考一下?？偟膩碚f，就是要先定義需要求解的方程f(x)，根據(jù)上文的介紹，我們需要求解的是太陽的地心黃經(jīng)和月亮的地心黃經(jīng)差值是0的時(shí)候的時(shí)間，《用天文方法計(jì)算二十四節(jié)氣》一文已經(jī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ù)仍然使用近似公式，也不解釋了，直接上迭代求解的代碼了：

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 }

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

        再次說明一下，以上算法中討論的時(shí)間都是力學(xué)時(shí)時(shí)間（TD），與國際協(xié)調(diào)時(shí)（UTC）以及各個(gè)時(shí)區(qū)的本地時(shí)間都有不同，以上計(jì)算出來的時(shí)間都需要調(diào)整成本地時(shí)間，比如中國的中原地區(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

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

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

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

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

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

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

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

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

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]);

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

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

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

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

m  365.2422 = n  29.5306

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

235 – 12  19 = 7

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

19  365.2422 - 235  29.5306  0.0892

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

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

        地球沿著一個(gè)近似橢圓軌道繞太陽公轉(zhuǎn)，這個(gè)公轉(zhuǎn)軌道所在的平面就是“黃道面”，黃道面向外延伸與天球的交線就是“黃道”。古人由于觀測條件限制，只能根據(jù)視覺感覺認(rèn)為是太陽沿著黃道繞地球運(yùn)轉(zhuǎn)，因此設(shè)定太陽從黃經(jīng)（黃道經(jīng)度）零度起（以春分點(diǎn)為起點(diǎn)自西向東度量），將太陽沿黃經(jīng)每運(yùn)行15度所經(jīng)歷的時(shí)日稱為“一個(gè)節(jié)氣”。太陽每年運(yùn)行360度，共經(jīng)歷二十四個(gè)節(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é)氣又細(xì)分為十二節(jié)氣和十二中氣，二十四節(jié)氣按照順序排在奇數(shù)位置上的就是節(jié)氣，排在偶數(shù)位置上的就是中氣。也就是說，立春、驚蟄、清明、立夏、芒種、小暑、立秋、白露、寒露、立冬、大雪和小寒就是十二個(gè)節(jié)氣，而雨水、春分、谷雨、小滿、夏至、大暑、處暑、秋分、霜降、小雪、冬至和大寒就是十二個(gè)中氣。二十四個(gè)節(jié)氣平分在公歷的12個(gè)月中，每月一節(jié)氣一中氣。二十四節(jié)氣反映了太陽的周年運(yùn)動(dòng)（以地球?yàn)閰⒄瘴锏囊曔\(yùn)動(dòng)），所以節(jié)氣在現(xiàn)行的公歷中日期基本固定，上半年在6日、21日，下半年在8日、23日，前后不差 1～2天。中國民間流傳的《二十四節(jié)氣歌》就是為了方便記憶這些節(jié)氣：

春雨驚春清谷天，

夏滿芒夏暑相連，

秋處露秋寒霜降，

冬雪雪冬小大寒，

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

最多相差一兩天。

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

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

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

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

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

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

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

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

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

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

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

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

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

 

節(jié)氣名	與上一節(jié)氣之間的時(shí)間差	與小寒節(jié)氣的累積時(shí)間差
小寒	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é)氣時(shí)間間隔表（單位：秒鐘）

 

已知1900年小寒時(shí)刻為1月6日2:05:00，以這個(gè)節(jié)氣時(shí)刻為基準(zhǔn)，推算其它年份節(jié)氣的算法實(shí)現(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 };    15     16 //已知1900年小寒時(shí)刻為1月6日02:05:00    17 const double base1900_SlightColdJD = 2415025.5868055555;    18     19 double CalculateSolarTermsByExp(int year, int st)    20 {    21     if((st < 0) || (st > 24))    22         return 0.0;    23     24     double stJd = 365.24219878 * (year - 1900) + s_stAccInfo[st] / 86400.0;    25     26     return base1900_SlightColdJD + stJd;    27     28 }

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

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

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

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

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

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

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

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

 

朔日編號(hào)	合朔時(shí)間	對應(yīng)公歷日期	月長	月名
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年朔望月與公歷日期關(guān)系表

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

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

5 double CalculateSolarTerms(int year, int angle);

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

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

 

8 double CalculateMoonShuoJD(double tdJD);

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

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

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

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

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é)氣在黃道上對應(yīng)的度數(shù)。GetNewMoonJDs()函數(shù)從指定時(shí)間開始連續(xù)計(jì)算15個(gè)朔日時(shí)間，從第一個(gè)冬至節(jié)氣前的第一個(gè)朔日開始。15個(gè)朔日可以形成14個(gè)完整的朔望月，保證在有閏月的情況下也能包含兩個(gè)冬至節(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; /*轉(zhuǎn)到下一個(gè)最接近朔日的時(shí)間*/   145     }   146 }

 

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

170 bool CChineseCalendar::BuildAllChnMonthInfo()   171 {   172     CHN_MONTH_INFO info; //一年最多可13個(gè)農(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é)氣和朔日時(shí)間判斷在兩個(gè)冬至節(jié)氣之間的農(nóng)歷年是否有閏月，判斷的依據(jù)就是看第十四個(gè)朔日是否在第二個(gè)冬至節(jié)氣之前，如果第十四個(gè)朔日發(fā)生在第二個(gè)冬至節(jié)氣之前，就說明在兩個(gè)冬至節(jié)氣之間發(fā)生了十三次朔日，需要置閏月。因?yàn)檗r(nóng)歷中十二個(gè)中氣屬于哪個(gè)農(nóng)歷月是固定的，因此置閏月的過程就是依次判斷十二個(gè)中氣是否在對應(yīng)的農(nóng)歷月中，如果本應(yīng)該屬于某個(gè)農(nóng)歷月的中氣卻沒有落在這個(gè)農(nóng)歷月中，則這個(gè)農(nóng)歷月就是閏月，需要設(shè)置閏月標(biāo)志，同時(shí)調(diào)整這個(gè)月之后的月名。調(diào)整農(nóng)歷月名的方法就是月名減一，比如原來是八月就要調(diào)整為七月，這樣就將十三個(gè)月對應(yīng)上了十二個(gè)月名（其中多出來的一個(gè)農(nóng)歷月被命名為閏某月）。如果節(jié)氣和朔日發(fā)生在同一天，CalcLeapChnMonth()函數(shù)采用的是民間歷法的規(guī)則，與現(xiàn)行歷法一致：

194 void CChineseCalendar::CalcLeapChnMonth()   195 {   196     assert(m_ChnMonthInfo.size() > 0); /*陰歷月的初始化必須在這個(gè)之前*/   197    198     int i;   199   200     if(int(m_NewMoonJD[13] + 0.5) <= int(m_SolarTermsJD[24] + 0.5)) //第13月的月末沒有超過冬至，說明今年需要閏一個(gè)月   201     {   202         //找到第一個(gè)沒有中氣的月   203         i = 1;   204         while(i < (NEW_MOON_CALC_COUNT - 1))   205         {   206    207             /*m_NewMoonJD[i + 1]是第i農(nóng)歷月的下一個(gè)月的月首，本該屬于第i月的中氣如果比下一個(gè)月   208               的月首還晚，或者與下個(gè)月的月首是同一天（民間歷法），則說明第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)歷月調(diào)整月名*/   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èi)可以計(jì)算前后幾千年的農(nóng)歷年歷，但是對古代的農(nóng)歷計(jì)算需要小心。首先是“平朔”和“定朔”的問題，唐代以前使用的是平朔方法定月首，本文介紹的計(jì)算方法采用的是“定朔”方法，因此計(jì)算出的年歷與唐代以前的歷史會(huì)不一致。另外，即是在唐代以后采用“定朔”的歷法，因?yàn)楣糯煳挠^測和計(jì)算受條件限制，可能不夠精確，因此與現(xiàn)在用天文算法計(jì)算出的結(jié)果可能并不一致。所以對歷史農(nóng)歷的計(jì)算應(yīng)該以歷史事實(shí)為主，天文計(jì)算為輔，當(dāng)計(jì)算與歷史不一致時(shí)，要根據(jù)歷史數(shù)據(jù)進(jìn)行校正。Calendar.exe是根據(jù)本文介紹的算法編寫的日歷小程序，沒有太多的功能，主要是為了驗(yàn)證算法，因?yàn)闆]有歷史數(shù)據(jù)用于修正結(jié)果，因此不支持1601年以前的農(nóng)歷計(jì)算（也就是說按照天文算法計(jì)算出來的結(jié)果可能和實(shí)際歷史上的歷法不符）。

圖（5）演示程序的界面

小知識(shí)1：民間歷法和歷理歷法

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

小知識(shí)2：通式壽星公式

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

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

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

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

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日。

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

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

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

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

計(jì)算從1900年開始第m個(gè)朔日的公式是：

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

小知識(shí)4：平朔和定朔

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

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

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

小知識(shí)1：公歷的閏年

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

實(shí)際上公歷的置閏還有一條規(guī)則，就是對于數(shù)值很大的年份，如果能整除3200，同時(shí)能整除172800則是閏年。這是因?yàn)榍懊婕词顾陌倌暌婚c，仍然多了0.12天，平均就是每天多0.0003天，于是每3200年就又多出0.96天，也就是說每3200年還要減少一個(gè)閏年，于是能被3200整除的年就不是閏年了。然而誤差并沒有終結(jié)，每3200年減少一個(gè)閏年（減少一天）實(shí)際上多減了0.04天，這個(gè)誤差還要繼續(xù)累計(jì)計(jì)算，這已經(jīng)超出了本文的范圍，有興趣的讀者可以自己計(jì)算。

小知識(shí)2：儒略歷和格里歷

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

小知識(shí)3：約化儒略日

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

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

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

#cal 9 1752

你會(huì)看到這樣一個(gè)奇怪的月歷輸出：

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

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

• --noryes

Copyright © 2022 朔方
Powered by .NET 6 on Kubernetes