第一章 关于huangdaojiri.org天文历(万年历)

一、huangdaojiri.org天文历概述

huangdaojiri.org天文历是一款采用现代天文算法制作的天文、历算程序,可以方便地进行公历、农历、回历三历之间的转换;提供了历谱数据导出功能;提供精确的日月食等天象的计算功能;含有行星、恒星星历表计算。它提供公元-4712年到公元9999年的公历日期查询功能以及-3000年至+3000年的农历查询功能,其中-721年到1960年的农历数据已经与张培瑜的《三千五百年历日天象》、陈垣的《二十史朔闰表》、方诗铭的《中国史历日和中西历日对照表》核对;它还含有公元前2000多年以前到今的基本年号;含有二千多个国内城市的经纬度,并且用户可根据自已的需要扩展经纬度数据。

本软件是一款精准的年代跨度大的日历工具,可作为一般的实用日历工具,对史学家、考古学家、历算工作者、天文爱好者均有较大的参考价值。

如果以力学时作为时间标尺,与DE406星历表比对,太阳坐标的最大可能误差为0.1角秒,月亮坐标的最大可能误差为0.8角秒,平均误差为最大可能误差的1/6,这使得精确的日月食计算成为可能,

是目前世界上唯一一款整合了精确日月食过程计算的万年历。

在最近的几百年范围内,本软件的节气时刻及日月合朔时刻的平均误差小于1秒。

huangdaojiri.org还提供日出日没、月出月没时间查询等功能,与国内外著名的天文软件或天文年历比较,结果完全一致。此外,天文历中还即时显示了“月亮被照亮部分比例”、太阳和月亮的地心黄道坐标、地心赤道坐标、站心地平坐标等。其中日月的地心黄道、赤道坐标采用高精度算法,“月亮被照亮部分比例”的计算采用低精度的算法,但可以满足一般需求;地平坐标的计算已适当考虑了大气折射,由于大气的气压及气温变化具有许多不确定的因素,所以只考虑平均情况下的大气折射修正,在本软件中,地平坐标被描述为方位角与高度角,它是以观测点为中心的坐标,所以已经考虑了周日视差修正。

屏幕显示的方位角从正南向西测量

方位角与高度角已转换到是站心地平坐标,并且做了视差修正,同时在地平真纬度大于0时进行大气折射修正。

为了方便爱好者查阅数据,huangdaojiri.org提供了常数系统、日食概略计算等。

[历谱说明]:

1645年农历八月开始使用《时宪历》,七月及之前为《授时历》

[天文算法]:

huangdaojiri.org天文历在进行月历推算时采用了先进的现代天文算法,它是基于2002年巴黎天文台西尔特处Jean CHAPRONT 和 Gerard FRANCOU发表的ELP/MPP02月球运动理论,在几百年范围内,该理论与当今世界上公认的最精密的LE405/406系列星历表仅相差3毫角秒,在几千年范围相差不超过4角秒。由于农历的计算不需要这么高的精度,我们做了一些截断处理及实质性的简化处理,但仍可以满足高精度历算的要求。程序中采用了1988年法国巴黎天文台Pierre Bretagnon与Gerard Francou发表的VSOP87行星运动理论计算太阳坐标。由于该理论与DE405/406星历表存在少量差异,我们对VSOP87进行少量修正,提高了计算精度。节气及合朔时刻问题在天文学层面看来,就是已知天体的坐标反求时间的问题。在VSOP87或ELP中,天体的坐标是时间的函数z = f(t),所谓的求节气时刻或月相时刻就是已知z求t,显然这是在求解一个关于t的方程。伟大的英国天文学家物理学家牛顿给出了一种非常有效的迭代算法:牛顿求根法。用这种方法,求t所花费的时间仅是求f(t)花费时间的1.2——1.3倍,并且迭代引起的额外时间误差只零点零几秒。牛顿求根法在一般的《数值方法》书籍中均有介绍。

[日食月食]:

天文历中提供的图表,可以精确的确定日月食情况。计算结果与《2008年中国天文年历》比对,日食误差在1秒左右,月食误差在3秒以内,但应注意到,《2008年中国天文年历》提供的月食精度只有0.1分,即+-3秒。事实上,受到地球大气的影响,月食的精度无法计算得很高,所以计算时把地球看成了圆球,而不是椭球。《huangdaojiri.org天文历》日月食预报及计算的本质就是把日月坐标以图表的形式展示出来,并给出日食中心线、地方食表、日食南北界及月球本影半影的计算。通过日食图表与用户互动观察,有利于学习一些比较专业的科普知识,体验一下DIY的乐趣。

另外:

如果您利用huangdaojiri.org天文历得到的日食计算结果与媒体上公布的或《中国天文年历》上的数据相差2秒钟以上,多半是地理坐标数据的采用值不同或是媒体的错误造成的,请认真核对。比如《2008年中国天文年历》第526页,给出了2007年3月19日阿勒泰市(新疆)的日食计算及结果,其实这里有个低级错误,日期应改为2008年8月1日。

[日月出没]:

受到各地的地形地貌、大气状态、海拔高度等的影响,日出日没、月出月没时间无法计算得很准确,通常只能精确到1分钟,因此,软件中时间显示精确到到了秒数量级,但这并不表明日出日没时间的计算达到了这个精度。日月中天时刻受大气等因素的影响要小得多,所以可以精确到秒。日月的“出、中天、没”的计算不需要高精度日月位置坐标(因为一天中,日月在视野中的位置主要由地球自转决定),程序中使用非常低精度的方法计算日月坐标仍可以把中天时刻精确到秒。

[古历处理]:

总体上说农历是天文历(天象日历),理想情况下古人计算的农历与今人计算的农历是相同的,但是由于种种原因(当时的天体运行理论局限、计算工具落后、观测精度有限、时间系统不同、动力学理论不完善等原因),古人计算的农历与今人计算的农历有点出入。通常,古人的日历是皇家日历,不能随便更改。所以,我们以古代的农历为标准,重新历算,具体方法详见以下的“古代农历与huangdaojiri.org天文历”章节。有些气朔数据,使用古算法的结果与其相应的颁行历法有冲突,本天文历仍使用颁行历为准,极个别使用古算法的结果(不多,只有公元320年以前几个),当然还存在一些无法考证的数据。

[地标数据]:

城市经纬度数据库做了压缩处理。要扩展经纬数据库,可使用附件中的【城市经纬压缩器】生成代码并替换程序中的JWv数组。程序中现有的经纬度数据本来只有几百个,后来网友郑彬给我一份3000多个城市的数据,比较全面,因此本版使用他的数据。

[时区数据]:

各地时区数据库做了小量压缩处理。如要扩展,可使用附件中的【时区合成器】生成代码并替换程序中的SQv数组。时区数据库是“中华农历网”netghost先生提供。

[年号纪年]:

天文历中的纪年数据参考万国鼎《中国历史纪年表》,有些朝代存在多组帝王,程序中不一定全部纳入。年号纪年问题比较罗嗦,本人史学功底肤浅,无法严格的解决此问题,仅在程序中加入南北朝(约公元420年)至今的年号数据,未能保证其完全正确。在V2.04版本以后,补上了由“中华农历论坛”ldlcau(亦是“国学数典”论坛的ldlcau)先生提供的公元420年以前的年号数据。从V2.04版开始,程序中提供了年号数据管理工具,详见附件中的【纪年管理器】。

[命理八字]:

鼠标经过日历表,显示普通的八字信息。屏幕下方的八字是严格的,以立春时刻界定年,以节气时刻界定月,以本地真太阳时23点界定日(这里的真太阳时使用低精度算法,误差可能达到1秒),以真太阳时23点、1点、3点……界定时辰。注意,这时屏幕上方输入的时间为当地时间(电脑时间,UTC)。

[日数标尺]:

本天文历中使用两种日数标尺。其一是公历,其二是儒略日数。为了简化历算过程,本天文历的历算代码全部使用“儒略日数”作为日数标尺。日历数据的显示则使用公历,它也叫格里历。严格的说,格里历始于1582年10月15日,1582年10月4日及以前的公历是4年1闰的儒略历。公历是在最近一百两年受到西方政治、经济优势的影响才逐渐广泛使用,显然本天文历把“公历”外延到千年以前只是为了建立一个相对统一的连惯的日数标尺,而不是说很久以前全世界就在使用这种历法。再者,“公历”中的儒略历也不完全是古罗马颁行“儒略历”,因为古罗马儒略历曾有过置闰失误。

[开源软件]:

huangdaojiri.org天文历是开放源代码的,我并不打算作任何加密处理。月历生成部分涉及多种历法,而且还要处理纪年、节日信息、古历算等,所以相应的程序比较长。天文算法部分的程序比较短,只要算法已知,编写起来是比较容易的,困难的是如何设计出我们所需高精度、高速的算法。要设计算良好的天文算法,需参考较多资料并作认真的推导、验证等,工作量比较大,所以最好不要更改这部分代码,以免发生错误。

二、古代农历与huangdaojiri.org天文历

平气平朔日推算:只要知道当时采用的平朔望月(或相邻平气)长度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,huangdaojiri.org天文历选择了一种快捷方法,即在现有的历谱数据中提取各历法的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年)

2008年9月5日,我在中华农历网发贴“求:中国古代农历算法”,许多热心的网友给我回复,并提供了大量资料或思想方法。在此感谢各位的同好,正是这些热心朋友的帮助才得以进一步可靠的前推历谱。郑先生(网名zb9003)发合我一本“古代历法计算”pdf书,主要讲三统历、四分历、乾象历,使我对古历计算有了信心。粗略看了以后,发现平气朔法的历理非常简单,只是古人把它算得很复杂,多半在处理分数及公倍数、公约数等初等数论相关的东西。我不主张使用古方法计算,在本软件中,所有古代平气平朔法均使用线性方程来解决问题。记得尤先生(网名“易子”)曾说过“大约在唐代”开始使用平气定朔,围绕这一个要点,准确得知平朔使用的所代范围。

(1)先从徐先生(网名“春光”)发给我的《古代历法摘要(上)》获得“中国历法表”,在这个表中得知中国曾经版发的历法名称及适用范围。春光的一句话给我很大的帮助:“主要参考律历志”。

(2)在网上下载《三千五百年历日天象》pdf书,把古代日历部分传换为bmp格式。通用ocr软件无法准确的将bmp日历转换为文本格式日历数据表,所以利用C++Builder编写一个图片转文字的专用程序,将这本书扫描成文本格式的日历表。在扫描过程中,遇到识别困难的,在识别后的文字前加一个问号。这项工作整整历时5天。误码率为1%,误码的部分通常会有个问号,所以真正误码的不到千分一。当初转化时,只从-100年开始,所以其中的-103年到-101的这3年的日历采用手工输入。

(3)编写javascript程序,对扫描后的数据进行校验,以消除误码并找出原书的错误。校验方法:a.大小月与朔日的日期进行校验;b.1645年以后的一些关键数据还与huangdaojiri.org天文历早期版本比对;c.节气的日期与节气干支比对;d.节气干支与朔日干支校验;e.相邻节气长度粗大误差校验。所有校验不通过的给出报告并进行人工校对,直到不发生错误报告为止。校验过程中,发现原书也有不少错误。

(4)最强有力的校验:分段线性拟合校检!1645年以前的日历采用平气法,也就是说在某一个历法实行期间,相邻节气交接时刻之差(或朔时刻之差)为常数,各交接时刻可以连成直线。由于原书没有给出交接时刻,所以采用了一些变换方法。

设拟合公式为 D=k*n+b,n=0,1,2,...,用D(n)的整数部份与历表比对,如果不能通过则调整k和b并重新比对(k和b的反复调整方法需要算法支持,如果算法适当,只需调整1至4次就可得到最佳的k和b),直到通过,如果不能通过,则说明书中的的数据有问题或古代历算家的计算结果有问题或是扫描出错,那么就减小拟合的日历范围,找出拟合能够通过与不能通过的临界节气,然后人工查看一下原书中该日期到底错在哪里,为什么错。当然找出为什么出错是需要一些技巧,为此我花了好多时间。

分段拟合所用的时间分段范围使用第(1)所述的“中国历法表”。并不是“中国历法表”中所有的历法颁发时段完全与《三千五百年历日天象》吻合,其中的律历志、新唐书等相关的是吻合的,出现这些问题的原因是,国家分裂时期可能有多部历法同时颁行。

(5)古代定气的算法:本天文历使用开普勒椭圆轨道定气,同时考虑了岁差、章动、轨道要素的长期变化。但必竞使用了现代数据,与古代历算参数有点差异,所以程序中对定气朔进行了订正,使之与古历相符。

也许有人会问,为什么不考虑使用古代原来的算法?实际上,为了提高历算精度,清代的《历象考成》已经使用了第谷的宇宙体系,在乾隆七年(公元1742年)以后的《历象考成后编》,开始应用了开普勒行星运动第一、第二定律(但椭圆焦点上的是地球而不是太阳)。因此大约在1734年以后,古代定气与现代计算结果的已经相差无几。这就造成简单的加减乘除根本不能解决问题,所以宁愿采用“现代算法+订正”思路来解决问题。

(6)古代定朔的算法:从公元619年开始,我国就已经使用定朔法,不同年代的定朔算法差异很大。早期使用一些修正表并插值来推算朔日,晚期已经使用牛顿运动定律求解,本人没有时间也没有能力逐一追朔古人的算法,翻阅历代文献,对我来说简直是天书。于是,利用现代推算的月球轨道要素(已考虑了长期变化),把月球运动近似为椭圆进而算出朔日,再与古历比对,生成朔日订正表。接下来,在正式的程序中使用“椭圆运动+朔日订正表”得到精准的古代定朔日历。

反向推算过程2:

(-721年至-104年)

(1)根据《三千五百年历日天象》手工输入秦汉气日朔日表、战国朔日表、春秋朔日表。花费了一个下午的时间。

(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是历法颁行期间朔或气的总个数。

//“朔”直线拟合参数
1457698.231017,29.53067166, // -721-12-17 h=0.00032 古历·春秋
1546082.512234,29.53085106, // -479-12-11 h=0.00053 古历·战国
1640640.735300,29.53060000, // -221-10-31 h=0.01010 古历·秦汉
1642472.151543,29.53085439, // -216-11-04 h=0.00040 古历·秦汉
1683430.509300,29.53086148, // -104-12-25 h=0.00313 汉书·律历志(太初历)平气平朔
1752148.041079,29.53085097, //   85-02-13 h=0.00049 后汉书·律历志(四分历)
1807665.420323,29.53059851, //  237-02-12 h=0.00033 晋书·律历志(景初历)
1883618.114100,29.53060000, //  445-01-24 h=0.00030 宋书·律历志(何承天元嘉历)
1907360.704700,29.53060000, //  510-01-26 h=0.00030 宋书·律历志(祖冲之大明历)
1936596.224900,29.53060000, //  590-02-10 h=0.01010 随书·律历志(开皇历)
1939135.675300,29.53060000, //  597-01-24 h=0.00890 随书·律历志(大业历)
1947168.00//  619-01-21
//“气”直线拟合参数
1640650.479938,15.21842500, // -221-11-09 h=0.01709 古历·秦汉
1642476.703182,15.21874996, // -216-11-09 h=0.01557 古历·秦汉
1683430.515601,15.218750011,// -104-12-25 h=0.01560 汉书·律历志(太初历) 回归年y=365.25000
1752157.640664,15.218749978,//   85-02-23 h=0.01559 后汉书·律历志(四分历)     y=365.25000
1807675.003759,15.218620279,//  237-02-22 h=0.00010 晋书·律历志(景初历)       y=365.24689
1883627.765182,15.218612292,//  445-02-03 h=0.00026 宋书·律历志(何承天元嘉历) y=365.24670
1907369.128100,15.218449176,//  510-02-03 h=0.00027 宋书·律历志(祖冲之大明历) y=365.24278
1936603.140413,15.218425000,//  590-02-17 h=0.00149 随书·律历志(开皇历)       y=365.24220
1939145.524180,15.218466998,//  597-02-03 h=0.00121 随书·律历志(大业历)       y=365.24321
1947180.798300,15.218524844,//  619-02-03 h=0.00052 新唐书·历志(戊寅元历)     y=365.24460
1964362.041824,15.218533526,//  666-02-17 h=0.00059 新唐书·历志(麟德历)       y=365.24480
1987372.340971,15.218513908,//  729-02-16 h=0.00096 新唐书·历志(大衍历,至德历)y=365.24433
1999653.819126,15.218530782,//  762-10-03 h=0.00093 新唐书·历志(五纪历)       y=365.24474
2007445.469786,15.218535181,//  784-02-01 h=0.00059 新唐书·历志(正元历,观象历)y=365.24484
2021324.917146,15.218526248,//  822-02-01 h=0.00022 新唐书·历志(宣明历) y=365.24463
2047257.232342,15.218519654,//  893-01-31 h=0.00015 新唐书·历志(崇玄历) y=365.24447
2070282.898213,15.218425000,//  956-02-16 h=0.00149 旧五代·历志(钦天历) y=365.24220
2073204.872850,15.218515221,//  964-02-16 h=0.00166 宋史·律历志(应天历) y=365.24437
2080144.500926,15.218530782,//  983-02-16 h=0.00093 宋史·律历志(乾元历) y=365.24474
2086703.688963,15.218523776,// 1001-01-31 h=0.00067 宋史·律历志(仪天历,崇天历) y=365.24457
2110033.182763,15.218425000,// 1064-12-15 h=0.00669 宋史·律历志(明天历) y=365.24220
2111190.300888,15.218425000,// 1068-02-15 h=0.00149 宋史·律历志(崇天历) y=365.24220
2113731.271005,15.218515671,// 1075-01-30 h=0.00038 李锐补修(奉元历)     y=365.24438
2120670.840263,15.218425000,// 1094-01-30 h=0.00149 宋史·律历志         y=365.24220
2123973.309063,15.218425000,// 1103-02-14 h=0.00669 李锐补修(占天历)     y=365.24220
2125068.997336,15.218477932,// 1106-02-14 h=0.00056 宋史·律历志(纪元历) y=365.24347
2136026.312633,15.218472436,// 1136-02-14 h=0.00088 宋史·律历志(统元历,乾道历,淳熙历)365.24334
2156099.495538,15.218425000,// 1191-01-29 h=0.00149 宋史·律历志(会元历) y=365.24220
2159021.324663,15.218425000,// 1199-01-29 h=0.00149 宋史·律历志(统天历) y=365.24220
2162308.575254,15.218461742,// 1208-01-30 h=0.00146 宋史·律历志(开禧历) y=365.24308
2178485.706538,15.218425000,// 1252-05-15 h=0.04606 淳祐历               y=365.24220
2178759.662849,15.218445786,// 1253-02-13 h=0.00231 会天历               y=365.24270
2185334.020800,15.218425000,// 1271-02-13 h=0.00520 宋史·律历志(成天历) y=365.24220
2187525.481425,15.218425000,// 1277-02-12 h=0.00520 本天历               y=365.24220
2188621.191481,15.218437484,// 1280-02-13 h=0.00013 元史·历志(郭守敬授时历) y=365.24250
2321919.49// 1645-02-04

值得注意的是:儒略日数的小数部分为0.5才对应晚上0点;新唐书·历志(戊寅元历)公元619年开始使用平气定朔,之前使用平气平朔,1645年及以后使用定气定朔。

举例说明:

求公元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颁行历的朔日与今算得天象不相符的数据:

以下古历部分数据来自《二十史朔闰表》

天象日期 =》 古历日期

天象日期 =》 古历日期

1501-06-16 六月 =》06-151508-01-02 十二 =》01-03
1513-10-29 十月 =》10-281516-10-25 十月 =》10-26
1521-10-01 九月 =》09-301526-07-10 六月 =》07-09
1527-06-29 六月 =》06-281534-06-12 五月 =》06-11
1535-08-29 八月 =》08-281535-10-26 十月 =》10-27
1544-05-22 五月 =》05-211546-01-02 十二 =》01-03
1546-07-28 七月 =》07-271571-08-21 八月 =》08-20
1572-08-09 七月 =》08-081581-10-27 十月 =》10-28
1582-07-20 七月 =》07-191588-04-26 四月 =》04-25
1589-01-16 十二 =》01-171591-09-18 八月 =》09-17
1599-01-26 正月 =》01-271600-02-15 正月 =》02-14
1612-03-02 二月 =》03-031616-05-16 四月 =》05-15
1622-07-09 六月 =》07-081627-09-10 八月 =》09-09
1628-01-06 十二 =》01-071630-04-12 三月 =》04-13
1634-08-24 七月 =》08-231643-02-18 正月 =》02-19
1649-05-12 四月 =》05-111650-11-23 十一 =》11-24
1652-02-09 正月 =》02-101653-09-22 八月 =》09-21
1654-01-18 十二 =》01-191662-02-19 正月 =》02-18
1673-11-08 十月 =》11-091685-02-04 正月 =》02-03
1687-03-14 二月 =》03-131694-06-23 五月 =》06-22
1704-10-28 十月 =》10-291708-02-22 二月 =》02-21
1720-07-06 六月 =》07-051759-03-29 三月 =》03-28
1763-09-08 八月 =》09-071778-03-29 三月 =》03-28
1779-07-14 六月 =》07-131787-12-10 十一 =》12-09
1789-07-23 六月 =》07-221796-06-06 五月 =》06-05
1804-08-06 七月 =》08-051821-06-29 六月 =》06-30
1831-04-13 三月 =》04-121842-01-12 十二 =》01-11
1863-01-20 十二 =》01-191880-11-02 十月 =》11-03
1896-02-14 正月 =》02-131914-11-18 十月 =》11-17
1916-02-04 正月 =》02-031920-11-11 十月 =》11-10

《二十史朔闰表》中有错误的数据:

[书中错误]: 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-27 23:18:45 原著为1821-08-28,当时的测算水平已经很高,误差应小于半小时(一般小于10分钟)。 1588-03-27 10:46:16 原著为1588-03-26,竞然相差半天左右,当时出现如此误差十分费解,原书可能有误。

《三千五百年历日天象》中有错误的数据:

带*表示改后数据与方诗铭《中国史历日和中西历日对照表》相同,带#号表示改后与方诗铭的不同。

[朔的修改]

: -221年以年以前的不易考证的数据 -701年 06.02 改 06.01 # -654年 12.15 改 12.16 # -618年 06.14 改 06.13 * -579年 09.29 改 09.28 * -577年 12.05 改 12.15 原书笔误 -571年 12.27 改 12.28 * -535年 06.26 改 06.25 # -491年 01.22 改 01.23 # -491年 03.22 改 03.23 # -487年 01.08 改 01.09 # -328年 02.20 改 02.19 # -257年 05.13 改 05.14 * -249年 07.14 改 07.13 * -149年 03.21 改 02.21 原书笔误

[朔的错误]

: -26年 正月和二月的大小有误 18年 十、十一月的大小有误 747年 十一月日期错误,应是12.7 1095年 十二月日期错误 1732年 八、十月的大小有误

[古代历书计算出错]

(为了算法统一改了古历): 278年 三月 0410,暂改为 0409 314年 正月 0201,暂改为 0202 314年 三月 0401,暂改为 0402

[同一年代不同历法的差异]

: 270年 四月 0507,暂改为 0506 430年 二月 0311,暂改为 0310

[节气订正]:

-214年 九月份 廿六 改 廿七 -213年 十一月 十五 改 十四 -209年 九月份 八日 改 七日 -174年 六月份 二日 改 上月三十(减2天) -166年 十一月 廿三 改 廿四 -136年 后九月 十九 改 二十 -26年 二月份 十四 改 十五 -26年 二月份 廿九 改 三十 18年 十一月 十三 改 十二 18年 十一月 廿九 改 廿八 28年 九月份 十七 改 十八 44年 十二月 十六 改 十七 81年 九月份 三日 改 四日 89年 五月份 十二 改 十三 91年 十二月 十一 改 十二 134年 二月份 廿八 改 廿九 365年 十二月 九日 改 十日 466年 二月份 十八 改 十九 616年 四月份 十四 改 十三 627年 四月份 二七 改 十七 876年 八月份 十四 改 十日 1115年 十月份 十三 改 十二 1149年 八月份 十三 改 十二 1427年 十月份 九日 改 十日 1719年 十一月 十一 改 十二 1728年 十二月 六日 癸改壬 1732年 九月份 五日 戊改己, 二十 癸改甲 1732年 十月份 五日 戊改己, 二十 癸改甲 1754年 二月份 廿八 壬改戊 698年 十一正 第一个节气减2变为壬日 十月份 第一个节气减1变为壬日 724年 十二月 十三 改 十二 725年 正月份 第一节气减2天变为甲日 1280年 十一月 第一节气减1变甲日 九月份 第一节气减1变壬日 正月份 第一节气减1变己日 (1280年的气本年仅是可能有误)

置闰问题:

从huangdaojiri.org万年历2.0版本(后期改名为huangdaojiri.org天文历)开始,处理古代农历置闰问题时,取消以往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) mod 12,式中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天。

三、《huangdaojiri.org天文历》开发过程

2000年,为了解决公历与农历转换问题,我买了一本王纪卿的《民俗通书万年历》。拿到这本书以后,对里面的节气交接时刻产生了浓厚的举趣。书中介绍了节气交接时刻的简易计算方法,可是我算来算去,忙了一整天就是和书中的节气时刻不相符,误差可达0.5小时。不知为什么,从那以后我没有再研究万年历了。后来,我买的那本《万年历》成了妈妈的工具书,用来查找“吉日”。08年春节,学校网站上的农历显示不正常,于是花了一天时间,利用查表法解决了网站上的农历显示问题,从那以后,我下定决心要弄明白农历到底是怎么回事。

为此,这年春节我花了5天时间,利用“科威尔”数值积分方法,创建了8大行星的数值积分程序,用于求解日月及行星的位置。程序写完了以后,输入某些初始数据进行计算,程序可以在几千年内超高精度的计算。但是,用于计算实际的8大行星时却遇到极大的困难,主要原因是很难找到合理的精确的初始数据。另外,与大地形状等有关的问题也很难解决。为了寻找精确的初始数据,必须有星历表,在网络上查找“星历表”的时候,发现了几种现成的星历计算方法(有些是中华农历网的春光介绍的),于是,花费了1个多月的时间翻译、整理了这些资料,当我翻译完了这些资料,也就基本明白经典天文学在研究些什么,农历计算问题自然也就解决了。

我在网络上找了很久,参阅的天文资料(全部在网络上找的):行星运动VSOP87方法(外文)、月球运动ELP2000-82B(外文)、月球运动ELP/MPP02方法(外文)、行星运动IPS2000方法(外文)、DE系列星历表(外文)、瑞士星历表(外文)、《天文算法》(除了第41、42、43、44章,本书已全部译完,有需要的我可以发一个)、北师大高健的《球面天文学讲议》ppt、《球面三角基本公式》pdf、《有限差分法》doc、《数值方法》、《天体力学引论》(前几年在超星下载打印的)、《IAU1982岁差章动》《IAU2000岁差章动》(外文pdf)、陈垣《二十史朔闰表》pdf、万国鼎《中国历史纪年表》pdf、唐汉良和余宗宽《日月食及其计算概要》pdf、《fortran教程》(现代天文算法一般是用fortran写的)等。要想利用天文算法实现农历计算,这些书籍、资料或多或少是要读的。为了方便验算比对,在农历网又托人购买了《2008年中国天文年历》。

2008年7月7日,正式放假了,我开始至力于万年历的编写,到了8月初,我已花费近25天时间(每天至少9个小时)。汗水没有白费,总算初步完成了。我曾在“中华农历论坛”中发表了一些关于“天文算法”的文章,本天文历的天文算法内核就是使用那里论述的。

在寻找星历表的过程中,无意中找到了“日梭万年历2006测试版”,让我初次领教了“利用天文方法”计算农历的威力。相对于查表法,有几个明显优势:(1)查表法的较长时间跨度的原始数据不易获得,即使有也是海量数据,而天文算法所需的数据量要少得多。(2)精度很高。(3)除错、验证比较方便。正因如此,我很想获得这个程序的原代码。遗憾的是,刘安国先生并没有开放源代码,所以只好自己编写一个,前前后后花费了半年多的业余时间,而且还一直在改进升级之中。虽然程序写好了,错误再所难免,欢迎大家指正。

当我分析了月球的ELP/MPP02的fortran程序之后,就用c++写出相同的程序,并利用这个程序导出javascript所需月球数据及算法,所以本程序所用的算法与原版的ELP/MPP02有所不同:是对数据序列进行转换处理,使得序列的表达形式与VSOP87的序列的形式相同,大大降底了算法的难度,速度也有所提升,缺点是数据增加了20%——30%

四、《huangdaojiri.org天文历》的精度、可靠性及性能

历算精度是一个复杂的问题,以下从多个方面分析历算的精度问题。

世界时:以子午圈作为时间指针(与地球自转同步),以恒星(太阳也是恒星)作为刻度,如此构成了天然的手表。这个手表比我们的石英表准确很多(目前,一年内的误差不超过1秒)。由于地球自转时快时慢,也不太有规律,所以这种时间是不均匀的。对于将来,地球自转有变慢的趋势,这部特殊的手表也将变慢,逐渐落后于原子时。

原子时:是一种随时间均匀计时的高级手表。目前,这种手表在全世界范围内为数不多。

力学时:利用太阳系行星运动规律推导出来的时间,这种时间是均匀。通常,在计算时,原子时与力学时没有太大区别(主要注意一下起算时间即可)。严格上说,力学时还分为地心力学时与太阳系质心力学时,如果我们在定义力学时起点时,有意让起这两个力学时的起算点相同,那么,由于地球公转的相对论效应造成这两种力学时相差的时间不到2毫秒,所以程序中不区分是那一种力学时。

协调世界时(UTC):秒长使用原子时,所以它是力学时。可是,为了与世界时同步,UTC在每年的年底做了一个小动作,即闰秒。所以,从长远看UTC属世界时范畴,在一年内看UTC属力学时范畴。

力学时与世界时之间存在一个差值,叫作ΔT,随着时间的增长ΔT不断变大。在过去的几百年中,ΔT是已知的,但是对于遥远的过去与未来,ΔT是未知的,但我们可以粗略表达为:
t = (y-1820)/100
ΔT = -20+jsd*t2
   式中y是年份(如2200),jsd是ΔT变化的加速度,结果的单位是秒。在瑞士星历表中jsd取值为31,在skymap11中jsd取值为29。这两个软件都是天文大师的精品,但是他们对未来的jsd的估计不尽相同。jsd的估计是凭感觉的,没有严格的科学依据。如果jsd误差为2,1000年(公元3000年)后ΔT的误差为:2*t*t = 2*11.8*11.8 =278秒,3000年后(公元5000年) 2*t*t = 2022秒

2008年以后ΔT的加速度估计值取31(采用瑞士星历表的),如果要使用skymap 11的,请把程序中的jsd改为29。

与LE405/406比较,公元-3000年到公元3000年,时间(力学时)可能发生的理论最大误差小于6秒钟,平均误差约为1秒,由于作了截断处理,实际误差会大一些,实际精度详见星历表中的即时帮助。日月位置精度比上次写的那个javascript精度高5倍,速度也快了3倍以上。

在公元+3000到+5000,拟合瑞士星历表,在这一时间范围内月球黄经与瑞士星历表相差10角秒以内,平均5角秒以内(对应力学时相差约10秒钟)

可见,对于遥远的过去与遥远的将来,历算的误差主要来自ΔT。对于过去的几百年,误差主要来自于日月坐标。对于日历的准确性,主要取决于核对时否认真以及历书著作的可靠性等。

算法的性能:由于程序主体使用javascript,运行速度是很慢的,所以必需考虑高效的算法。要实现高效运算,必须考虑天文算法问题及程序实现的诸多细节,需要一定的程序设计经验,尽管以前设计软件等工作积累下来的经验对这个软件的设计提供了很大的帮助,但对我来说,设计这个软件仍然十分费工夫。在P4/2.4G(FSB 800M)计算机中,本程序的高精度定朔速度约为500个/秒,低精度定朔速度为5000个/秒,高精度定气速度900个秒,底精度定气速度为9000个/秒,这些速度数据可作为基于eph.js和eph0.js的日历设计的参考,当然,后期版再次提高精度,所以速度下降了很多。为了高速有效迭代计算,本人另外推导了一整组日月运动相关的中等或较低精度的数据或方法,它们与高精度算法结合,实现高速计算,测试代码中的“粗定气误差”等就是专门用来检验这些方法的可靠性的。

程序的大小:核心程序大约30K,其它是注释、说明、城市经纬数据、纪年数据、页面等。也就是说,如果做成精简版,只需30几K,后期版本增加了九行星和日月食、地图数据,程序变大了很多,约500k。

月球视坐标比对

V4.13版本及以后 时间 2008.1.6 00:00:00 TT (JED = 2926.5+2451545) 本程序 视黄经 256°54'36.31" 视黄纬 -4°52'14.12" 距离 401817.73千米 DE406 视黄经 256°54'36.319" 视黄纬 -4°52'14.134" 距离 401817.6711 中国 视黄经 256°54' 36.32" 视黄纬: -4°52'14.09" 《2008年中国天文年历》它与JPL不一致 swiss 视黄经: 256°54'36.319" 视黄纬 -4°52'14.090" 距离* 401798.6263 本程序 视赤经 17h 00m58.06s 视赤纬 -27°38'33.64" 距离 401817.73千米 DE406 视赤经 17h 00m58.061s 视赤纬 -27°38'33.691" 距离:401817.6711 中国 视赤经 17h 00m58.061s 视赤纬:-27°38'33.69" 《2008年中国天文年历》它与JPL一致 swiss 视赤经 17h 00m58.061s 视赤纬 -27°38'33.648" 距离*:401798.6263 JPL 视赤经 17h 00m58.0575s视赤纬:-27°38'33.689" 距离*:401798.6270 JPL网站查询(http://ssd.jpl.nasa.gov/horizons.cgi) 比较的最后结果: (1)“ELP/mpp02”、“本人利用DE406计算”、“JPL网站查询”三者结果基本一致。  存在几个毫角秒的差异主要是岁差参数以及黄赤交角参数的选用值不同造成的。  以上三者坐标都在J2000惯性黄道&赤道坐标系中计算并比对。 (2)《中国天文年历》的月亮视赤经、视赤纬与JPL的一致。 (3)《中国天文年历》的月亮视黄纬与JPL不一致,它采用的黄道与J2000惯性黄道存在0.047角秒的夹角。 (4)《2008中国天文年历》的月亮坐标参考系存在前后不统一的现象! (5)swiss的月亮视黄纬、视赤纬与JPL不一致,它采用的黄道与J2000惯性黄道存在0.047角秒的夹角。 (6)swiss与《2008中国天文年历》的月亮视黄纬一致,其黄道与J2000惯性黄道存在0.047角秒的夹角。 (7)上面视坐标中的“距离*”与“距离”的含义是不相同的。设月亮光行时为T,前者是t-T时刻月球  与t时刻地球的距离,后者为t-T时刻地月距。 (8)此刻本程序误差为几十毫角秒

月球黄经计算结果与《2008年中国天文年历》等权威数据的比较

V4.12版本及以前 2008年01月01日0h TD +197°19' 24.43" 中国天文年历 +197°19’24.91" 本程序 2008年01月06日0h TD +256°54' 36.32" 中国天文年历 +256°54' 36.11" 本程序 2008年01月18日0h TD + 56°04' 29.83" 中国天文年历 + 56°04' 29.68" 本程序 2100年01月01日0h TD +157°24' 01.183" 瑞士星历表 +157°24' 01.96" 本程序 2100年01月18日0h TD + 22°14' 39.400" 瑞士星历表 + 22°14' 40.47" 本程序 2200年01月02日0h TD +108°26' 45.916" 瑞士星历表 +108°26' 46.12" 本程序

V5.03以上版本的月亮精度

V5.03版本及以后,月亮星历再提高一倍 1.黄经精度(角秒):0.5 2.黄纬精度(角秒):0.5 3.距离精度(千米):1 4.J2000.0时间最大偏移(世纪数):30 5.小数点增补位数:3 误差检测 1.测试点数N;2月球方位角与高度角 站点坐标: L=-116°23' φ=+ 39°54' 日期时间: 2000年1月1日12h TD(即力学时=0) SkyMap 方位角 348°13' 16" 高度角 -60°58' 19" 本程序 方位角 168°13' 14" 高度角 -60°58' 19"
黄历查询
1234567
891011121314
15161718192021
22232425262728
2930