|
节气及农历计算程序 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)公历是太阳历,农历是太阳历与月亮历的合成,农历中的太阳历非常严格的表达天文学上的太阳位置,而公历不严格的。公历比较简
单,而且太反映了与气候直接相关的太阳运动,所以公历在世界范围内广范使用。
|
加好友 
发短信
Post By：2011/7/19 22:53:00
Post By：2011/7/22 6:27:00
