寿星天文历是一款采用现代天文算法制作的天文、历算程序,可以方便地进行公历、农历、回历三历之间的转换;提供了历谱数据导出功能;提供精确的日月食等天象的计算功能;含有行星、恒星星历表计算。它提供公元-4712年到公元9999年的公历日期查询功能以及-3000年至+3000年的农历查询功能,其中-721年到1960年的农历数据已经与张培瑜的《三千五百年历日天象》、陈垣的《二十史朔闰表》、方诗铭的《中国史历日和中西历日对照表》核对;它还含有公元前2000多年以前到今的基本年号;含有二千多个国内城市的经纬度,并且用户可根据自已的需要扩展经纬度数据。
本软件是一款精准的年代跨度大的日历工具,可作为一般的实用日历工具,对史学家、考古学家、历算工作者、天文爱好者均有较大的参考价值。
为了方便爱好者查阅数据,寿星提供了常数系统、日食概略计算等。
[历谱说明]:1645年农历八月开始使用《时宪历》,七月及之前为《授时历》
[日食月食]:天文历中提供的图表,可以精确的确定日月食情况。计算结果与《2008年中国天文年历》比对,日食误差在1秒左右,月食误差在3秒以内,但应注意到,《2008年中国天文年历》提供的月食精度只有0.1分,即+-3秒。事实上,受到地球大气的影响,月食的精度无法计算得很高,所以计算时把地球看成了圆球,而不是椭球。《寿星天文历》日月食预报及计算的本质就是把日月坐标以图表的形式展示出来,并给出日食中心线、地方食表、日食南北界及月球本影半影的计算。通过日食图表与用户互动观察,有利于学习一些比较专业的科普知识,体验一下DIY的乐趣。另外:如果您利用寿星天文历得到的日食计算结果与媒体上公布的或《中国天文年历》上的数据相差2秒钟以上,多半是地理坐标数据的采用值不同或是媒体的错误造成的,请认真核对。比如《2008年中国天文年历》第526页,给出了2007年3月19日阿勒泰市(新疆)的日食计算及结果,其实这里有个低级错误,日期应改为2008年8月1日。
[地标数据]:城市经纬度数据库做了压缩处理。要扩展经纬数据库,可使用附件中的【城市经纬压缩器】生成代码并替换程序中的JWv数组。程序中现有的经纬度数据本来只有几百个,后来网友郑彬给我一份3000多个城市的数据,比较全面,因此本版使用他的数据。
[时区数据]:各地时区数据库做了小量压缩处理。如要扩展,可使用附件中的【时区合成器】生成代码并替换程序中的SQv数组。时区数据库是“中华农历网”netghost先生提供。
[年号纪年]:天文历中的纪年数据参考万国鼎《中国历史纪年表》,有些朝代存在多组帝王,程序中不一定全部纳入。年号纪年问题比较罗嗦,本人史学功底肤浅,无法严格的解决此问题,仅在程序中加入南北朝(约公元420年)至今的年号数据,未能保证其完全正确。在V2.04版本以后,补上了由“中华农历论坛”ldlcau(亦是“国学数典”论坛的ldlcau)先生提供的公元420年以前的年号数据。从V2.04版开始,程序中提供了年号数据管理工具,详见附件中的【纪年管理器】。
[日数标尺]:本天文历中使用两种日数标尺。其一是公历,其二是儒略日数。为了简化历算过程,本天文历的历算代码全部使用“儒略日数”作为日数标尺。日历数据的显示则使用公历,它也叫格里历。严格的说,格里历始于1582年10月15日,1582年10月4日及以前的公历是4年1闰的儒略历。公历是在最近一百两年受到西方政治、经济优势的影响才逐渐广泛使用,显然本天文历把“公历”外延到千年以前只是为了建立一个相对统一的连惯的日数标尺,而不是说很久以前全世界就在使用这种历法。再者,“公历”中的儒略历也不完全是古罗马颁行“儒略历”,因为古罗马儒略历曾有过置闰失误。
[开源软件]:寿星天文历是开放源代码的,我并不打算作任何加密处理。月历生成部分涉及多种历法,而且还要处理纪年、节日信息、古历算等,所以相应的程序比较长。天文算法部分的程序比较短,只要算法已知,编写起来是比较容易的,困难的是如何设计出我们所需高精度、高速的算法。要设计算良好的天文算法,需参考较多资料并作认真的推导、验证等,工作量比较大,所以最好不要更改这部分代码,以免发生错误。
二、古代农历与寿星天文历
平气平朔日推算:只要知道当时采用的平朔望月(或相邻平气)长度k及历算起始参数b,就可推算出古历的气朔数据,当然,不同的历法所采用的k和b是不同的。例如:假设公元100年1月1日为朔日,相应的平朔时刻定为这天的12点(100年1月1.5日),平朔的k=29.53。那么下一个平朔时刻=100年1月1.5日+29.53=100年2月1.03日,则相应的朔日发生在2月1日。
显然,推算平气平朔的代数表达是:D=[k*n+b],式中D是气朔的日期(日数),k为气朔长度,n为气朔积累计数值,b为初始值,中括号表示取整。为了得到各个颁行历所使用的k和b,寿星天文历选择了一种快捷方法,即在现有的历谱数据中提取各历法的k、b采用值。在数十万的历谱数据中找k和b,简值是大海捞针。因此本人借助计算机强大的数据处理能力,反向推算古历算诸参数。本程序附带的反向推算的算法是通用的,不仅适用已知的连惯历谱,也适用于不连续的片段历谱。
反向推算的结果可能与古代历算的实际采用值存在少量差异。比如,对“四分历”历谱(公元85年以后1百多年的历谱数据)反推得到平朔望月长度值是29.53085097。历书记载的值是29+499/940=29.53085106。但本天文历使用的是前者,而不是后者。因为历算时全部采用浮点数,存在计算机固有的截断误差,直接使用499/940是危险的。如,对于分数499/940乘940得正确值499,而计算机得到的可能是489.9999999999,取整数日以后变为489,误差1天。相反,k值取前者是安全的,因为k和b的值在反推过程中已在历法适用范围内得到完全验证。关于反向推算的算法详见主目录下的《离散序列的直线拟合》章节。
反向推算过程1:(-103年至1960年)
(1)先从徐先生(网名“春光”)发给我的《古代历法摘要(上)》获得“中国历法表”,在这个表中得知中国曾经版发的历法名称及适用范围。春光的一句话给我很大的帮助:“主要参考律历志”。
(2)在网上下载《三千五百年历日天象》pdf书,把古代日历部分传换为bmp格式。通用ocr软件无法准确的将bmp日历转换为文本格式日历数据表,所以利用C++Builder编写一个图片转文字的专用程序,将这本书扫描成文本格式的日历表。在扫描过程中,遇到识别困难的,在识别后的文字前加一个问号。这项工作整整历时5天。误码率为1%,误码的部分通常会有个问号,所以真正误码的不到千分一。当初转化时,只从-100年开始,所以其中的-103年到-101的这3年的日历采用手工输入。
(3)编写javascript程序,对扫描后的数据进行校验,以消除误码并找出原书的错误。校验方法:a.大小月与朔日的日期进行校验;b.1645年以后的一些关键数据还与寿星天文历早期版本比对;c.节气的日期与节气干支比对;d.节气干支与朔日干支校验;e.相邻节气长度粗大误差校验。所有校验不通过的给出报告并进行人工校对,直到不发生错误报告为止。校验过程中,发现原书也有不少错误。
(4)最强有力的校验:分段线性拟合校检!1645年以前的日历采用平气法,也就是说在某一个历法实行期间,相邻节气交接时刻之差(或朔时刻之差)为常数,各交接时刻可以连成直线。由于原书没有给出交接时刻,所以采用了一些变换方法。
(5)古代定气的算法:本天文历使用开普勒椭圆轨道定气,同时考虑了岁差、章动、轨道要素的长期变化。但必竞使用了现代数据,与古代历算参数有点差异,所以程序中对定气朔进行了订正,使之与古历相符。
也许有人会问,为什么不考虑使用古代原来的算法?实际上,为了提高历算精度,清代的《历象考成》已经使用了第谷的宇宙体系,在乾隆七年(公元1742年)以后的《历象考成后编》,开始应用了开普勒行星运动第一、第二定律(但椭圆焦点上的是地球而不是太阳)。因此大约在1734年以后,古代定气与现代计算结果的已经相差无几。这就造成简单的加减乘除根本不能解决问题,所以宁愿采用“现代算法+订正”思路来解决问题。
反向推算过程2:(-721年至-104年)
(2)对气朔表做线性拟合。对闰年表做线性拟合。拟合的过程中同时校验查错。
(3)拟合以后发现,该书中的春秋、战国分别使用一个平朔线性算式。实际上,原数据个别数据错误造成一个线性算式无法拟合通过,我已对其改正。从原理上分析,上元取值不同造成的的“错误”不会是个别的,所以《三千五百年历日天象》中的这些错误多半是原作者算错或笔误了或颁行历法的错误。不过,颁行历法发生错误的可能性要小很多。
反向推算的结果:
在没用足够的古代历书的帮助下,直接使用历谱数据,通过以上方法就可重现古代不同历法平气平朔的历算参数及定气朔的修正表。所有参数确定后,还特别校对了一下1645年以后定气定朔颁行历与现代算法计算结果的差异。经过以上拟合校验,最后得到古代历谱计算所需的数据只有4k,并使得计算结果的准确性超过了《三千五百年历日天象》本身,可利用它辅助校验历书数据。
上述平气、平朔的计算公式表达为:D=k*n+b,式中n=0,1,2,3,...,N-1;下面给出各朝代颁发历法的k与b的参数,每行第1个参数为b,第2参数为k,结果D表达为儒略日数;h表示k不变b允许的误差,如果b不变则k许可误差为h/N;N是历法颁行期间朔或气的总个数。
举例说明:求公元88年2月15(儒略日数为1753245)附近的朔日的具体公历日期。
88年的颁行历为后汉书·律历志的四分历,它的平朔参数b=1752148.041079,k=29.53085097。先把朔日D估计为88年2月15日。又因D=k*n+b,所以n=(D-b)/k=37.15,四舍五入取整得n=37。易得准确日期为D=k*n+b=1753240.68,转为公历得公元88年2月11号4时,朔日在11日。
1501至1940颁行历的朔日与今算得天象不相符的数据:
以下古历部分数据来自《二十史朔闰表》
《二十史朔闰表》中有错误的数据:[书中错误]:1808-12-17原书为1808-12-11误差多达6天,可能是作者笔误。1770-06-23原书为1770-06-22造成月小于29天。1505-12-25原书为1505-12-23造成月小于29天。1579-05-25原书为05-23造成月小于29天。[书中如下数据可能有误]:1821-08-2723:18:45原著为1821-08-28,当时的测算水平已经很高,误差应小于半小时(一般小于10分钟)。1588-03-2710:46:16原著为1588-03-26,竞然相差半天左右,当时出现如此误差十分费解,原书可能有误。《三千五百年历日天象》中有错误的数据:
带*表示改后数据与方诗铭《中国史历日和中西历日对照表》相同,带#号表示改后与方诗铭的不同。
从寿星万年历2.0版本(后期改名为寿星天文历)开始,处理古代农历置闰问题时,取消以往1.x版本的置闰修正算法,而是在准确算出古代气朔表的基础上统一使用“无中置闰法”(-103年至今)或“19年7闰法”(-721年至-104年)得到。
十九年七闰法:7个闰年均匀安插在19个年整数中,闰年的末月置为闰月。
平气无中气日置闰法:不含平中气日的月份置为闰月。这样,平气平朔历法中,7闰月均匀安插在235个月中。定气无中气日置闰法:使用定气法之后,把包含冬至日的月份定为一年中的首月。定年首的后果是有的年份有13个月,有的只有12个月。在含有13个月的年份中,采用无中气日置闰法,并且只对首个无中气的月份置闰。
“定气无中置闰”兼容“平气无中置闰”,反之不兼容。也就是说,平气无中置闰也可加入“含有13个月”这个条件,不会影响历算结果。因为,在平气历法中不可能出现一个月包含两个中气的现象,被置闰的年份中必含有13个月。这使算法变得比较统一,降底了程序的复杂性。本天文历严格执行上述置闰法,得到与实闰相同结果。
月建问题:
史记·天官书:“值斗杓所指,以建时节”。在古代,常以北斗七星的斗柄方向判断四季,历法与天象存在直接的联系,这里的“建时节”可用“建何月”进行历书表达。换句话说,斗柄指向某一角度,当前月定为“建某月”。往后的历法以标准的太阳历定四季(以12中气为准,斗柄只是参考),即太阳黄经达到某一位置角度(中气定义的角度),当前月定为“建某月”。
古历法中,农历1个月中含冬至的是建子月,含大寒的建丑月,其余类推。没有中气的月分,就没有独立月建。因此,在今人进行农历的历月计算时,月建实际上就是一个冬至起算的中气序数(不含中气的农历月没有独立的序数,其序数与前一个月相同)。
月建与农历月名称:
我们不妨把“正月、二月、三月等”看作月建的别名。并且建子为十一,建丑为十二,建寅为正……由于某些原因,有时古代帝王随意更改了别名。本软件也相应作了更改:-721年12月17日至-479年11月春秋战国历采用19年7闰,闰年的末月置闰并取名闰十三。年首为正月。-221年10月31日至-104年11月秦汉历同样采用19年7闰,闰年的末月置闰并取名后九月。年首为十月。以上两行所述的年代,其月建别名计算方法:月建序数=(月积累数-置闰积累数+C)mod12,式中C为某一常数。008年01月15日至023年12月02日:建子为十二,建丑为正月,建寅为二月,其它顺推。237年04月12日至239年12月13日:建子为十二,建丑为正月,建寅为二月,其它顺推。689年12月18日至701年01月14日:建子为正,建寅为一月,其它不变。761年12月02日至762年03月30日:建子为正,其它顺推。
由于上述月建别名的更改,造成排算的结果是239年12.13为十二月,240年01月12日变回原来的月建别名,这样240年01月12日建丑十二月,显然连续的出现两个十二月,本软件中把上个月表示为“拾贰月”。同样的问题也出现在23年底。
237年正月开始使用新的历法,造成前一月只有28天。
三、《寿星天文历》开发过程
我在网络上找了很久,参阅的天文资料(全部在网络上找的):行星运动VSOP87方法(外文)、月球运动ELP2000-82B(外文)、月球运动ELP/MPP02方法(外文)、行星运动IPS2000方法(外文)、DE系列星历表(外文)、瑞士星历表(外文)、《天文算法》(除了第41、42、43、44章,本书已全部译完,有需要的我可以发一个)、北师大高健的《球面天文学讲议》ppt、《球面三角基本公式》pdf、《有限差分法》doc、《数值方法》、《天体力学引论》(前几年在超星下载打印的)、《IAU1982岁差章动》《IAU2000岁差章动》(外文pdf)、陈垣《二十史朔闰表》pdf、万国鼎《中国历史纪年表》pdf、唐汉良和余宗宽《日月食及其计算概要》pdf、《fortran教程》(现代天文算法一般是用fortran写的)等。要想利用天文算法实现农历计算,这些书籍、资料或多或少是要读的。为了方便验算比对,在农历网又托人购买了《2008年中国天文年历》。
当我分析了月球的ELP/MPP02的fortran程序之后,就用c++写出相同的程序,并利用这个程序导出javascript所需月球数据及算法,所以本程序所用的算法与原版的ELP/MPP02有所不同:是对数据序列进行转换处理,使得序列的表达形式与VSOP87的序列的形式相同,大大降底了算法的难度,速度也有所提升,缺点是数据增加了20%——30%
四、《寿星天文历》的精度、可靠性及性能
历算精度是一个复杂的问题,以下从多个方面分析历算的精度问题。
协调世界时(UTC):秒长使用原子时,所以它是力学时。可是,为了与世界时同步,UTC在每年的年底做了一个小动作,即闰秒。所以,从长远看UTC属世界时范畴,在一年内看UTC属力学时范畴。
2008年以后ΔT的加速度估计值取31(采用瑞士星历表的),如果要使用skymap11的,请把程序中的jsd改为29。
可见,对于遥远的过去与遥远的将来,历算的误差主要来自ΔT。对于过去的几百年,误差主要来自于日月坐标。对于日历的准确性,主要取决于核对时否认真以及历书著作的可靠性等。
程序的大小:核心程序大约30K,其它是注释、说明、城市经纬数据、纪年数据、页面等。也就是说,如果做成精简版,只需30几K,后期版本增加了九行星和日月食、地图数据,程序变大了很多,约500k。
本程序是开源的,你可以使用其中的任意部分代码,但不得随意修改“天文算法(eph.js)”及“农历算法(lunar.js)中古历部分的数据及算法”。一旦修改可能影响万年历的准确性,如果你对天文学不太了解而仅凭对历法的热情,请不要对此做任何修改,以免弄巧成拙。
特别感谢中华农历网浪淘沙等人帮助我调试软件,作为天文历法,有幸结识这么多同好,一生难忘。
作者:许剑伟,2008年11月于家里。xunmeng04@163.com,13850262218