农历24节气算法

2016-12-06 • panger • 1条评论 • 2,170 次浏览

农历24节气算法

[摘要] 古老而又现代的中国农历，是一种天文学性质的阴阳历。本文采用VSOP87星历算法并结合运动学方法以及牛顿求根法得到24节气的精确时间，揭开农历计算的神秘面纱。

[关键字] 农历算法、星历、节气

[正文] 计算中国农历，首先要计算出二十四节气时刻。在计算机问世之前，二十四节气的许算是非常复杂的。随着计算机及互联网的普及，美国航空航天局、法国巴黎天文台各自在网络上发布了精密星历表的计算方法，这使得民间计算农历成为可能。本文以法国巴黎天文台的VSOP87算法为基础，给出中国农历的二十四节气算法。

　　在农历中，太阳黄经为0度时，对应春分节气。相邻节气对应的太阳黄经相差15度。一周年内，太阳黄经从0度变化到360度，共有24个节气。

一、时间标尺——儒略日数计算

　　计算星历之前首先要解决时间尺问题。公历规定平年365日，闰年366日。1582年10月4日以前，公历规定每4年设置一个闰年，平均年长度365.25天，这期间的公历称为儒略历。在1582年10月15日之后实行格里高利历，规定每400年97闰，平均年长度为365.2425天。
由于儒略历存在严重的“多闰”问题，到了1582年，公历跑快了10天左右，当时就人为调整了10天，并从此实行格里历。因此务必注意1582年10月4日（儒略历）的下一日为1582年10月15日（格里历）。就是说1582年10月份少了10天。

　　在儒略历中，能被4整除的年份为闰年，这一年有366天，其它年份为平年（365天）。如900年和1236年为闰年，而750年和1429年为平年。

　　格里高利历法也采用这一规则，但下列年份除外：不能被100整除的年份为平年，如1700年，1800年，1900年和2100年。其余能被400整除的年份则为闰年，如1600年，2000年和2400年。
儒略日数（简称儒略日）：

　　儒略日数是指从公元 -4712 年开始连续计算日数得出的天数及不满一日的小数，通常记为 JD (** )。传统上儒略日的计数是从格林尼治平午，即世界时12点开始的。若以力学时（或历书时）为标尺，这种计数通常表达为“儒略历书日”，即JDE (**)，其中E只是一种表征，即按每天86400个标准秒长严格地计日。例如：

1977年4月26.4日 UT = JD 2443259.9

1977年4月26.4日 TD = JDE 2443259.9

儒略日的计算：

设Y为给定年份，M为月份，D为该月日期（可以带小数）。

若M > 2，Y和M不变，若 M =1或2，以Y–1代Y，以M+12代M，换句话说，如果日期在1月或2月，则被看作是在前一年的13月或14月。

对格里高利历有：A = INT（Y/100） B = 2 – A + INT(A/4)

对儒略历，取 B = 0

儒略日即为：

JD = INT(365.25(Y+4716))+INT(30.6001(M+1))+D+B-1524.5

使用数值30.6取代30.6001才是正确的，但我们仍使用30.6001，以确保总能取得恰当的整数。事实上可用30.601甚至30.61来取代30.6001。例如，5乘30.6精确等于153，然而大多数计算机不能精确表示出30.6，这导致得出一个152.999 9998的结果，它的整数部分为152，如此算出的JD就不正确了。

由儒略日推算历日：

将JD加上0.5，令 Z 为其整数部分，F 为尾数（小数）部分。

若 Z < 2299161，取A = Z

若 Z 大于等于2299 161，计算

α=INT((Z-1867216.25)/36524.25)

A=Z+1+α-INT(α/4)

然后计算

B = A+1524

C = INT((B-122.1)/365.25)

D = INT(365.25C)

E = INT((B-D)/30.6001)

该月日期（带小数部分）则为：

d = B – D – INT(30.6001E) + F

月份m为：

IF E < 14 THEN m = E – 1

IF E=14 or E=15 THEN m = E – 13

年份为y：

IF m >2 THEN y = C – 4716

IF m =1 or m=2 THEN y = C – 4715

这个公式里求E时用的数30.6001不能代之以30.6，哪怕计算机没有先前所说的问题。否则，你得到的结果会是2月0日而不是1月31日，或者4月0日而不是3月31日。

值得记住的一个常数是：2000年1月1日12:00:00的儒略日数是J2000 = 2451545

二、力学时与世界时的差值（deltat T）计算

一般的，可以把手表时（UTC）近似看作世界时（UT），二者的主要差别在于时区。如北京手表时8点对应世界时0点。世界时与地球自转严格同步，但有趣的是，我们的手表时实际上称为协调世界时，它的秒长是原子钟的秒长，由于地球自转速度不均匀，时快时慢，这就注定手表时与地球自转不完全同步。现在，地球自转速度正在变慢，我们不得不在某些年份的年末把手表拨慢1秒，使得手表时更好的与地球自转同步，并美言为“跳秒”。力学时是根据太阳系的动力学原理导出的，是一种均匀的时间系统，其秒长与原子钟的秒长相同。因此，协调世界时（UTC）与世界时（记为UT）其本同步，但力学时（记作TD）与世界时不太同步，二者的差值记作deltat T或记作△T。利用直接的天文观测可以得知每年的△T，利用古代的日月食观测资料可以反推古代的△T。所有年份的△T计算出来后，可以拟合出以下多项式表达，使得△T的计算更快捷，计算结果的单位是秒。

我们利用下表可以严格计算△T（即△T =TD – UT）

年份 a b c d

-4000,108371.7,-13036.80,392.000, 0.0000

-500, 17201.0, -627.82, 16.170,-0.3413

-150, 12200.6, -346.41, 5.403,-0.1593

150, 9113.8, -328.13, -1.647, 0.0377

500, 5707.5, -391.41, 0.915, 0.3145

900, 2203.4, -283.45, 13.034,-0.1778

1300, 490.1, -57.35, 2.085,-0.0072

1600, 120.0, -9.81, -1.532, 0.1403

1700, 10.2, -0.91, 0.510,-0.0370

1800, 13.4, -0.72, 0.202,-0.0193

1830, 7.8, -1.81, 0.416,-0.0247

1860, 8.3, -0.13, -0.406, 0.0292

1880, -5.4, 0.32, -0.183, 0.0173

1900, -2.3, 2.06, 0.169,-0.0135

1920, 21.2, 1.69, -0.304, 0.0167

1940, 24.2, 1.22, -0.064, 0.0031

1960, 33.2, 0.51, 0.231,-0.0109

1980, 51.0, 1.29, -0.026, 0.0032

2000, 63.87, 0.1, 0, 0,

2005

表中每一行适用一定的年代范围，如第1行适用于公元-4000年到-500年，第2行适用于公元-500到-1500年，其它类推。每行的起始年份记作Y1，终止年份记作Y2，如果年份y在Y1到Y2之间，那么该年的deltat T表达为：

△T = a + bt1 + ct2 + d*t3，单位是秒

其中t1 = （y-Y1）/（Y2-Y1）10， t2 = t1t1， t3 = t1t1t1

对于2005年以后的deltat T是未知的，要做外推计算：

2005至2014年建议使用1995到2005年期间△T的平均增速计算，即：

△T = F(y) = 64.7 + (y-2005) * b，其中速度 b = 0.4

2114年以后可以使用二次曲线外推

△T = f(y) = -20+ a * [(y-1820)/100]^2 ，其中加速度a = 31

2114年到2014年之间的外推，可以在上面两个外推算式的基础上做一次的曲线连接，使之连续即可。比如可
以这么计算：

△T = f(y) + (y-2114) * [f(2014) – F(2014)] /100

以下数值可供程序验证参考

2008年△T = 66.0秒

1950年△T = 29秒

500年 △T = 5710秒

三、太阳视黄经（真分点视坐标）

算法基于VSOP87半解析法。

力学时t为J2000起算的儒略世纪数，t2 = tt，t3 = t2t，t4 = t3*t

A、低精度算法

L0(t) = 48950621.66 + 6283319653.318*t 弧度
B、中精度算法

L1(t) = [ 48950621.66 + 6283319653.318t + 53t*t

334116cos( 4.67+628.307585t)
2061cos( 2.678+628.3076t)*t ] / 10000000 弧度

C、高精度算法

L2(t) = [ 48950621.66 + 6283319653.318*t

52.9674t2 + 0.00432t3 – 0.001124*t4
+334166 * cos( 4.669257+ 628.307585*t)

+3489 * cos( 4.6261 + 1256.61517*t )
- 350 * cos( 2.744 + 575.3385*t)
- 342 * cos( 2.829 + 0.3523*t)
- 314 * cos( 3.628 + 7771.3771*t)
- 268 * cos( 4.418 + 786.0419*t)
- 234 * cos( 6.135 + 393.021*t )
- 132 * cos( 0.742 + 1150.677*t )
- 127 * cos( 2.037 + 52.9691*t)
- 120 * cos( 1.11 + 157.7344*t)
- 99 * cos( 5.23 + 588.493*t )
- 90 * cos( 2.05 + 2.63*t )
- 86 * cos( 3.51 + 39.815*t )
- 78 * cos( 1.18 + 522.369*t )
- 75 * cos( 2.53 + 550.755*t )
- 51 * cos( 4.58 + 1884.923*t )
- 49 * cos( 4.21 + 77.552*t )
- 36 * cos( 2.92 + 0.07*t )
- 32 * cos( 5.85 + 1179.063*t )
- 28 * cos( 1.9 + 79.63*t )
- 27 * cos( 0.31 + 1097.71*t )
+2060.6 * cos( 2.67823 + 628.307585*t ) * t

+43.0 * cos( 2.635 + 1256.6152*t ) * t

+8.72 * cos( 1.072 + 628.3076*t ) * t2

-994 – 834 * sin(2.1824-33.75705*t)
64 * sin(3.5069+1256.66393*t) ] / 10000000 弧度

最后两行分别为光行差和章动