sm1.htm

寿星万年历是一款采用现代天文算法制作的农历历算程序

 698年 十一正 第一个节气减2变为壬日
       十月份 第一个节气减1变为壬日
 724年 十二月 十三 改 十二
 725年 正月份 第一节气减2天变为甲日
1280年 十一月 第一节气减1变甲日
       九月份 第一节气减1变壬日
       正月份 第一节气减1变己日

(1280年的气本年仅是可能有误)
</pre></td>
</tr></table>


<p><b>置闰问题：</b></p>
<p>从寿星万年历2.0版本开始，处理古代农历置闰问题时，取消以往1.x版本的置闰修正算法，而是在准确算出古代气朔表的基础上统一使用“无中置闰法”(-103年至今)或“19年7闰法”(-721年至-104年)得到。</p>
<p><font color=red>十九年七闰法</font>：7个闰年均均安插在19个年整数中，闰年的末月置为闰月。</p>
<p><font color=red>平气无中气日置闰法</font>：不含平中气日的月份置为闰月。这样，平气平朔历法中，7闰月均匀安插在235个月中。
<p><font color=red>定气无中气日置闰法</font>：使用定气法之后，把包含冬至日的月份定为一年中的首月。定年首的后果是有的年份有13个月，有的只有12个月。在含有13个月的年份中，采用无中气日置闰法，并且只对首个无中气的月份置闰。</p>
<p>“定气无中置闰”兼容“平气无中置闰”，反之不兼容。也就是说，平气无中置闰也可加入“含有13个月”这个条件，不会影响历算结果。因为，在平气历法中不可能出现一个月包含两个中气的现象，被置闰的年份中必含有13个月。这使算法变得比较统一，降底了程序的复杂性。<b><font color=red>本万年历严格执行上述置闰法，得到与实闰相同结果。</font></b></p>

<p><b>月建问题：</b></p>
<p>史记·天官书：“值斗杓所指，以建时节”。在古代，常以北斗七星的斗柄方向判断四季，历法与天象存在直接的联系，这里的“建时节”可用“建何月”进行历书表达。换句话说，斗柄指向某一角度，当前月定为“建某月”。往后的历法以标准的太阳历定四季(以12中气为准，斗柄只是参考)，即太阳黄经达到某一位置角度(中气定义的角度)，当前月定为“建某月”。</p>
<p>古历法中，农历1个月中含冬至的是建子月，含大寒的建丑月，其余类推。没有中气的月分，就没有独立月建。因此，在今人进行农历的历月计算时，月建实际上就是一个冬至起算的中气序数(不含中气的农历月没有独立的序数，其序数与前一个月相同)。</p>
<p><b>月建与农历月名称：</b></p>
<p>我们不妨把“正月、二月、三月等”看作月建的别名。并且建子为十一，建丑为十二，建寅为正……由于某些原因，有时古代帝王随意更改了别名。本万年历也相应作了更改：<br><br>
-721年12月17日至-479年11月春秋战国历采用19年7闰，闰年的末月置闰并取名闰十三。年首为正月。<br>
-221年10月31日至-104年11月秦汉历同样采用19年7闰，闰年的末月置闰并取名后九月。年首为十月。<br>
以上两行所述的年代，其月建别名计算方法：<br>
　　<i>月建序数 = (月积累数 - 置闰积累数 + C) mod 12，式中C为某一常数。</i><br>
<br>
008年01月15日至023年12月02日：建子为十二，建丑为正月，建寅为二月，其它顺推。<br>
237年04月12日至239年12月13日：建子为十二，建丑为正月，建寅为二月，其它顺推。<br>
689年12月18日至701年01月14日：建子为正，建寅为一月，其它不变。<br>
761年12月02日至762年03月30日：建子为正，其它顺推。</p>

<p>由于上述月建别名的更改，造成排算的结果是239年12.13为十二月，240年01月12日变回原来的月建别名，这样240年01月12日建丑十二月，显然连续的出现两个十二月，本万年历中把上个月表示为“拾贰月”。同样的问题也出现在23年底。</p>
<p>237年正月开始使用新的历法，造成前一月只有28天。</p>

<p><b>三、《寿星万年历》开发过程</b></p>
<p>2000年，为了解决公历与农历转换问题，我买了一本王纪卿的《民俗通书万年历》。拿到这本书以后，对里面的节气交接时刻产生了浓厚的举趣。书中介绍了节气交接时刻的简易计算方法，可是我算来算去，忙了一整天就是和书中的节气时刻不相符，误差可达0.5小时。不知为什么，从那以后我没有再研究万年历了。后来，我买的那本《万年历》成了妈妈的工具书，用来查找“吉日”。今年春节，学校网站上的农历显示不正常，于是花了一天时间，利用查表法解决了网站上的农历显示问题，从这以后，我下定决心要弄明白农历到底是怎么回事。</p>
<p>为此，今年春节我花了5天时间，利用“科威尔”数值积分方法，创建了8大行星的数值积分程序，用于求解日月及行星的位置。程序写完了以后，输入某些初始数据进行计算，程序可以在几千年内超高精度的计算。但是，用于计算实际的8大行星时却遇到极大的困难，主要原因是很难找到合理的精确的初始数据。另外，与大地形状等有关的问题也很难解决。为了寻找精确的初始数据，必须有星历表，在网络上查找“星历表”的时候，发现了几种现成的星历计算方法(有些是中华农历网的春光介绍的)，于是，花费了1个多月的时间翻译、整理了这些资料，当我翻译完了这些资料，也就基本明白经典天文学在研究些什么，农历计算问题自然也就解决了。</p>
<p>我在网络上拼命找，参阅的天文资料(全部在网络上找的)：行星运动VSOP87方法(外文)、月球运动ELP2000-82B(外文)、月球运动ELP/MPP02方法(外文)、行星运动IPS2000方法(外文)、DE系列星历表(外文)、瑞士星历表(外文)、《天文算法》(除了第41、42、43、44章，本书已全部译完，有需要的我可以发一个)、北师大高健的《球面天文学讲议》ppt、《球面三角基本公式》pdf、《有限差分法》doc、《数值方法》、《天体力学引论》(前几年在超星下载打印的)、《IAU1982岁差章动》《IAU2000岁差章动》(外文pdf)、陈垣《二十史朔闰表》pdf、万国鼎《中国历史纪年表》pdf、唐汉良和余宗宽《日月食及其计算慨要》pdf、《fortran教程》(现代天文算法一般是用fortran写的)等。要想利用天文算法实现农历计算，这些书籍、资料或多或少是要读的。</p>
<p>2008年7月7日，正式放假了，我开始至力于万年历的编写，时至今日，我已花费近25天时间(每天至少9个小时)。汗水没有白费，总算初步完成了。我曾在“中华农历论坛”中发表了一些关于“天文算法”的文章，本万年历的天文算法内核就是使用那里论述的。</p>
<p>在寻找星历表的过程中，无意中找到了“日梭万年历2006测试版”，让我初次领教了“利用天文方法”计算农历的威力。相对于查表法，有几个明显优势：(1)查表法的较长时间跨度的原始数据不易获得，即使有也是海量数据，而天文算法所需的数据量要少得多。(2)精度很高。(3)除错、验证比较方便。正因如此，我很想获得这个程序的原代码。遗憾的是，刘安国先生并没有开放源代码，所以只好自己编写一个，前前后后花费了半年多的业余时间。虽然程序写好了，错误再所难免，欢迎大家指正。</p>
<p>当我分析了月球的ELP/MPP02的fortran程序之后，就用c++写出相同的程序，并利用这个程序导出javascript所需月球数据及算法，所以本程序所用的算法与原版的ELP/MPP02有所不同：是对数据序列进行转换处理，使得序列的表达形式与VSOP87的序列的形式相同，大大降底了算法的难度，速度也有所提升，缺点是数据增加了20%——30%</p>

<p><b>四、《寿星万年历》的精度、可靠性及性能</b></p>
<p>历算精度是一个复杂的问题，以下从多个方面分析历算的精度问题。</p>
<p>世界时：以子午圈作为时间指针(与地球自转同步)，以恒星(太阳也是恒星)作为刻度，如此构成了天然的手表。这个手表比我们的石英表准确很多(目前,一年内的误差不超过1秒)。由于地球自转时快时慢，也不太有规律，所以这种时间是不均匀的。对于将来，地球自转有变慢的趋势，这部特殊的手表也将变慢，逐渐落后于原子时。</p>
<p>原子时：是一种随时间均匀计时的高级手表。目前，这种手表在全世界范围内只有数台。</p>
<p>力学时：利用太阳系行星运动规律推导出来的时间，这种时间是均匀。通常，在计算时，原子时与力学时没有太大区别(主要注意一下起算时间即可)。严格上说，力学时还分为地心力学时与太阳系质心力学时，如果我们在定义力学时起点时，有意让起这两个力学时的起算点相同，那么，由于地球公转的相对论效应造成这两种力学时相差的时间不到2毫秒，所以程序中不区分是那一种力学时。</p>
<p>协调世界时(UTC)：秒长使用原子时，所以它是力学时。可是，为了与世界时同步，UTC在每年的年底做了一个小动作，即闰秒。所以，从长远看UTC属世界时范畴，在一年内看UTC属力学时范畴。</p>

<p>力学时与世界时之间存在一个差值，叫作ΔT，随着时间的增长ΔT不断变大。在过去的几百年中，ΔT是已知的，但是对于遥远的过去与未来，ΔT是未知的，但我们可以粗略表达为：<br>
 t = (y-1820)/100<br>
 ΔT = -20+jsd*t<sup>2</sup><br>
&nbsp;&nbsp;  式中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秒
</p>

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

<p>与DE405/406比较,公元-3000年到公元3000年，时间(力学时)可能发生的理论最大误差小于6秒钟，平均误差约为1秒。日月位置精度比上次写的那个javascript精度高5倍，速度也快了3倍以上。</p>
<p>在公元+3000到+5000，拟合瑞士星历表,在这一时间范围内月球黄经与瑞士星历表相差10角秒以内,平均5角秒以内(对应力学时相差约10秒钟)</p>
<p>可见，对于遥远的过去与遥远的将来，历算的误差主要来自ΔT。对于过去的几百年，误差主要来自于日月坐标。对于日历的准确性，主要取决于核对时否认真以及历书著作的可靠性等。</p>

<p>算法的性能：由于程序主体使用javascript，运行速度是很慢的，所以必需考虑高效的算法。要实现高效运算，必须考虑天文算法问题及程序实现的诸多细节，需要一定的程序设计经验，尽管以前设计软件等工作积累下来的经验对本万年历的设计提供了很大的帮助，但对我来说，设计这个软件仍然十分费工夫。在P4/2.4G(FSB 800M)计算机中，本程序的高精度定朔速度约为500个/秒，底精度定朔速度为5000个/秒，高精度定气速度900个秒，底精度定气速度为9000个/秒，这些速度数据可作为基于eph.js的日历设计的参考。为了高速有效迭代计算，本人另外推导了一整组日月运动相关的中等或较低精度的数据或方法，它们与高精度算法结合，实现高速计算，测试代码中的“粗定气误差”等就是专门用来检验这些方法的可靠性的。</p>
<p>程序的大小：核心程序大约30K，其它是注释、说明、城市经纬数据、纪年数据、页面等。也就是说，如果做成精简版，只需30几K</p>


<pre>
<b>月球黄经计算结果与《2008年中国天文年历》等权威数据的比较</b>
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"  本程序
</pre>
<pre>
<b>月球方位角与高度角</b>
站点坐标: L=-116°23' φ=+ 39°54'
日期时间: 2000年1月1日12h TD(即力学时=0)
SkyMap  方位角 348°13' 16"  高度角 -60°58' 19"  
本程序  方位角 168°13' 14"  高度角 -60°58' 19"
</pre>

<p><b>五、版权问题</b><br></p>
<p>本程序是开源的，你可以使用其中的任意部分代码，但不得随意修改“天文算法(eph.js)”及“农历算法(lunar.js)中古历部分的数据及算法”。一旦修改可能影响万年历的准确性，如果你对天文学不太了解而仅凭对历法的热情，请不要对此做任何修改，以免弄巧成拙。</p>
<p>如果在你自己开发的软件中使用了本程序的核心算法及数据，你可以在你的软件中申明“数据或算法来源于寿星万年历”，也可以不申明，但不可以申明为它其它来源。如有异义，可与我共内探讨。</p>
<p>作者：许剑伟，2008年11月于家里。[e-mail]，13850262218</p>


</td></tr>
</table>
</center>
</body></html>