农历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 + b*t1 + c*t2 + d*t3,单位是秒 其中t1 = (y-Y1)/(Y2-Y1)*10, t2 = t1*t1, t3 = t1*t1*t1 对于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 = t*t,t3 = t2*t,t4 = t3*t A、低精度算法 L0(t) = 48950621.66 + 6283319653.318*t 弧度 B、中精度算法 L1(t) = [ 48950621.66 + 6283319653.318*t + 53*t*t + 334116*cos( 4.67+628.307585*t) + 2061*cos( 2.678+628.3076*t)*t ] / 10000000 弧度 C、高精度算法 L2(t) = [ 48950621.66 + 6283319653.318*t + 52.9674*t2 + 0.00432*t3 - 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 弧度 最后两行分别为光行差和章动 四、太阳黄经速度 平速度: v0 = 628.3319653318 即时速度:v1 = 628.332 +21 * sin(1.527+628.307585*t) 速度的单位是“弧度/儒略世纪”即“弧度/36525天” 注意,平速度比即时速度的精度要高得多,务必保留足够的有效数字,否则将带来严重的计算误差。 五、节气时刻计算 以上天体黄经时间的函数,即L = f(t),所谓的求节气时刻就是已知L求t,显然这是在求解一个关于t的方程。伟大的英国天文学家物理学家牛顿给出了一种非常有效的迭代算法:牛顿求根法。用这种方法,求t所花费的时间仅是求f(t)花费时间的1.2——1.3倍。设某个节气对应的黄经为W,那么算法如下。 牛顿迭代算法设计: 第1步迭代:t = 0 第2步迭代:t = t + ( W – L0(t) ) / v0 第3步迭代:t = t + ( W – L1(t) ) / v1(t) 第4步迭代:t = t + ( W – L2(t) ) / v1(t) 误差:算法误差2分钟以内,实际找到的误差一般在30秒以内,平均15秒 注意:W指的是太阳黄经。1999年春分对应W=0,以后每W每增加15度对应下一个节气。迭代的的结果是力学时,单位是儒略世纪数。最后结果还应转换为北京时间,即:JD = J2000 + t*36525 - △T/86400 + 8/24 最后使用“儒略日数转公历”所述方法得到节气的日期和时间。 六、计算结果比较 为了进行误差比较,下文列出2007年的24节气,并与《寿星天文历》比对。《寿星天文历》是笔者制作的一款精度优于1秒的农历工具,已发布于互联网上,其算法与本文类似。
[参考文献] 1、Pierre Bretagnon与Gerard Francou 《VSOP87行星运动理论》 2、Jean Meeus《Astronomical.Algorithms》 3、NASA下属JPL实验室的《DE405/406星历表》 |
|