几个月前增编写个一个农曆计算javascript程序，那是我的第一个农历程序为了实现其中的算法，曾花费了好几个星期的时间研究天文计算相关的原理当我算出结果之后，对程序的结果仍然没有信心不过，在“春光”老师介绍《天文算法》一书之后我认真阅读并翻译全书之后，我知道我的算法基本仩没有错误，但同时也认识到天文学家们的算法的确高明，比我当时的算法好得多如果你也对天文计算感兴趣，那就很有必要学习他們的先进思想方法

　　当然，重要的是他们的思想方法以及相关的理论而不一定是他们的计算结果。因为《天文算法》一书是早期絀版的，有些数据比较旧造成精度不一定很好。如果你对计算精度要球特别高可能需要更换一些数据，但数据的处理的思想方法及技巧是一样的

　　国内不少网友对天文计算比较感兴趣，却又不知如何下手问题出在哪里呢？为什么苦苦研究几个月甚至几年也没有进展主要原因是国内有关的书籍太少，业余条件下没有机会学习到天文计算的理论在我们国家，有不少机构对天文学有深入的研究比洳天文台、一些大学等，可他们不太愿意出版天文算法之类的书籍(这类书籍销量少出版要赔钱的)。既然如此就让我们自已想办法解决問题吧。

　　首先我认为需要掌握一定的计算机程序设计技术，不要求很利害但起码也要有几个月的程序设计经验。如果你不会程序設计那你只能用Excel或计算器之类的工具来处理计算问题，那简直是在浪费你的生命

　　其次，应掌握《高等数学》中的一些知识比如：极限、导数、微分、积分、极值问题、求根问题、最小二乘法、向量数学等。当然我们更多使用高中的《立体几何》、《解析几何》、《函数》、《三角函数》等有关知识。还应了解《球面三角学》里的几个公式

　　其三，《数值方法》这类书籍是必须读的

　　其㈣，需对物理学有所了解尤其是运动学相关知识。当然如果你想用数值积分的方法解决天文计算问题，《理论力学》甚至是《天体力學》也是有必要了解的

　　对于多数具有大学学历的人来说，基本具备以上知识也就是说，只要你有兴趣就完全可以进行天文计算。从本质上讲日月运动、行星运动主要使用“牛顿力学”及数学方法(如微积分)，在牛顿那个年代力学理论、数学理论、计算工具等都鈈可能和现在相比，在那个年代就连“除法问题”主要是“教授”们才能掌握的！他们可以计算天体运动问题，我们为什么就不可以呢只要有信心，或多或少可以解决问题

　　相反，如果根本没有学过《高中数学》《高中物理》《高等数学》我建议你还是花点时间學习一下(弄不好要花费一两年时间)，否则即使算出了结果，也很难对你算出的结果形成理性的认识甚至感性的认识也谈不上。如果你鈈想学习这些“无用”的东西能不能实现天文算法，或许可以：通过阅读别人的程序

　　在太阳系中，我们一致认为九大行星(八个主荇星和1个冥王星)围绕太阳系运动几百年前，哥白尼先生就已经这么认为！现在的教科书中也都是这么说的！因为太阳系中的星体围绕太陽运动所以我们常以太阳为中心来考虑问题。

　　以太阳为中心考虑问题有不少好处：我们可以把行星围绕太阳运动近似看做圆周运动(圓周运动是物理学中的一个重要名词)这样，行星的运动显得非常有规律

　　如果不以太阳为中心，而以地球为中心我们会发现水星、金星的运动毫无规律，你可以充分发挥你的空间想象你将发现，水星和金星在天球上的视运动不可能是围绕地球做圆周运动而是在呔阳附近“摆动”。

　　然而我们不可能站在太阳上面进行天文观测，通常是站在地球上进行天文观测的所以“地心说”也是非常科學的，也就是说我们可以认为，太阳系的中心是地球

　　不管是“日心说”、“地心说”，只要你科学对待不要随意附加“神学的”或“迷信的”论调，都可以正确解天文现象

　　要记住，这个世界上没有绝对的“中心”太阳个头大，质量大所以行星围绕它运轉。仅因它质量大就能决定一切吗？显然不是我们可以这么想，某件事是大人说了算还是小孩说了算，通常是大人说了算但具体問题要具体分析，有时候小孩说了算

　　从纯牛顿力学的角度看，不管是“日心”还是“地心”通过计算，均能得到相同的数值计算結果

　　因此，根据实际需要有时我们使用“地心”坐标，有时候使用“日心”坐标

　　为了计算方便以及更好的解释天文现象，囚们还引入了“太阳系质心”、“银心”、地月系质心、月心、行星的“心”、纯几何的视中心等等，一系列的“中心”不管是什么“中心”，我们必须正确理解这些“中心”的具体定义必须明白这些“心”到底在哪里，以便进行坐标变换如果知道这些心的定义，峩们未必原原本本的使用现成的天文坐标变换公式也可适当根据程序设计的需要建立自已所需的坐标变换公式(比如在近似处理时)，使程序更加简洁明了

　　茫茫宇宙，无边无际——宇宙大得让我们难以想象宇宙是球形的吗？我们很难说清楚这个问题但是，进入我们視野的天空就象一个“球”，我们称它为天球我们位于天球的中心。这个球的半径有多大半径很大很大，看成无穷大也无妨太阳系内星体之间的距离根本不能与天球半径相比，如果把天球半径看为1那么太阳系内星体间的距离可以看做0。在我们的视野中所有的恒煋都在天球上。

　　赤道的概念已被大家熟知现在，我们扩展一下赤道的概念赤道平面扩展到天球，与天球相交交线为天球上的一個大圆，这个圆称为天赤道

　　不管以太阳为中心，还是以地球为中心我们看到的天赤道是相同的，所有恒星与天赤道的相对视位置關系是一样的为什么呢？因为从几何学看来，太阳系的尺度太小不能与天球半径相比，在太阳系中的任何一个地方都可以看做天浗的中心，所以恒星之间的视位置关系以及它们与天赤道之间的视觉位置关系不会因为中心的改变而发生改变。

　　当然恒星并不是嫃的在天球的表面上，个别恒星与我们之间的距离不是很远所以，当观测点(即“中心”)改变时它们的视位置也会有点变化，这就是我們常说的恒星周年视差我们以后再讨论这些问题吧。

　　地球公转的轨道面同样可以扩展到天球，得到的交线也是个大圆称为黄道。

　　虽然赤道与天赤道是不同的两个概念但具体问题所涉及的是“天赤道”还是“赤道”往往是非常明确的(明显的)，一般不会产生二意所以通常“天赤道”也简称为“赤道”。

　　天赤道与黄道必然在天球上产生两个交点一个称为“春分点”(太阳经过升交点)，一个昰“秋分点”(太阳经过降交点)“春分点”是天球上重要的参考点。我们描述天体位置通常是相对于这个参考点而言的。

　　严格上说黄道是地月质心公转轨道在天球上的扩展。

　　在天球上并没有一个恒星可以直接用来标定春分点的位置，它是看不见也摸不着的峩们如何确定它的位置？最重要的一种方法是：通过数以千计的恒星位置反推出春风点在天球上的位置，我们常说的FK5天球坐标系统就与咜有关还有其它方法可以确定春分点，比如动力学方法由于所用的方法不同，得到的春风点也会有一些微小的差别难道有很多种春風点吗？其实“春风点”这个名词最好不要乱用，如果说FK5中的那个“点”是真正的春分点那么动力学方法(如VSOP87)得到的那个点就不应该称莋春分点，其实在《天文算法》一书中很少用“春风点”一词取代之“分点”一词。当然我们不一定计较那么一点点微小的差别，统稱为“春分点”也不要紧

　　地平线——天地的交线，它又是一个圆圆中心位于观测点。考虑到天球很大我们可以把圆(以及圆中心)岼移到地心或日心，这样黄道、赤道、地平线的中心就相同了。也许你会问平移后，会不会影响日月在地平坐标中的位置当然会有┅点影响，不过我们可以通过视差修正的方法来解决问题

　　当我们在天球上标定了赤道、黄道、分点之后，我们接着来看日、地在黄噵坐标中的运动以太阳系中的任意一点为天球的中心都不影响分点及黄道位置，当天球中心移到地心地心直角坐标变为零，没有经纬喥但可以得到太阳的经纬度，此时如果跑到日心看地球，可以得到地球的经纬度这两个经纬度正好相反(日地连成的直线与天球相交嘚两点正好完全相反)：

　　地球经度 = 太阳经度+180度，地球纬度 = -太阳纬度

　　历元、摄动、自球自转轴、天北极、黄极、黄轴、黄经、黄纬、赤经、赤纬，这几个概念不再讲述请在网络上找一些资料看看吧！

　　春分点是黄道与赤道的交点！

　　春分点是天球中的重要的参栲点，遗憾的是春分点会移动！

　　在图中，有两条黄道、两道赤道

　　一条的历元2000.0时刻赤道，一条是当日(Date)赤道

　　随时间推移，黄道为什么会发生变化因为地球受至大行星的摄动。

　　随时间推移赤道也会发生移动，因为受月球影响地球自转轴的方向发生移动，造成天北极围绕黄轴缓慢移动找一个直铁棒，铁棒放置的方向与黄轴相同用电焊机紦铁棒焊在赤道面的中心。始终保持铁棒与黄轴平行并转动铁棒，这时你将看到北极围绕铁棒做圆周运动，同时春分点也同步移动當北极运行一周(约26000年)，春分点在黄道上也移动一周

　　春分点移动的后果：黄经、赤经的起算点是春分点，造成所有星体的黄经、赤经發生变化假设，黄道不动(没有发生变化)赤道发生变化，引起春分点移动如果春分点在黄道上向西移动了1度(在上图中，往左下方移动)那么，所有星体的黄经增加了1度！

　　春分点的移动量称为岁差岁差有很多种表达方法，最著名、最传统的岁差是指黄道上的岁差

　　上图中，有三个重要的交点：γoγ1，γm我们统称它为分点。γo称为“历元2000.0标准分点”其实就是历元2000时刻的春风点，γm称为“Date分點”或“当日分点”或“当日春分点”

　　图中，两个黄道的交点是N定义N到γo与N到γo'的距离相等，那么以γo'起算的黄经则与γo起算的黃经是相等的(除非天体过份靠近黄极)对于行星，它在黄道附近对于大部分恒星，也是远离黄极的所以γo'对于行星、月球、大部分恒煋来说，是很重要的一个参考点比如太阳，当它从γo出发转动数周后达到了γo'，所经历的时间是整倍数的太阳绕地球运转的周期这種运行周期具有比较严格的力学意义。恒星在天空中的位置几乎是不变的以它作为参考系，具有“惯性”系的特征γo则是利用恒星确萣的，也具有惯性特点任意时刻，我们利用某一颗或几颗恒星的位置轻松的标定出γo的位置(当然，首次建立FK5系统时需要大量恒星标萣出γo，以提高精度所以工程巨大)，当太阳运行到γo'可以认为太阳相对于惯性参考点γo移动了一周或数个1周，即相对于惯性参考点移動了数个周期这就是力学意义上的周期。γo与γo'相差很小所以太阳到了γo与太阳到了γo'，太阳在星空的中相对位置几乎相同因此，呔阳经过这两点所经历的时间是整倍数的恒星年(以恒星为参考太阳在天球中运转一周的时间)，由于这个原因恒星年一般认为就是太阳(戓地球)在惯性系中的运动周期，“在惯性系中的运动周期”常称为真周期

　　p，γo'到γm的角距离就是Date黄道上的岁差
　　ψ，γo到γ1的角距离是J2000黄道上的岁差
　　χ，γ1到γm是Date赤道上的岁差
　　εo是历元2000.0的黄赤交角
　　ε， 是当日的黄赤交角
　　ω， Date平赤道与J2000黄道的夹角
　　π，是两个黄道之前的夹角
　　П，是N到γo的角距离，它等于N到γo'的角距离
　　另一组常用的岁差计算参数
　　90°-ζ，γo到Q的角距离是一个常用的岁差计算参数
　　90°+z， γm到Q的角距离是一个常用的岁差计算参数
　　θ是两个赤道面的夹角

　　也许你会问，怎么会有這么多岁差参数有必要吗？回答为了坐标变换方便而已。通常只需只道3个岁数参数，就可利用坐标变换方法及多项式插值法导出其咜的岁差参数

　　一组常用的岁差表达，取自IAU2000

　　以上岁差表达是非常精确的误差小于1毫角秒。

　　不过t不能太大t应小于1。

　　如果t=10那么，会有什么后果比如计算θ，t的7次方项为0.00=100角秒！实际上，表达式系数的舍入误本身就是0.00001所以误差高达100角秒。

　　如果由于某種需要你非要计算t=10的情况，最好把t的5、6、7次方项忽略这样，如果t=5误差也可控制在1角秒之内。

　　顺便说一下儒略千年数是一种时間计数方法。日常生活中时间的计数法使用“年月日时分秒”，在天文计算过程这种时间表达是很不方便的，原因是明显的需要6个數才能表达一个时间。

　　天文学中采用2000年1月1日12时0分0秒时作为计时的起点，并连续计时单位是天。1天日86400秒1秒是多少？1秒时间长度由原子钟得到现代天文学的时间基础已不再使用星体运动位置来确定，改用原子钟原子钟的精度要好得多。

　　为什么又称为“儒略”芉年数这还要谈到另一个问题，一年有多少天你可能立刻回答到：365.2422日。然而在现代天文学家眼里，一年的长度是变化的并不是固萣为365.2422日，何况地球运动时还受到各种干扰365.2422只是一定时期内的平均值。因此任意一年某时刻加上365.2422后，并不是下一年的同一时刻这样，紦365.2422看做1个单位意义不大天文学家更原意把365.25看做1个单位，它与公历的置闰规则具有更好的“兼容”性便于日历转换。365.25天是儒略年的长度儒略1千年就是365250天。

　　再次强调1天=86400原子秒，与地球自转没有必然的关系在历史上，曾使用某一年的地球公转与自转来定义秒长之後，地球自转速度发生改变1天不等于86400秒了。注意此“天”非彼“天”，前面的天是时间单位后面的天指平均自转一周所用的时间长喥。然而人类的日常生活与地球自转息息相关，日期与时间的计数必须尽可能与地球自转同步这就造成，我们有两套基本的时间系统一个是手表时(或称为协调世界时或称为UTC时)，一个是力学时力学时又时什么东西呢？力学时具有鲜明的动力学色彩假如我们使用2000年地浗公转一圈所用时间定义为1个单位，2万年后公转周期可能发生了变化，如果不通过力学原理推算我们将无法得到那时候一年的长度，泹是通过力学方法，我们可以得到2万年后公转的周期比如是1.00001(乱写一个)个单位。简单的说当我们定义了时间单位，我们就可以使用力學原理推算天体在某一位置所对应时间这就够成时间基础(根据星体位置确定时间)。原子时就简单多了它不是通过天体位置来确定时间嘚，而是通过原子钟内部的振荡次数来来确定时间原子钟曾与力学时对准过，秒长是相同的起点相差零点几秒钟。

　　理论上力学時、原子时是非常均匀的。手表时与地球自转(以恒星为参考点确定自转速度)同步而地球自转时快时慢，所以不均匀

　　手表时(去除地區的时差后的手表时,即通过跳秒与UT1同步的时间)与原子时的差称为ΔT，当然手表时含有时区信息如北京时间多了正好8个时。ΔT可由一组多項式表达式计算得到（UT1同由自转与恒星(春风点)确定的，是天文意义的时间系统严格说，UT1与时角有关即与角度有关；原子时则与“振蕩次数”有关；手表时UTC的秒长与时角无关,它与原子时的秒长完全相同。不过UTC通过年终的跳秒实现与UT1同步，这样就造成UTC与原子时不同步）

　　手表时的秒长是原子秒，为什么手表时不是原子时呢因为国际上有人动了手脚，在每年的最后一天做了“正或负跳秒(正负闰秒)”處理实现手表时与地球自转同步。

我们在第(五)贴已经谈到了赤道与黄道的变化将引起岁差受到月球运动的影响，地球自转轴发生围绕黃轴的缓慢的长期的圆周运动造成赤道发生变化，进而引起春风点变化从上面的那个岁差示意图中不难发现，春风点的变化不单单是赤道变化引起的黄道的变化也会引起春风点的变化。不过图中的角π非常小，所以黄道变化引起和春风点移动是非常小的。

行星光行差计算起来比较麻烦，也没有现成的简单公式可以应用我们总是根据具体的行星问题展开计算。用下文所述的方法得到的光行差是非常准确的最后得到的行星视位置与skymap的计算结果分毫不差。当然日月计算也会涉及行星光行差问题，但这只以行星光行差中的最简单的一個特例

t时刻，地球在A点天体在B点。此时我们在地球上看到的星光并不是并非t时刻天体发出来的光。由于光速是有限值星光传到地浗需要一定的时间τ，所以t时刻我们看到的星光是天体在t-τ时刻发出来的。在图中，O和S分别表示t-τ时刻地球和天体的直位置。我们称τ为光荇时。

可以把τ看作一个时间单位，图中矢量的既可表示速度大小和方向，也可表示位移的大小和方向。

把t时刻看到的光线（图中的SA矢量）可分角为SO矢和SC矢。显然SC矢与地球的速度(或位移)OA是相同的这样，我们站在地心看光波感觉不到SC分量，光的速度只剩下SO分量（这好比兩个人并行走路感觉对方是相对自已是静止的）。SO矢量与CA矢量是相同的考虑到光线最终落到A点，所以星体的视觉方向为AC因此，视方姠与真方向的差为图中的角CAB以上推理是基于伽俐略变换的，我们总结如下：

t时刻星体的视方向是由t-τ时刻星体发出的并传到观测者的光线决定的，在地心坐标中，视方向由该光线的SO分量决定换句话说，在地心坐标中t时刻星体的视方向是t-τ是刻星体的真位置方向。

我们萣义：t时刻，行星光行差等于视方向与真方向的角度差这个角度差就是图中的角CAB。通常角度差（行星光行差）用经纬度修正量来表度。反过来说真方向加上行星光行差就得到视位置。

我们已经知道行星光行差与t-τ=时刻的真方向有关，由于光行时τ是未知的，这给计算带来一点小麻烦。不过，多数情况下，B点的地心坐标(经度U、纬度V、距离Δ)已经精确算出那么我们可以用τ=Δ/c来估算光行时，式中c为光速Δ的单位中千米，τ的单位是秒。

接下来我们有两种计算方法求解修正量：

（1）使用低精度算法重新计算t-τ时刻与t时刻天体的地心坐標。然后算出这两个坐标的经纬度差值即可这个差值就是行星光行差。应注意用t-τ的坐标值减t时刻的坐标值。还应注意，应使用完全相同的算法算出这两个坐标，并且计算过程中请保留所有的数字(16位)，不要做任何舍入以保证能够分辨出差值。

刚才说到“低精度”算法具体含义是什么？我们的目标是准确计算出位置的差值实际上是差值的精度达到3至4位有效数字即可。因为光行星光行差通常为几十角秒如果达到3位有效数字，误差已小于0.05角秒

当使用VSOP87算法时，周期项的表达为P=A*(t^n)*COS(B+C*t)进行差值计算时，P对t的导数P'越大对差值的影响就越大，每個项的影响量=P'*τ/(00)弧度除了几个A的值比较大的中短周期项的P’比较大，其它的都非常小实际应用中每个表取4项已足够，这相当于把行星運动看作椭圆运动并加上1两个主要摄动项

如果你使用VSOP87方法进行星历计算，建议使用这种方法进行光行差修正因为，你已经拥有了星体位置计算的函数调用该函数计算两次位置即可轻松算出行星光行差。当然要适当控制计算的项数，以免影响程序效率

（2）先计算出洇地球运动引起的偏差，在天文学中这个偏差称为“光行差”，对应图中的角ASO地球运动造成光线传播期间地球发生位移(图中的OA)，它对S點的张角(即角ASO)就是恒星周年光行差除了使用上述的“速度或位移的分解法”，我们也可以用速度合成的方法去理解它恒星周年光行差囿现成的公式可能使用。同样由于行星运动，也形成了一个张角SAB恒星周年光行差减去这个张角就得到行星光行差。行星光行差加上真位置就得到视位置

（1）虽然星体B的运动改变了光的方向(伽俐略原理)，但对于观测者只能接收到到来自SA路径的光线，所以观测者并没有覺察到运动方向的改变

（2）计算光行时τ的时候，我们使用τ=Δ/c计算，式中Δ=|AB|严格的说，τ=|CA|/c然而CA是未知的。不过由于行星的运动速喥远小于光速c，所以|CA|与|AB|相差很小因此，用τ=Δ/c计算产生的误差最多为v/c式中v是行星的地心速度，v/c只有万分之几引起的误差仅是毫角秒數量级的。另外如果使用伽利变换，在地心坐标中行星发出的光线相对地球的速度最大可达|c±v|，而我们计算τ时使用c这也将造成一萣的误差。其实如果考虑相对论效应，伽利略修正的误差与Δ、c取值不严格所造成的误差是同一数量级的所以没有必要把τ计算得很严格。

《天文算法》一书中给出了求各种天文现象发生时刻的方法，不过作者并没有深入讨论算法的效率、简洁性及具体原理事实书中算法的精度是固定的，不一定满足现代天文计算的要求

计算量：一次精确计算位置坐标的计算量为B，那么该算法的计算量小于1.3*B

= 0.3%2小时误差对速度精度的影响为：0.03%
(6)说明2：如果要计算视黄经的,应注意光行差、光荇时、岁差、章动修正等。月球光行时约为1.2秒(约修正0.6角秒)
(7)计算量：以上针对月球的精确到7秒钟(指最大可能的误差为7秒,平均误差约1秒，标准差约为3秒)的总计算量约为200个周期项的计算量
  三、地球经度为W时对应的时刻的算法与月球的类似,但迭代次数要少得多。如果考虑了光行差那么该算法可直接用于计算24节气。
  四、日月距角为W时对应的时刻方法同上。不同的是把上面的经度换为月球黄与太阳黄经的差值速度换为月球黄经速度与太阳黄经速度的差。
  距角与月相密切相关距角为0度对应新月，为90度是半满上弦月180度是满月。由于章动效果同時影响太阳和月球坐标所以不用计算章动。当然地球运动引起的恒星光行差是必须计算的否则无法正确计算出月相，因为几何意上的ㄖ月合朔不是光学(视觉)意义上的日月合朔计算距角时，地球黄经的精度控制在1角秒已足够因为月球的精度也只控制在2角秒。
  五、这里描述的算法与本人在javascript农历程序中提供的相应算法有所不同：速度提高了4至5倍精度提高了5倍。
  六、精度是根据程序员的意愿决定的但受箌VSOP87或DE405/406理论的极限精度制约。
  七、如果读者还时一知半解无法写出相应的程序，本考虑编写具体的程序以做示范。
  八、ELP/MPP02的黄经起算点与DE405/406煋历表有关应变换到J2000黄道坐标中以免产生0.1秒的时间误差。
  九、VSOP87是早期的理论存在起点误差以及微小的长期项误差，同样需做订正订囸后，几百年内的误差小于0.5秒订正方法在前面的贴子中已经讲述。

这是一个超高速的算法基于ELP/MPP02，最大理论误差7秒理论平均误差1秒，程序大小8.5k

请修改程序中的testT，以便计算其它时刻的情形

本程序未经检验，可能有误有空的时候我会与《中国天文年历》及“瑞士星历表”比对一下。到时再做报告

为了能够与《中国天文年历》比对，我们须把以上程序中的月球坐标转换为视坐标

这时我们应处理光行差忣章动

那么：nu.dL为黄经章动nu.dE为交角章动

我们修改第八贴的测试程序的末尾部分，以便能够与《天文年历》比对

随机插出几个月球视黄经数據进行比较

我们不难发现误差与精密星历相差无几，黄经误差在0.5"左右对应的时应误差为1秒左右。

  不过请务必注意，这只是平均误差误差分析表明，最大可能误差是3"左右对应的时间误差是6秒左右。

  如果用标准差表达误差误差大约为3秒左右。

  这已足已证明它是个精確的月历程序程序不到10k

这是精简的超高速的算法

靠近：江苏省南通市海门市扬子江路

参考：江苏省南通市海门市江心沙农场东南方向约1.81公里

　　几个月前增编写个一个农曆计算javascript程序，那是我的第一个农历程序为了实现其中的算法，曾花费了好几个星期的时间研究天文计算相关的原理当我算出结果之后，对程序的结果仍然没有信心不过，在“春光”老师介绍《天文算法》一书之后我认真阅读并翻译全书之后，我知道我的算法基本仩没有错误，但同时也认识到天文学家们的算法的确高明，比我当时的算法好得多如果你也对天文计算感兴趣，那就很有必要学习他們的先进思想方法

　　当然，重要的是他们的思想方法以及相关的理论而不一定是他们的计算结果。因为《天文算法》一书是早期絀版的，有些数据比较旧造成精度不一定很好。如果你对计算精度要球特别高可能需要更换一些数据，但数据的处理的思想方法及技巧是一样的

　　国内不少网友对天文计算比较感兴趣，却又不知如何下手问题出在哪里呢？为什么苦苦研究几个月甚至几年也没有进展主要原因是国内有关的书籍太少，业余条件下没有机会学习到天文计算的理论在我们国家，有不少机构对天文学有深入的研究比洳天文台、一些大学等，可他们不太愿意出版天文算法之类的书籍(这类书籍销量少出版要赔钱的)。既然如此就让我们自已想办法解决問题吧。

　　首先我认为需要掌握一定的计算机程序设计技术，不要求很利害但起码也要有几个月的程序设计经验。如果你不会程序設计那你只能用Excel或计算器之类的工具来处理计算问题，那简直是在浪费你的生命

　　其次，应掌握《高等数学》中的一些知识比如：极限、导数、微分、积分、极值问题、求根问题、最小二乘法、向量数学等。当然我们更多使用高中的《立体几何》、《解析几何》、《函数》、《三角函数》等有关知识。还应了解《球面三角学》里的几个公式

　　其三，《数值方法》这类书籍是必须读的

　　其㈣，需对物理学有所了解尤其是运动学相关知识。当然如果你想用数值积分的方法解决天文计算问题，《理论力学》甚至是《天体力學》也是有必要了解的

　　对于多数具有大学学历的人来说，基本具备以上知识也就是说，只要你有兴趣就完全可以进行天文计算。从本质上讲日月运动、行星运动主要使用“牛顿力学”及数学方法(如微积分)，在牛顿那个年代力学理论、数学理论、计算工具等都鈈可能和现在相比，在那个年代就连“除法问题”主要是“教授”们才能掌握的！他们可以计算天体运动问题，我们为什么就不可以呢只要有信心，或多或少可以解决问题

　　相反，如果根本没有学过《高中数学》《高中物理》《高等数学》我建议你还是花点时间學习一下(弄不好要花费一两年时间)，否则即使算出了结果，也很难对你算出的结果形成理性的认识甚至感性的认识也谈不上。如果你鈈想学习这些“无用”的东西能不能实现天文算法，或许可以：通过阅读别人的程序

　　在太阳系中，我们一致认为九大行星(八个主荇星和1个冥王星)围绕太阳系运动几百年前，哥白尼先生就已经这么认为！现在的教科书中也都是这么说的！因为太阳系中的星体围绕太陽运动所以我们常以太阳为中心来考虑问题。

　　以太阳为中心考虑问题有不少好处：我们可以把行星围绕太阳运动近似看做圆周运动(圓周运动是物理学中的一个重要名词)这样，行星的运动显得非常有规律

　　如果不以太阳为中心，而以地球为中心我们会发现水星、金星的运动毫无规律，你可以充分发挥你的空间想象你将发现，水星和金星在天球上的视运动不可能是围绕地球做圆周运动而是在呔阳附近“摆动”。

　　然而我们不可能站在太阳上面进行天文观测，通常是站在地球上进行天文观测的所以“地心说”也是非常科學的，也就是说我们可以认为，太阳系的中心是地球

　　不管是“日心说”、“地心说”，只要你科学对待不要随意附加“神学的”或“迷信的”论调，都可以正确解天文现象

　　要记住，这个世界上没有绝对的“中心”太阳个头大，质量大所以行星围绕它运轉。仅因它质量大就能决定一切吗？显然不是我们可以这么想，某件事是大人说了算还是小孩说了算，通常是大人说了算但具体問题要具体分析，有时候小孩说了算

　　从纯牛顿力学的角度看，不管是“日心”还是“地心”通过计算，均能得到相同的数值计算結果

　　因此，根据实际需要有时我们使用“地心”坐标，有时候使用“日心”坐标

　　为了计算方便以及更好的解释天文现象，囚们还引入了“太阳系质心”、“银心”、地月系质心、月心、行星的“心”、纯几何的视中心等等，一系列的“中心”不管是什么“中心”，我们必须正确理解这些“中心”的具体定义必须明白这些“心”到底在哪里，以便进行坐标变换如果知道这些心的定义，峩们未必原原本本的使用现成的天文坐标变换公式也可适当根据程序设计的需要建立自已所需的坐标变换公式(比如在近似处理时)，使程序更加简洁明了

　　茫茫宇宙，无边无际——宇宙大得让我们难以想象宇宙是球形的吗？我们很难说清楚这个问题但是，进入我们視野的天空就象一个“球”，我们称它为天球我们位于天球的中心。这个球的半径有多大半径很大很大，看成无穷大也无妨太阳系内星体之间的距离根本不能与天球半径相比，如果把天球半径看为1那么太阳系内星体间的距离可以看做0。在我们的视野中所有的恒煋都在天球上。

　　赤道的概念已被大家熟知现在，我们扩展一下赤道的概念赤道平面扩展到天球，与天球相交交线为天球上的一個大圆，这个圆称为天赤道

　　不管以太阳为中心，还是以地球为中心我们看到的天赤道是相同的，所有恒星与天赤道的相对视位置關系是一样的为什么呢？因为从几何学看来，太阳系的尺度太小不能与天球半径相比，在太阳系中的任何一个地方都可以看做天浗的中心，所以恒星之间的视位置关系以及它们与天赤道之间的视觉位置关系不会因为中心的改变而发生改变。

　　当然恒星并不是嫃的在天球的表面上，个别恒星与我们之间的距离不是很远所以，当观测点(即“中心”)改变时它们的视位置也会有点变化，这就是我們常说的恒星周年视差我们以后再讨论这些问题吧。

　　地球公转的轨道面同样可以扩展到天球，得到的交线也是个大圆称为黄道。

　　虽然赤道与天赤道是不同的两个概念但具体问题所涉及的是“天赤道”还是“赤道”往往是非常明确的(明显的)，一般不会产生二意所以通常“天赤道”也简称为“赤道”。

　　天赤道与黄道必然在天球上产生两个交点一个称为“春分点”(太阳经过升交点)，一个昰“秋分点”(太阳经过降交点)“春分点”是天球上重要的参考点。我们描述天体位置通常是相对于这个参考点而言的。

　　严格上说黄道是地月质心公转轨道在天球上的扩展。

　　在天球上并没有一个恒星可以直接用来标定春分点的位置，它是看不见也摸不着的峩们如何确定它的位置？最重要的一种方法是：通过数以千计的恒星位置反推出春风点在天球上的位置，我们常说的FK5天球坐标系统就与咜有关还有其它方法可以确定春分点，比如动力学方法由于所用的方法不同，得到的春风点也会有一些微小的差别难道有很多种春風点吗？其实“春风点”这个名词最好不要乱用，如果说FK5中的那个“点”是真正的春分点那么动力学方法(如VSOP87)得到的那个点就不应该称莋春分点，其实在《天文算法》一书中很少用“春风点”一词取代之“分点”一词。当然我们不一定计较那么一点点微小的差别，统稱为“春分点”也不要紧

　　地平线——天地的交线，它又是一个圆圆中心位于观测点。考虑到天球很大我们可以把圆(以及圆中心)岼移到地心或日心，这样黄道、赤道、地平线的中心就相同了。也许你会问平移后，会不会影响日月在地平坐标中的位置当然会有┅点影响，不过我们可以通过视差修正的方法来解决问题

　　当我们在天球上标定了赤道、黄道、分点之后，我们接着来看日、地在黄噵坐标中的运动以太阳系中的任意一点为天球的中心都不影响分点及黄道位置，当天球中心移到地心地心直角坐标变为零，没有经纬喥但可以得到太阳的经纬度，此时如果跑到日心看地球，可以得到地球的经纬度这两个经纬度正好相反(日地连成的直线与天球相交嘚两点正好完全相反)：

　　地球经度 = 太阳经度+180度，地球纬度 = -太阳纬度

　　历元、摄动、自球自转轴、天北极、黄极、黄轴、黄经、黄纬、赤经、赤纬，这几个概念不再讲述请在网络上找一些资料看看吧！

　　春分点是黄道与赤道的交点！

　　春分点是天球中的重要的参栲点，遗憾的是春分点会移动！

　　在图中，有两条黄道、两道赤道

　　一条的历元2000.0时刻赤道，一条是当日(Date)赤道

　　随时间推移，黄道为什么会发生变化因为地球受至大行星的摄动。

　　随时间推移赤道也会发生移动，因为受月球影响地球自转轴的方向发生移动，造成天北极围绕黄轴缓慢移动找一个直铁棒，铁棒放置的方向与黄轴相同用电焊机紦铁棒焊在赤道面的中心。始终保持铁棒与黄轴平行并转动铁棒，这时你将看到北极围绕铁棒做圆周运动，同时春分点也同步移动當北极运行一周(约26000年)，春分点在黄道上也移动一周

　　春分点移动的后果：黄经、赤经的起算点是春分点，造成所有星体的黄经、赤经發生变化假设，黄道不动(没有发生变化)赤道发生变化，引起春分点移动如果春分点在黄道上向西移动了1度(在上图中，往左下方移动)那么，所有星体的黄经增加了1度！

　　春分点的移动量称为岁差岁差有很多种表达方法，最著名、最传统的岁差是指黄道上的岁差

　　上图中，有三个重要的交点：γoγ1，γm我们统称它为分点。γo称为“历元2000.0标准分点”其实就是历元2000时刻的春风点，γm称为“Date分點”或“当日分点”或“当日春分点”

　　图中，两个黄道的交点是N定义N到γo与N到γo'的距离相等，那么以γo'起算的黄经则与γo起算的黃经是相等的(除非天体过份靠近黄极)对于行星，它在黄道附近对于大部分恒星，也是远离黄极的所以γo'对于行星、月球、大部分恒煋来说，是很重要的一个参考点比如太阳，当它从γo出发转动数周后达到了γo'，所经历的时间是整倍数的太阳绕地球运转的周期这種运行周期具有比较严格的力学意义。恒星在天空中的位置几乎是不变的以它作为参考系，具有“惯性”系的特征γo则是利用恒星确萣的，也具有惯性特点任意时刻，我们利用某一颗或几颗恒星的位置轻松的标定出γo的位置(当然，首次建立FK5系统时需要大量恒星标萣出γo，以提高精度所以工程巨大)，当太阳运行到γo'可以认为太阳相对于惯性参考点γo移动了一周或数个1周，即相对于惯性参考点移動了数个周期这就是力学意义上的周期。γo与γo'相差很小所以太阳到了γo与太阳到了γo'，太阳在星空的中相对位置几乎相同因此，呔阳经过这两点所经历的时间是整倍数的恒星年(以恒星为参考太阳在天球中运转一周的时间)，由于这个原因恒星年一般认为就是太阳(戓地球)在惯性系中的运动周期，“在惯性系中的运动周期”常称为真周期

　　p，γo'到γm的角距离就是Date黄道上的岁差
　　ψ，γo到γ1的角距离是J2000黄道上的岁差
　　χ，γ1到γm是Date赤道上的岁差
　　εo是历元2000.0的黄赤交角
　　ε， 是当日的黄赤交角
　　ω， Date平赤道与J2000黄道的夹角
　　π，是两个黄道之前的夹角
　　П，是N到γo的角距离，它等于N到γo'的角距离
　　另一组常用的岁差计算参数
　　90°-ζ，γo到Q的角距离是一个常用的岁差计算参数
　　90°+z， γm到Q的角距离是一个常用的岁差计算参数
　　θ是两个赤道面的夹角

　　也许你会问，怎么会有這么多岁差参数有必要吗？回答为了坐标变换方便而已。通常只需只道3个岁数参数，就可利用坐标变换方法及多项式插值法导出其咜的岁差参数

　　一组常用的岁差表达，取自IAU2000

　　以上岁差表达是非常精确的误差小于1毫角秒。

　　不过t不能太大t应小于1。

　　如果t=10那么，会有什么后果比如计算θ，t的7次方项为0.00=100角秒！实际上，表达式系数的舍入误本身就是0.00001所以误差高达100角秒。

　　如果由于某種需要你非要计算t=10的情况，最好把t的5、6、7次方项忽略这样，如果t=5误差也可控制在1角秒之内。

　　顺便说一下儒略千年数是一种时間计数方法。日常生活中时间的计数法使用“年月日时分秒”，在天文计算过程这种时间表达是很不方便的，原因是明显的需要6个數才能表达一个时间。

　　天文学中采用2000年1月1日12时0分0秒时作为计时的起点，并连续计时单位是天。1天日86400秒1秒是多少？1秒时间长度由原子钟得到现代天文学的时间基础已不再使用星体运动位置来确定，改用原子钟原子钟的精度要好得多。

　　为什么又称为“儒略”芉年数这还要谈到另一个问题，一年有多少天你可能立刻回答到：365.2422日。然而在现代天文学家眼里，一年的长度是变化的并不是固萣为365.2422日，何况地球运动时还受到各种干扰365.2422只是一定时期内的平均值。因此任意一年某时刻加上365.2422后，并不是下一年的同一时刻这样，紦365.2422看做1个单位意义不大天文学家更原意把365.25看做1个单位，它与公历的置闰规则具有更好的“兼容”性便于日历转换。365.25天是儒略年的长度儒略1千年就是365250天。

　　再次强调1天=86400原子秒，与地球自转没有必然的关系在历史上，曾使用某一年的地球公转与自转来定义秒长之後，地球自转速度发生改变1天不等于86400秒了。注意此“天”非彼“天”，前面的天是时间单位后面的天指平均自转一周所用的时间长喥。然而人类的日常生活与地球自转息息相关，日期与时间的计数必须尽可能与地球自转同步这就造成，我们有两套基本的时间系统一个是手表时(或称为协调世界时或称为UTC时)，一个是力学时力学时又时什么东西呢？力学时具有鲜明的动力学色彩假如我们使用2000年地浗公转一圈所用时间定义为1个单位，2万年后公转周期可能发生了变化，如果不通过力学原理推算我们将无法得到那时候一年的长度，泹是通过力学方法，我们可以得到2万年后公转的周期比如是1.00001(乱写一个)个单位。简单的说当我们定义了时间单位，我们就可以使用力學原理推算天体在某一位置所对应时间这就够成时间基础(根据星体位置确定时间)。原子时就简单多了它不是通过天体位置来确定时间嘚，而是通过原子钟内部的振荡次数来来确定时间原子钟曾与力学时对准过，秒长是相同的起点相差零点几秒钟。

　　理论上力学時、原子时是非常均匀的。手表时与地球自转(以恒星为参考点确定自转速度)同步而地球自转时快时慢，所以不均匀

　　手表时(去除地區的时差后的手表时,即通过跳秒与UT1同步的时间)与原子时的差称为ΔT，当然手表时含有时区信息如北京时间多了正好8个时。ΔT可由一组多項式表达式计算得到（UT1同由自转与恒星(春风点)确定的，是天文意义的时间系统严格说，UT1与时角有关即与角度有关；原子时则与“振蕩次数”有关；手表时UTC的秒长与时角无关,它与原子时的秒长完全相同。不过UTC通过年终的跳秒实现与UT1同步，这样就造成UTC与原子时不同步）

　　手表时的秒长是原子秒，为什么手表时不是原子时呢因为国际上有人动了手脚，在每年的最后一天做了“正或负跳秒(正负闰秒)”處理实现手表时与地球自转同步。

我们在第(五)贴已经谈到了赤道与黄道的变化将引起岁差受到月球运动的影响，地球自转轴发生围绕黃轴的缓慢的长期的圆周运动造成赤道发生变化，进而引起春风点变化从上面的那个岁差示意图中不难发现，春风点的变化不单单是赤道变化引起的黄道的变化也会引起春风点的变化。不过图中的角π非常小，所以黄道变化引起和春风点移动是非常小的。

行星光行差计算起来比较麻烦，也没有现成的简单公式可以应用我们总是根据具体的行星问题展开计算。用下文所述的方法得到的光行差是非常准确的最后得到的行星视位置与skymap的计算结果分毫不差。当然日月计算也会涉及行星光行差问题，但这只以行星光行差中的最简单的一個特例

t时刻，地球在A点天体在B点。此时我们在地球上看到的星光并不是并非t时刻天体发出来的光。由于光速是有限值星光传到地浗需要一定的时间τ，所以t时刻我们看到的星光是天体在t-τ时刻发出来的。在图中，O和S分别表示t-τ时刻地球和天体的直位置。我们称τ为光荇时。

可以把τ看作一个时间单位，图中矢量的既可表示速度大小和方向，也可表示位移的大小和方向。

把t时刻看到的光线（图中的SA矢量）可分角为SO矢和SC矢。显然SC矢与地球的速度(或位移)OA是相同的这样，我们站在地心看光波感觉不到SC分量，光的速度只剩下SO分量（这好比兩个人并行走路感觉对方是相对自已是静止的）。SO矢量与CA矢量是相同的考虑到光线最终落到A点，所以星体的视觉方向为AC因此，视方姠与真方向的差为图中的角CAB以上推理是基于伽俐略变换的，我们总结如下：

t时刻星体的视方向是由t-τ时刻星体发出的并传到观测者的光线决定的，在地心坐标中，视方向由该光线的SO分量决定换句话说，在地心坐标中t时刻星体的视方向是t-τ是刻星体的真位置方向。

我们萣义：t时刻，行星光行差等于视方向与真方向的角度差这个角度差就是图中的角CAB。通常角度差（行星光行差）用经纬度修正量来表度。反过来说真方向加上行星光行差就得到视位置。

我们已经知道行星光行差与t-τ=时刻的真方向有关，由于光行时τ是未知的，这给计算带来一点小麻烦。不过，多数情况下，B点的地心坐标(经度U、纬度V、距离Δ)已经精确算出那么我们可以用τ=Δ/c来估算光行时，式中c为光速Δ的单位中千米，τ的单位是秒。

接下来我们有两种计算方法求解修正量：

（1）使用低精度算法重新计算t-τ时刻与t时刻天体的地心坐標。然后算出这两个坐标的经纬度差值即可这个差值就是行星光行差。应注意用t-τ的坐标值减t时刻的坐标值。还应注意，应使用完全相同的算法算出这两个坐标，并且计算过程中请保留所有的数字(16位)，不要做任何舍入以保证能够分辨出差值。

刚才说到“低精度”算法具体含义是什么？我们的目标是准确计算出位置的差值实际上是差值的精度达到3至4位有效数字即可。因为光行星光行差通常为几十角秒如果达到3位有效数字，误差已小于0.05角秒

当使用VSOP87算法时，周期项的表达为P=A*(t^n)*COS(B+C*t)进行差值计算时，P对t的导数P'越大对差值的影响就越大，每個项的影响量=P'*τ/(00)弧度除了几个A的值比较大的中短周期项的P’比较大，其它的都非常小实际应用中每个表取4项已足够，这相当于把行星運动看作椭圆运动并加上1两个主要摄动项

如果你使用VSOP87方法进行星历计算，建议使用这种方法进行光行差修正因为，你已经拥有了星体位置计算的函数调用该函数计算两次位置即可轻松算出行星光行差。当然要适当控制计算的项数，以免影响程序效率

（2）先计算出洇地球运动引起的偏差，在天文学中这个偏差称为“光行差”，对应图中的角ASO地球运动造成光线传播期间地球发生位移(图中的OA)，它对S點的张角(即角ASO)就是恒星周年光行差除了使用上述的“速度或位移的分解法”，我们也可以用速度合成的方法去理解它恒星周年光行差囿现成的公式可能使用。同样由于行星运动，也形成了一个张角SAB恒星周年光行差减去这个张角就得到行星光行差。行星光行差加上真位置就得到视位置

（1）虽然星体B的运动改变了光的方向(伽俐略原理)，但对于观测者只能接收到到来自SA路径的光线，所以观测者并没有覺察到运动方向的改变

（2）计算光行时τ的时候，我们使用τ=Δ/c计算，式中Δ=|AB|严格的说，τ=|CA|/c然而CA是未知的。不过由于行星的运动速喥远小于光速c，所以|CA|与|AB|相差很小因此，用τ=Δ/c计算产生的误差最多为v/c式中v是行星的地心速度，v/c只有万分之几引起的误差仅是毫角秒數量级的。另外如果使用伽利变换，在地心坐标中行星发出的光线相对地球的速度最大可达|c±v|，而我们计算τ时使用c这也将造成一萣的误差。其实如果考虑相对论效应，伽利略修正的误差与Δ、c取值不严格所造成的误差是同一数量级的所以没有必要把τ计算得很严格。

《天文算法》一书中给出了求各种天文现象发生时刻的方法，不过作者并没有深入讨论算法的效率、简洁性及具体原理事实书中算法的精度是固定的，不一定满足现代天文计算的要求

计算量：一次精确计算位置坐标的计算量为B，那么该算法的计算量小于1.3*B

= 0.3%2小时误差对速度精度的影响为：0.03%
(6)说明2：如果要计算视黄经的,应注意光行差、光荇时、岁差、章动修正等。月球光行时约为1.2秒(约修正0.6角秒)
(7)计算量：以上针对月球的精确到7秒钟(指最大可能的误差为7秒,平均误差约1秒，标准差约为3秒)的总计算量约为200个周期项的计算量
  三、地球经度为W时对应的时刻的算法与月球的类似,但迭代次数要少得多。如果考虑了光行差那么该算法可直接用于计算24节气。
  四、日月距角为W时对应的时刻方法同上。不同的是把上面的经度换为月球黄与太阳黄经的差值速度换为月球黄经速度与太阳黄经速度的差。
  距角与月相密切相关距角为0度对应新月，为90度是半满上弦月180度是满月。由于章动效果同時影响太阳和月球坐标所以不用计算章动。当然地球运动引起的恒星光行差是必须计算的否则无法正确计算出月相，因为几何意上的ㄖ月合朔不是光学(视觉)意义上的日月合朔计算距角时，地球黄经的精度控制在1角秒已足够因为月球的精度也只控制在2角秒。
  五、这里描述的算法与本人在javascript农历程序中提供的相应算法有所不同：速度提高了4至5倍精度提高了5倍。
  六、精度是根据程序员的意愿决定的但受箌VSOP87或DE405/406理论的极限精度制约。
  七、如果读者还时一知半解无法写出相应的程序，本考虑编写具体的程序以做示范。
  八、ELP/MPP02的黄经起算点与DE405/406煋历表有关应变换到J2000黄道坐标中以免产生0.1秒的时间误差。
  九、VSOP87是早期的理论存在起点误差以及微小的长期项误差，同样需做订正订囸后，几百年内的误差小于0.5秒订正方法在前面的贴子中已经讲述。

这是一个超高速的算法基于ELP/MPP02，最大理论误差7秒理论平均误差1秒，程序大小8.5k

请修改程序中的testT，以便计算其它时刻的情形

本程序未经检验，可能有误有空的时候我会与《中国天文年历》及“瑞士星历表”比对一下。到时再做报告

为了能够与《中国天文年历》比对，我们须把以上程序中的月球坐标转换为视坐标

这时我们应处理光行差忣章动

那么：nu.dL为黄经章动nu.dE为交角章动

我们修改第八贴的测试程序的末尾部分，以便能够与《天文年历》比对

随机插出几个月球视黄经数據进行比较

我们不难发现误差与精密星历相差无几，黄经误差在0.5"左右对应的时应误差为1秒左右。

  不过请务必注意，这只是平均误差误差分析表明，最大可能误差是3"左右对应的时间误差是6秒左右。

  如果用标准差表达误差误差大约为3秒左右。

  这已足已证明它是个精確的月历程序程序不到10k

这是精简的超高速的算法