一、《寿星万年历》概述

寿星万年历是一款采用现代天文算法制作的农历历算程序，含有公历与回历信息，可以很方便的进行公、农、回三历之间的转换。提供公元-4712年到公元9999年的日期查询功能。其中1500年到1940农历数据已经与陈垣的《二十史朔闰表》核对；含有从公420元(南北朝/宋武帝元年)到今的基本年号。在过去几百年中，寿星万年历的误差是非常小的，节气时刻计算及日月合朔时刻的平均误差小于1秒，太阳坐标的最大可能误差为0.2角秒，月亮坐标的最大可能误差为3角秒，平均误差为误差的1/6。万年历中含有几百个国内城市的经纬度，并且用户可根据自已的需要扩展经纬度数据。

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

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

寿星万年历在进行月历推算时采用了最先进的现代天文算法，它是基于2002年巴黎天文台西尔特处Jean CHAPRONT 和 Gerard FRANCOU发表的ELP/MPP02月球运动理论，在几百年范围内，该理论与当今世界上公认的最精密的DE405/406系列星历表仅相差3毫角秒，在几千年范围相差不超过4角秒。由于农历的计算不需要这么高的精度，我们做了一些截断处理及实质性的简化处理，但仍可以满足高精度历算的要求。

程序中采用了1988年法国巴黎天文台Pierre Bretagnon与Gerard Francou发表的VSOP87行星运动理论计算太阳坐标。由于该理论与DE405/406星历表存在少量差异，我们对VSOP87进行少量修正，大大提高了计算精度。

目前的农历是严格的天文历，按理说古人计算的农历应与今人计算的农历是相同的，但是由于种种原因(当时计算工具落后、观测精度有限、时间系统不同、动力学理论不完善)，古人计算的农历与今人计算的农历有点出入。通常，古人的日历是皇家日历，不能随便更改。我们以古人的农历为标准，订正如下：
  [积月:指从公元-4712年1月10日(也可能是9日)那个朔月开始，连续计算的朔月的个数。-4712年1月10日那个朔月的积月是1。]
积月 日期 农历月 =》古历日, 积月 日期 农历月 =》古历日
76852 1501-06-16 六月 =》06-15, 76933 1508-01-02 十二 =》01-03
77005 1513-10-29 十月 =》10-28, 77042 1516-10-25 十月 =》10-26
77103 1521-10-01 九月 =》09-30, 77162 1526-07-10 六月 =》07-09
77174 1527-06-29 六月 =》06-28, 77260 1534-06-12 五月 =》06-11
77275 1535-08-29 八月 =》08-28, 77277 1535-10-26 十月 =》10-27
77383 1544-05-22 五月 =》05-21, 77403 1546-01-02 十二 =》01-03
77410 1546-07-28 七月 =》07-27, 77720 1571-08-21 八月 =》08-20
77732 1572-08-09 七月 =》08-08, 77846 1581-10-27 十月 =》10-28
77855 1582-07-20 七月 =》07-19, 77926 1588-04-26 四月 =》04-25
77935 1589-01-16 十二 =》01-17, 77968 1591-09-18 八月 =》09-17
78059 1599-01-26 正月 =》01-27, 78072 1600-02-15 正月 =》02-14
78221 1612-03-02 二月 =》03-03, 78273 1616-05-16 四月 =》05-15
78349 1622-07-09 六月 =》07-08, 78413 1627-09-10 八月 =》09-09
78417 1628-01-06 十二 =》01-07, 78445 1630-04-12 三月 =》04-13
78499 1634-08-24 七月 =》08-23, 78604 1643-02-18 正月 =》02-19
78681 1649-05-12 四月 =》05-11, 78700 1650-11-23 十一 =》11-24
78715 1652-02-09 正月 =》02-10, 78735 1653-09-22 八月 =》09-21
78739 1654-01-18 十二 =》01-19, 78839 1662-02-19 正月 =》02-18
78984 1673-11-08 十月 =》11-09, 79123 1685-02-04 正月 =》02-03
79149 1687-03-14 二月 =》03-13, 79239 1694-06-23 五月 =》06-22
79367 1704-10-28 十月 =》10-29, 79408 1708-02-22 二月 =》02-21
79561 1720-07-06 六月 =》07-05, 80040 1759-03-29 三月 =》03-28
80095 1763-09-08 八月 =》09-07, 80275 1778-03-29 三月 =》03-28
80291 1779-07-14 六月 =》07-13, 80395 1787-12-10 十一 =》12-09
80415 1789-07-23 六月 =》07-22, 80500 1796-06-06 五月 =》06-05
80601 1804-08-06 七月 =》08-05, 80810 1821-06-29 六月 =》06-30
80931 1831-04-13 三月 =》04-12, 81064 1842-01-12 十二 =》01-11
81324 1863-01-20 十二 =》01-19, 81544 1880-11-02 十月 =》11-03
81733 1896-02-14 正月 =》02-13, 81965 1914-11-18 十月 =》11-17
81980 1916-02-04 正月 =》02-03, 82039 1920-11-11 十月 =》11-10

古历如下月份是闰月(只列出与现代算法不同的部分)，1500-1600年:
  76854,76922,76955,77057,77090,77156,77191,77224,77257,77292,77324,
77392,77425,77459,77493,77526,77559,77628,77661,77727,77762,77795,
77828,77863,77895,77929,77963,77996,78030

古历如下月份是闰月(只列出与现代算法不同的部分)，1600-1700年:
  78098,78130,78198,78231,78265,78333,78366,78399,78466,78500,78533,
78567,78634,78668,78701,78802,78833

古历如下月份是闰月(只列出与现代算法不同的部分)，1700-1800年:
  79645

以上参考书籍为陈垣《二十史朔闰表》，获取勘误数据的输辅助程为[勘误器]

《二十史朔闰表》书中的1500到1940年的农历数据也有错误，核对时请小心：
原书中如下数据肯定有误：
  1808-12-17 原书为1808-12-11误差多达7天，可能是原作者笔误。
  1770-06-23 原书为1770-06-22造成月小于29天。
  1505-12-25 原书为1579-12-23造成月小于29天。
  1579-05-25 原书为05-23造成月小于29天。
书中如下数据可能有误：
  1821-08-27 23:18:45 原书为1821-08-28，当时的测算水平已经很高，误差应小于半小时(一般只有误差10钟以内)，此处的误差达42分钟以上(如果被计算为28日，说明当时算出的朔月时刻应在0点已后)。
  1588-03-27 10:46:16 原书为1588-03-26，竟然相差半天左右，按当时的天文观测水平，出现如此误差的可能性不大，所以原书可能有误。

二、《寿星万年历》开发过程

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

为此，今年春节我花了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年7月7日，正式放假了，我开始至力于万年历的编写，时至今日，我已花费近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*t²
   式中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

与DE405/406比较,公元-3000年到公元3000年，时间(力学时)可能发生的理论最大误差小于6秒钟，平均误差约为1秒。日月位置精度比上次写的那个javascript精度高5倍，速度也快了3倍以上。

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

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

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

程序的大小：核心程序大约30K，其它是注释、说明、城市经纬数据、纪年数据、页面等。也就是说，如果做成精简版，只需30几K

城市经纬度数据库做了压缩处理。要扩展经纬数据库,可使用[压缩器]生成代码并替换程序中的JWv数组。

南北朝至今的年号参考万国鼎《中国历史纪年表》，有些朝代存在多组帝王，则只取其中一组，并制成数据总表。年号纪年问题比较罗嗦，本程序不打算严格的解决此问题。

屏幕显示的方位角从正南向西测量。方位角与高度角已转换到是站心地平坐标，并且做了视差修正，同时在地平真纬度大于0时进行大气折射修正。

节气及合朔时刻问题在天文学层面看来，就是已知天体的坐标反求时间的问题。在VSOP87或ELP中，天体的坐标是时间的函数z = f(t)，所谓的求节气时刻或月相时刻就是已知z求t，显然这是在求解一个关于t的方程。伟大的英国天文学家物理学家牛顿给出了一种非常有效的迭代算法：牛顿求根法。用这种方法，求t所花费的时间仅是求f(t)花费时间的1.2——1.3倍，并且迭代引起的额外时间误差只零点零几秒。牛顿求根法在一般的《数值方法》书籍中均有介绍。

月球黄经计算结果与《2008年中国天文年历》等权威数据的比较
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"  本程序

月球方位角与高度角
站点坐标: L=-116°23' φ=+ 39°54'
日期时间: 2000年1月1日12h TD(即力学时=0)
SkyMap  方位角 348°13' 16"  高度角 -60°58' 19"  
本年历  方位角 168°13' 14"  高度角 - 60°58' 19"

三、版权问题

你可以随意使用本程序中除天文算法之外的任意部分代码，如有异义，可与我共内探讨。xunmeng04@163.com，13850262218。 如果你对天文学不太了解，请不要对天文计算模块做任何修改，以免弄巧成拙。

作者：许剑伟，2008年8月31日于家里