huangdaojiri.org天文历是一款采用现代天文算法制作的天文、历算程序,可以方便地进行公历、农历、回历三历之间的转换;提供了历谱数据导出功能;提供精确的日月食等天象的计算功能;含有行星、恒星星历表计算。它提供公元-4712年到公元9999年的公历日期查询功能以及-3000年至+3000年的农历查询功能,其中-721年到1960年的农历数据已经与张培瑜的《三千五百年历日天象》、陈垣的《二十史朔闰表》、方诗铭的《中国史历日和中西历日对照表》核对;它还含有公元前2000多年以前到今的基本年号;含有二千多个国内城市的经纬度,并且用户可根据自已的需要扩展经纬度数据。
本软件是一款精准的年代跨度大的日历工具,可作为一般的实用日历工具,对史学家、考古学家、历算工作者、天文爱好者均有较大的参考价值。
如果以力学时作为时间标尺,与DE406星历表比对,太阳坐标的最大可能误差为0.1角秒,月亮坐标的最大可能误差为0.8角秒,平均误差为最大可能误差的1/6,这使得精确的日月食计算成为可能,
huangdaojiri.org还提供日出日没、月出月没时间查询等功能,与国内外著名的天文软件或天文年历比较,结果完全一致。此外,天文历中还即时显示了“月亮被照亮部分比例”、太阳和月亮的地心黄道坐标、地心赤道坐标、站心地平坐标等。其中日月的地心黄道、赤道坐标采用高精度算法,“月亮被照亮部分比例”的计算采用低精度的算法,但可以满足一般需求;地平坐标的计算已适当考虑了大气折射,由于大气的气压及气温变化具有许多不确定的因素,所以只考虑平均情况下的大气折射修正,在本软件中,地平坐标被描述为方位角与高度角,它是以观测点为中心的坐标,所以已经考虑了周日视差修正。屏幕显示的方位角从正南向西测量
。方位角与高度角已转换到是站心地平坐标,并且做了视差修正,同时在地平真纬度大于0时进行大气折射修正。
为了方便爱好者查阅数据,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的值在反推过程中已在历法适用范围内得到完全验证。关于反向推算的算法详见主目录下的《离散序列的直线拟合》章节。
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年开始,我国就已经使用定朔法,不同年代的定朔算法差异很大。早期使用一些修正表并插值来推算朔日,晚期已经使用牛顿运动定律求解,本人没有时间也没有能力逐一追朔古人的算法,翻阅历代文献,对我来说简直是天书。于是,利用现代推算的月球轨道要素(已考虑了长期变化),把月球运动近似为椭圆进而算出朔日,再与古历比对,生成朔日订正表。接下来,在正式的程序中使用“椭圆运动+朔日订正表”得到精准的古代定朔日历。
(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-06-16 六月 =》06-15 | 1508-01-02 十二 =》01-03 |
1513-10-29 十月 =》10-28 | 1516-10-25 十月 =》10-26 |
1521-10-01 九月 =》09-30 | 1526-07-10 六月 =》07-09 |
1527-06-29 六月 =》06-28 | 1534-06-12 五月 =》06-11 |
1535-08-29 八月 =》08-28 | 1535-10-26 十月 =》10-27 |
1544-05-22 五月 =》05-21 | 1546-01-02 十二 =》01-03 |
1546-07-28 七月 =》07-27 | 1571-08-21 八月 =》08-20 |
1572-08-09 七月 =》08-08 | 1581-10-27 十月 =》10-28 |
1582-07-20 七月 =》07-19 | 1588-04-26 四月 =》04-25 |
1589-01-16 十二 =》01-17 | 1591-09-18 八月 =》09-17 |
1599-01-26 正月 =》01-27 | 1600-02-15 正月 =》02-14 |
1612-03-02 二月 =》03-03 | 1616-05-16 四月 =》05-15 |
1622-07-09 六月 =》07-08 | 1627-09-10 八月 =》09-09 |
1628-01-06 十二 =》01-07 | 1630-04-12 三月 =》04-13 |
1634-08-24 七月 =》08-23 | 1643-02-18 正月 =》02-19 |
1649-05-12 四月 =》05-11 | 1650-11-23 十一 =》11-24 |
1652-02-09 正月 =》02-10 | 1653-09-22 八月 =》09-21 |
1654-01-18 十二 =》01-19 | 1662-02-19 正月 =》02-18 |
1673-11-08 十月 =》11-09 | 1685-02-04 正月 =》02-03 |
1687-03-14 二月 =》03-13 | 1694-06-23 五月 =》06-22 |
1704-10-28 十月 =》10-29 | 1708-02-22 二月 =》02-21 |
1720-07-06 六月 =》07-05 | 1759-03-29 三月 =》03-28 |
1763-09-08 八月 =》09-07 | 1778-03-29 三月 =》03-28 |
1779-07-14 六月 =》07-13 | 1787-12-10 十一 =》12-09 |
1789-07-23 六月 =》07-22 | 1796-06-06 五月 =》06-05 |
1804-08-06 七月 =》08-05 | 1821-06-29 六月 =》06-30 |
1831-04-13 三月 =》04-12 | 1842-01-12 十二 =》01-11 |
1863-01-20 十二 =》01-19 | 1880-11-02 十月 =》11-03 |
1896-02-14 正月 =》02-13 | 1914-11-18 十月 =》11-17 |
1916-02-04 正月 =》02-03 | 1920-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,竞然相差半天左右,当时出现如此误差十分费解,原书可能有误。
带*表示改后数据与方诗铭《中国史历日和中西历日对照表》相同,带#号表示改后与方诗铭的不同。
|
|
从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天。
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%
历算精度是一个复杂的问题,以下从多个方面分析历算的精度问题。
世界时:以子午圈作为时间指针(与地球自转同步),以恒星(太阳也是恒星)作为刻度,如此构成了天然的手表。这个手表比我们的石英表准确很多(目前,一年内的误差不超过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"