节气及农历计算程序 VSOP87+ELP/MPP02
版本:小鬼1.1(javascript), 设计:许剑伟, 2008.3, 电话:13375051698, 地址:莆田十中 本程序使用现代天文算法实现历算过程中最为困难的部分:日月位置的计算。当你得到任意时刻日月位置后,即可轻松推算万年日历。
公历年: 节气 手表时 中气 手表时 农历月 朔的手表时 大雪: 2010-12-07 13:38:22 冬至: 2010-12-22 07:38:29 12月大: 2011-01-04 17:02:24 小寒: 2011-01-06 00:54:41 大寒: 2011-01-20 18:18:38 正月大: 2011-02-03 10:30:33 立春: 2011-02-04 12:32:59 雨水: 2011-02-19 08:25:22 二月小: 2011-03-05 04:45:53 惊蛰: 2011-03-06 06:29:54 春分: 2011-03-21 07:20:40 三月大: 2011-04-03 22:32:10 清明: 2011-04-05 11:11:53 谷雨: 2011-04-20 18:17:24 四月大: 2011-05-03 14:50:30 立夏: 2011-05-06 04:23:11 小满: 2011-05-21 17:21:10 五月小: 2011-06-02 05:02:35 芒种: 2011-06-06 08:27:21 夏至: 2011-06-22 01:16:29 六月大: 2011-07-01 16:53:57 小暑: 2011-07-07 18:42:02 大暑: 2011-07-23 12:11:49 七月小: 2011-07-31 02:39:45 立秋: 2011-08-08 04:33:29 处暑: 2011-08-23 19:20:38 八月小: 2011-08-29 11:04:04 白露: 2011-09-08 07:34:12 秋分: 2011-09-23 17:04:36 九月大: 2011-09-27 19:08:49 寒露: 2011-10-08 23:19:01 霜降: 2011-10-24 02:30:17 十月小: 2011-10-27 03:55:57 立冬: 2011-11-08 02:34:52 小雪: 2011-11-23 00:07:49 11月大: 2011-11-25 14:09:44
说明: 1.算法及相关理论 ·太阳位置计算使用法国巴黎天文台的VSOP87行星理论,在前后数千年内可精确到分。 ·月球位置计算采用法国巴黎天文台的ELP/MPP02月球运动理论。 ·地球章动计算采用IAU1982的章动理论序列。 ·使用ELP2000-82B加VSOP87计算的出的农历及节气与刘国安先生所编写的"日梭万历"结果一置。 而本程序使用更精确的算法,远期的精度明显优于"日梭万年历" 2.计算时使用地心力学时(与原子时相当),它与协调世界时(即UTC时或手表时)存在一定的差值ΔT,ΔT随时间推移而发生距大变化。 ΔT的计算部分采用美国国家航空航天局网站上的算法。ΔT的远期估计从法国天文学家提供的skymap10.5星历程序中提取2150到 6000年之间的9个数据,并使用最小二乘法重新拟合出ΔT的多项式表达算式。 3.岁差计算采用IAU2000中关于岁差及章动中的B03表达式。程序中仅提供Date黄道坐标的及Date赤道坐标计算。并且这里黄道是地 月质心(EMB)瞬时平黄道,它是EMB公转轨道在天球上的扩展,黄经起算点是EMB的在天赤道上的升分点(地心坐标看太阳)或降分点 (日心坐标看地球),黄纬的起算是从平黄道开始,当地球过经度为零的黄经圈时,黄纬不为零,这是因受到月球扰动,地球在平黄道上 小幅上下波动,所以日过升降分点时,太阳黄纬不是零 4.程序中已对地球运动引起的光行差做了修正。 5.精度:ELP/MPP02理论结合DE405得到的远期6000千年范围内的精度为几个角秒(相当于数秒钟),ELP/MPP02理论结合LLR(激光测距系 统)在近期精度高达毫秒级(已经考虑了相对论修正等), 然而程序是javascript脚本,运行速度非常慢,所以不得不做截断处理,降低 一些精度以图提高计算速度。最后得到月球位置与原理论比较,在公元-1000至5000范围内月球计算平均误差为4角秒,最大13角秒。 行星计算也做了截断处理,最后得到太阳位置与原理论的千年误差:经度0.3角秒(最大0.8),纬度0.06角秒(最大0.2),距离万分之0.1 (最大0.3)。太阳每分移动2.5角秒左右,所以程序的精度仍非常高。当然,如果使用UTC时间表达数千年后的精度,误差是数分至数小 时。自转周期只要发生了毫秒级的变化,积累几千后就会变为几十个小时。没有人知道几千年后地球自转的速度,也许3000年后地 球自转周期是86400秒,也可能是86401秒,反正它时快时慢。所以用手表时表达远期的历算精度,误差巨大。 6.这里用到的大部分数据是从笔者用C++编写了高精度的星历程序中导出的 程序中的农历说明: 农历计算要点: ·规定1:日月合朔时刻所在的UTC日期为新月份初一,而不是上月的最后一天。 所以,在合朔那一天正好是节气或中气,那么该气属本月(的初一),不是上月的 例:冬至200x-12-22 10:00:00, K月的起止合朔时刻200x-11-23 xx:xx:xx和200x-12-22 11:00:00 K月的起止时刻范围已包含了冬至,但从日期看,冬至却不属于K月,冬至是下一个月的 ·规定2:农历年首始于冬至所在的朔望月。 古人通过冬至定年首,让月亮历也能表达太阳历的信息 可见冬至永远属于新年,新年首月份必含冬至这个中气,含有中气的月分必然不是闰月 根据此规定易知,上述的K月也不属于新年 ·规定3:一年内如果有13个朔望月,则该年是农历闰年 第一个月含有冬至的月份,第13个月不能含有冬至,那个冬至是明年的 ·规定4:闰年中首个无中气的月份定为闰月,闰月名称同上个月,闰月以后的几个月,月序号全部减1 ·算法描述: 1.先用行星运动理论计算出中气,首个中气时该的太阳地心视黄经为270度,第2个中气300度,第3个330... 所谓地心视黄经计算可按如下(注意有多种坐标变换方法,以下方法仅供参考): (1)先计算日心Date平黄道中地球位置坐标 (2)转到地心坐标,得到地心Date平黄道中太阳的位置,最好变换为球面坐标(如果未补岁差,此时补) (3)补上章动得到地心Date黄道的太阳位置坐标 (4)补光行差得到地心Date黄道视位置,视位置中的黄经就是这里说是地心视黄经 2.利用月球运动理论计算月球坐标,并想办法转换为地心视坐标,地月光行时仅1.3秒,无需处理光行差 3.冬至所在的UTC日期保存在A[0],根据"规定1"得知在A[0]之前(含A[0])的那个UTC朔日定为年首日期 冬至之后的中气分保存在A[1],A[2],A[3]...A[13],其中A[12]又回到了冬至,共计算14次中气 4.连续计算冬至后13个朔日,即起算时间时A[0]+1 14个朔日编号为0,1...12,保存在C[0],C[1]...C[12] 这14个朔日表示编号为0月,1月,...12月0月的各月终止日期,但要注意实际终止日是新月初一,不属本月 这14个朔日同样表示编号为1月,2月...的开始日期 设某月编号为n,那么开始日期为C[n-1],结束日期为C[n],如果每月都含中气,该月所含的中气为A[n] 注:为了全总计算出13个月的大小月情况,须算出14个朔日。 5.闰年判断:含有13个月的年份是闰年 当第13月(月编号12月)终止日期大于冬至日, 即C[12]〉A[12], 那么该月是新年,本年没月12月,本年共12个月 当第13月(月编号12月)终止日期小等于冬至日,即C[12]≤A[12],那么该月是本年的有效月份,本年共13个月 6.闰年中处理闰月: 13个月中至少1个月份无中气,首个无中气的月置闰,在n=1...12月中找到闰月,即C[n]≤A[n] 从农历年首的定义知道,0月一定含有中气冬至,所以不可能是闰月。 首月有时很贪心,除冬至外还可能再吃掉本年或前年的另一个中气 定出闰月后,该月及以后的月编号减1 7.以上所述的月编号不是日常生活中说的"正月","二月"等月名称: 如果"建子",0月为首月,如果"建寅",2月的月名"正月",3月是"二月",其余类推 程序中公历的说明: (1)公历的历月为十分混乱,1月31天,2月28或29天,3月31天.... (2)公历最历害的是:实现了400年97闰,平均年长为365.2425,已非常接近现代的365.2422,而且置闰规则简单。 (3)公历分为两个阶段1582年10月31日之前为儒略历,实行4年1闰制,1582年11月10日后实行400年97闰,是格里高利历。两阶段之间少了 10天,是因为以前的4年1闰不精确,闰过头了,造成日历时比真太阳慢了10天,到了1582年10月31日进行10天负闰。历史上,闰过头或闰 过少的的事件很多,各国也有所不同。这也说明,人们对太阳系运动规律的认识是不断提高的,而不是一步到位的。 上面讲至正负闰,正闰指给该年份加一些天数,负闰指给该年份减一些天数。 (4)公历是太阳历,农历是太阳历与月亮历的合成,农历中的太阳历非常严格的表达天文学上的太阳位置,而公历不严格的。公历比较简 单,而且太反映了与气候直接相关的太阳运动,所以公历在世界范围内广范使用。
|