现代历法（整理）

最新推荐文章于 2020-02-17 14:26:08 发布

songgz

最新推荐文章于 2020-02-17 14:26:08 发布

阅读量8.5k

点赞数 13

分类专栏： OOP 文章标签：算法 function date c javascript 测试

本文链接：https://blog.csdn.net/songgz/article/details/2680144

版权

OOP 专栏收录该内容

6 篇文章 0 订阅

订阅专栏

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

　　一条的历元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
　　P:(     0.0,        41.99604, 19.39715, -0.22350, -0.01035, 0.00019, 0.0,     0.0)
　　Q:(    0.0,      -468.09550,   5.10421, 0.52228, -0.00569, -0.00014, 0.00001, 0.0)
　　π:(     0.0,       469.97560, -3.35050, -0.12370, 0.00030, 0.0,      0.0,     0.0)
　　II:(629543.988,   -8679.218,    15.342,    0.005,   -0.037,   -0.001,    0.0,     0.0)
　　p:(     0.0,     50287.92262, 111.24406, 0.07699, -0.23479, -0.00178, 0.00018, 0.00001)
　　θ:(     0.0,     20041.90936, -42.66980,-41.82364, -0.07291, -0.01127, 0.00036, 0.00009)
　　ζ:(     2.72767, 23060.80472, 30.23262, 18.01752, -0.05708, -0.03040, -0.00013, 0.0)
　　z :(    -2.72767, 23060.76070, 109.56768, 18.26676, -0.28276, -0.02486, -0.00005, 0.0)
　　ε:( 84381.40880, -468.36051, -0.01667, 1.99911, -0.00523, -0.00248, -0.00003, 0.0)
　　ω:( 84381.40880,    -0.26501,   5.12769, -7.72723, -0.00492, 0.03329, -0.00031,-0.00006)
　　ψ:(     0.0,     50384.78750,-107.19530, -1.14366, 1.32832, -0.00940, -0.00350, 0.00017)
　　χ:(     0.0,       105.57686,-238.13769, -1.21258, 1.70238, -0.00770, -0.00399, 0.00016)
　　上表是关于t的多项式系数，t等于J2000起算的儒略千年数，为儒略千年=365250天
第一列对应t的0次方，第二列对应t的1次方，第三列对应t的2次方，其余类推。
如π=0.0*t^0 + 469.97560*t^1 - 3.35050*t^2 + ...
结果的单位是角秒。

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

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

　　如果t=10，那么，会有什么后果？比如计算θ，t的7次方项为0.00001*10000000=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与原子时不同步。）

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

//原子时与UT1的差
//部分参考了美国国家航空航天局网站里的算法,www.nasa.gov
//2050以后的外推参考了skymap10.5,在skymap的星历表中分段找几个数据,再用最小二乘法拟合
double deltatT(double y){ //输入年份(可带小数点),如2000.78
static double DTxs[15][12]={//{范围下,范围上,偏移,时间单位,多项各系数}
   {-9999,-500,1820,100,-20,      0,        32},
   {-500, 500, 0,   100, 10583.6,-1014.41, 33.78311, -5.952053, -0.1798452,   0.022174192, 0.0090316521},
   {500, 1600,1000,100, 1574.2, -556.01,   71.23472, 0.319781, -0.8503463, -0.005050998, 0.0083572073},
   {1600, 1700,1600,1,   120,    -0.9808, -0.015320, 1.0/7129},
   {1700, 1800,1700,1,   8.83,    0.1603, -0.0059285, 0.00013336,-1.0/1174000},
   {1800, 1860,1800,1,   13.72, -0.332447, 0.0068612, 0.0041116, -0.00037436, 1.21272e-5, -1.699e-7, 8.75e-10},
   {1860, 1900,1860,1,   7.62,    0.5737, -0.2517540, 0.01680668,-0.0004473624,1.0/233174},
   {1900, 1920,1900,1, -2.79,    1.494119,-0.0598939, 0.0061966, -0.000197},
   {1920, 1941,1920,1,   21.20,   0.84493, -0.0761000, 0.0020936},
   {1941, 1961,1950,1,   29.07,   0.407,   -1.0/233,   1.0/2547},
   {1961, 1986,1975,1,   45.45,   1.067,   -1.0/260, -1.0/718},
   {1986, 2005,2000,1,   63.86,   0.3345, -0.060374, 0.0017275, 0.000651814, 0.00002373599},
   {2005, 2150,2000,100, 64.723,-11.050,    232.182, -86.08},
   {2150, 9999,2000,100, 59.033, 103.4335, 28.980314, 0.000276}
};
//计算TD-UT,力学时与世界时之差
double *p,dt=0,u;
int i,j;
for(i=0;i<15;i++){
   p=DTxs[i];
   if(y<p[0]||y>p[1]) continue;
   y=(y-p[2])/p[3];
   for(j=4,u=1;j<12;j++) dt+=u*p[j],u*=y;
   return dt;
}
return 0;
}

我们在第(五)贴已经谈到了赤道与黄道的变化将引起岁差。受到月球运动的影响，地球自转轴发生围绕黄轴的缓慢的长期的圆周运动，造成赤道发生变化，进而引起春风点变化。从上面的那个岁差示意图中不难发现，春风点的变化不单单是赤道变化引起的，黄道的变化也会引起春风点的变化。不过，图中的角π非常小，所以黄道变化引起和春风点移动是非常小的。
从岁差示道意图中的几何关系：
    弧Nγ1 = 弧Nγm + χcos(ω)，注意，这是近似计算，但精度仍非常高
两边同时减去П得：
    Nγ1 - П = (Nγm - П)+ χcos(ω)，即 ψ = p+χcos(ω)
其实ψ就是：黄道不动，赤道变化引起的黄经岁差，我们称之是日月黄经岁差。
其实χ就是：Date赤道上的行星岁差。
如果把ψ = p+χcos(ω)改写为p = ψ-χcos(ω)
那么 -χcos(ω) 就是行星黄经岁差。
同样，从几何关系，还可轻松得到行星黄纬岁差：χsin(ω)
----------------
χ的岁差速度是：106"/千年，加速度为-238"/平方千年
χ岁差造成春风点东移：0".106*cos(23.5°)/年
日月岁差(黄经)造成春风点西移：50.384.78750/年
黄经总岁差速度约为二者之差：50".29/年
----------------
[章动]：
受到日月运动的影响，地球自转轴发生振动，这种振动称为章动。“章”字，从何讲起？地球自转的振周期主要与月球轨道升交点运动有关，升交点运动周期为18.6年，古人称之为一章。
在第(五)贴的图中，已给出受到章动后的赤道位置。图中的Δψ就是黄经章动值。
章动的后果：春风点在黄道上移动了Δψ(称为黄经章动)，黄赤交角变化了Δε(称为交角章动)
考虑章动以后的赤道称为真赤道，否则为平赤道。
黄道的振动是非常微弱，所以一般不再区分“真”与“平”
要得到真春点起算的黄经，必须加上Δψ，要继续旋转到真赤道坐标中，应使用真黄赤交角(即加上Δε)
-----------------
Δψ与Δε的计算
前面已发了一个关于IAU1980与IAU2000B章动序列比较的贴子，里面含有相关的所有计算，请阅读程序。
----------------
精度控制
完整的计算，计算量偏大，程序冗长。如果你只是为了计算农历，根本不需要那么高的精度。
为什么呢？
因为，现代的紫金山中国农历是基于北京手表时的，这种时间与地球自转有关，考虑到地球地转的不规律性，无法对几百年后的北京时间进行测算，误差0.5分钟是正常的，如果对5千年后的北京时间预测计算，误差可达几小时。0.5分钟的误差意味差什么呢？意味着太阳位置误差可达1".2，月亮位置误差可达15"，因此没有必要把Δψ算得很准，0".1已经足够。
以下javascript程序给出几千年内精度为0".05的黄经章动及交角章动。
<script language=javascript>
/*********************
对IAU1980章动计算的简化处理
许剑伟 2008年5月2日
*********************/
t=1; //J2000.0起算的儒略世纪数
var nuB=new Array(//章动表
2.1824,   -33.75705, 36e-6,-171996,-1742,92025, 89,
3.5069, 1256.66393, 11e-6, -13187, -16, 5736,-31,
1.3375, 16799.41822,-51e-6, -2274,   -2, 977, -5,
4.3649,   -67.51409, 72e-6,   2062,    2, -895, 5,
0.0431, -628.30196, 3e-6, -1426,   34,   54, -1,
2.3556, 8328.69143,155e-6,    712,    1,   -7, 0,
3.4638, 1884.96589, 8e-6,   -517,   12, 224, -6,
5.4383, 16833.17527,-87e-6,   -386,   -4, 200, 0,
3.6931, 25128.10965,103e-6,   -301,    0, 129, -1,
3.5501,   628.36198, 13e-6,    217,   -5, -95, 3
);
var dL3=0,dE3=0;
for(i=0;i<nuB.length;i+=7){
c=nuB[i] +nuB[i+1]*t +nuB[i+2]*t*t;
dL3+=(nuB[i+3]+nuB[i+4]*t/10)*Math.sin(c); //黄经章动
dE3+=(nuB[i+5]+nuB[i+6]*t/10)*Math.cos(c); //交角章动
}
dL3/=10000; //黄经章动
dE3/=10000; //交角章动
document.write("[精简的章动计算,精度0.05角秒,单位角秒]<br>");
document.write("黄经章动：" + dL3 + " 交角章动：" + dE3 + "<br>");

</script>

---------------------
[问题]：程序中使用了“J2000.0起算的儒略世纪数”，你会计算它吗？
[例题]：比如，2002年1月1日12:00:00 TD的儒略世纪数T是多少？TD表示力学时。
[解]：
   J2000.0对应2000年1月1日12:00:00 TD
   由于2000年是闰年，含有366天，也就是说：
     2001年1月1日12:00:00加上366天后变为2001年1月1日12:00:00
   又因为2001年是平年，所以再加上365天就是2002年1月1日12:00:00
   所以，所求时刻相对于J2000.0相差366+365=731天
   那么，T = 731/36525 = 0.020013689

行星光行差计算起来比较麻烦，也没有现成的简单公式可以应用，我们总是根据具体的行星问题展开计算。用下文所述的方法得到的光行差是非常准确的，最后得到的行星视位置与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'*τ/(365250*86400)弧度。除了几个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取值不严格所造成的误差是同一数量级的，所以没有必要把τ计算得很严格。

通过以下改正,可以大幅度提高VSOP87的精度
1、以下对VSOP87日心球面地球序列拟合到DE406:
   利用VSOP87计算出黄经L后,计算黄经的修正量dL如下:
      dL = 0.2 -0.22*t -0.054*t2 + 0.1*cos(t)*t - 0.11*cos(1.5*t)*t -0.13*cos(3*t)
   式中t是J2000.0起算的儒略千年数
   最后 L = L+dL
   修正前的最大误差:
     J2000+-3000年1.00"
          +-2000年0.67"
          +-1000年0.39"
          +- 500年0.28"
          +- 250年0.17"
   修正后的最大误差:
     J2000+-3000年0.6"
          +-2000年0.25"
          +-1000年0.12"
          +- 500年0.02"
          +- 250年0.012"
   距离与黄纬不再修正,其最大误差为:
     距离：+-250年3E-8AU, +-1000年1E-7AU, +-3000年1E-7AU
     黄纬：+-250年0.01" , +-1000年0.024", +-3000年0.08"
2、以下对VSOP87的date球面地球序列拟合到DE406:
   dL = 0.2 +2.95*t -0.15*t2 - 0.02*t3
        - 0.01*sin(5*t) -0.13*cos(3*t)- 0.1*cos(1*t)*t
   最后L = L + dL
   修正后的最大误差:
     J2000+-3000年0.8"
          +-2000年0.4"
          +-1000年0.1"
          +- 500年0.02"
          +- 250年0.014"
黄经拟合时已考虑了最新的岁差速度订正: -2.9965"/儒略千年
黄纬与距离不用修正。
3、使用J2000.0动力学分点,因此J2000黄道坐标中的黄经与DE405/406相差一个估定值0.04"。
   注意1，使用J2000.0分点(非常接近ICRF的赤经起标点)得到的星体坐标才能与最新的《中国天文年历》相同。
   注意2,公元3000以后的拟合了“瑞士星历表”
4、许剑伟,地震后第12天

《天文算法》一书中给出了求各种天文现象发生时刻的方法，不过作者并没有深入讨论算法的效率、简洁性及具体原理。事实书中算法的精度是固定的，不一定满足现代天文计算的要求。
比如试图把农历的精度计算到30秒,使用《天文算法》一书中的算法将不能实现。本贴讨论高精度的、超高速的算法。
问题假设：设经度L(t) = Lo + v*t + f(t)，对于方程L(t)=W，求t的值。
式中速度v值较大,f(t)对t求导数也得到速度量dv,如果它相对于v很小,那么:
解法(针对VSOP87或ELP2000-82或MPP02)：
(1)第一步,初步确定t的值：
t1 = (W-Lo)/v
这样，L的偏差为 dL = L(t1) - W
那么t1的偏差估计为dt = dL/v
对t1修正 t2 = t1 + dt
(2)第二步,用第一步的方法重新执行一次,不过计算过程中可以考虑v使用v+dv代替。
(3)第三步,重复第二步计算,但v一定要使用v+dv代替
注意事项1:计算L(t)时不一定需要完全精度,第一、二步计算时只需计算几个周期项即可。第三步则需多算几个项。如果精度还不满意，则重复第三步
注意事项2:v的解析式可以使用VSOP87或ELP/MPP02的序列求导数得到,只需考虑速度量较大的几个项即可。
注意事项3:如果使用DE405/DE406星历表或“瑞士星历表”,可以考虑使用求变差的办法求v+dv。“变差”法，指求相近两个时间的坐标差值，除以时间差后得到速度。
说明1：以上方法本质上就是《天文算法》中的迭代逼近思想方法。
说明2：上面讲到的v的运动学本质是“平角速度”，用2*3.1415926535/v得到“平周期”，如地球运动的“平周期”是365.2422天
说明3：对v+dv的取值还需与截断误差合并分析，这当中涉及“加速度”的分析，将在以后的贴子中讨论。
说明4：以上所述的是通用方法，对于实际问题，有必要进行更仔细的分析。

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

算法的极限精度取决于VSOP或ELP/MPP02理论本身
《农历天文算法完全解——第一部分(分两部分讲述)》
算法本质：迭代技术
一、日月黄经速度的表达
(1)地球黄经速度的表达(Date黄道坐标中的速度,含岁差速度)
以下表达式中t是千年数,误差小于万分之3,数据由VSOP87地球黄经序列求导数得到:
v = 6283.32 //误差小于万分3
+210 *Math.sin(1.527+f)
+4.4 *Math.sin(1.48+f*2)
+12.9*Math.sin(5.82+f)*t
+0.55*Math.sin(4.21+f)*t*t;
式中f=6283.07585*t;
(2)月球黄经速度的表达(Date黄道坐标中的速度,不含岁差速度)
岁差速度远小于月球黄经速度,所以下速度计算忽略了岁差速度。
平速度是(单位：单位弧度/世纪,式中的t是世经数)：
v=8400; //精确值是8399.68473007193
周期项速度(对ELP/MPP02黄经序列求导数得到)
dv=914*sin( 0.785 +8328.69142*t); //误差5%
dv=914*sin( 0.785 +8328.69142*t) //误差不超过0.8%
+179*sin( 2.54 +15542.7543*t )
+160*sin( 0.19 +7214.0629*t )
+ 62*sin( 3.14 +16657.383 *t )
+ 34*sin( 4.8 +16866.932 *t )
+ 22*sin( 4.9 +23871.45 *t )
+ 12*sin( 2.6 +14914.45 *t );
dv=914*sin( 0.7848 +8328.691425*t+ 1.523*t2 ) //误差小于0.3%
+179*sin( 2.543 +15542.7543*t )
+160*sin( 0.1874 + 7214.0629*t )
+62 *sin( 3.14   +16657.3828*t )
+34 *sin( 4.827 +16866.9323*t )
+22 *sin( 4.9    +23871.4457*t )
+12 *sin( 2.59   +14914.4523*t )
+7 *sin(0.23    + 6585.7609*t )
+5 *sin(0.9     +25195.624 *t )
+5 *sin(2.32    - 7700.3895*t )
+5 *sin(3.88    + 8956.9934*t )
+5 *sin(0.49    + 7771.3771*t )
+4 *sin(5.5     +24986.074 *t )
+4 *sin(4.3     +22756.817 *t )
+2 *sin(1       +32200.137 *t )
+2 *sin(1.9     +14428.126 *t )
+2 *sin(0.4     +31085.509 *t )
+2 *sin(1.528   + 628.302 *t );
dv的值可根据精度的需要从上面表达式中取一种。
月球真速度为：
v=v-dv;
(3)月球平黄经与黄经的表达式(ELP/MPP02)
w1 = 3.81034409083
    +8399.68473007193*t
    -3.31895204255e-05*t2
    +3.11024944911e-08*t3
    -2.03282376489e-10*t4
L=w1+周期项+泊松项，由于泊松项较小仅需考虑几项就可以了，除非有高精度计算的需要。
周期项与泊松项的具体算法待叙。
二、求月球经度为W时所对应的时间
经度计算函数为 Lon(t,n); t为J2000.0起算的世纪数,t不超过+-30世纪;n为序列截取数量，计算步骤如下:
数据准备,表达式准备:
to=0;
v0 = 8400; //平速度
v1(t) = 8400 + 914*sin(0.785+8328.69142*t); //低精度表达,5%精度
v2(t) = 中精度速度表达(见上文),0.3%精度
v = v0
(1)a.余项计算：dW = W - L(to,0)
   b.时间计算：to = to+ dW/v
   c.速度计算：v=v1(to); //为下一次计算准备速度,速度误差5%左右
   以上两行相当于解方程 W=3.81034 +8400*t,解为to = (W - 3.81034)/8400
   to是月球黄经为W是的“平”时间
   to与真时间的差值不超过18小时
   以上计算已经考虑了黄经的常数项(3.81034)和速度项(8400)。因为to的超不超过30，所以平黄经的高次项(加速度项)很小，不超过2度，黄经周期项不超过8度，所以黄经的的高次项及周期项总和不超过10度。因此如果把to再次回代到完整的序列中计算，上述方程最多只有10度(约0.17弧度)的误差：
W=3.81034 +8400*t+0.17，得to = (W -3.81034 -0.17)/8400，与原来的to比较不难发现，最大误差就是月球运动0.17(即10度)所用的时间,约18小时。这18个小时的误差主要是截断引起的。
(2)a.余项计算：dW = W - Lon(to,3);
   b.时间计算：to=to+dW/v;
   c.速度计算：v=v2(to); //为下一次计算准备速度,精度为0.3%，如果考虑由于时间不准确引起的额外误差，精度仍可达到0.35%(取0.4%)
   在步骤(1)中已经说到,误差为18小时,我们称之为余项(或称之为截断误差),而本次计算已经截取了3个主要周期项，可修正这18小时的余项。由于速度的精度只有5%，所以残留误差为18*5%=0.9小时。不过，实际误差不止这些，因为本次计算也存在截断误差。
   截断部分误差估计: 把被截断的最大10项的振幅取和作为最大误差估计,以后的项不用算(因为存在正负相消的作用),这10项取和为3000角秒=50角分，对应的最大时间误差为1.7小时。与刚才讨到的0.9小时综合考虑，不妨把误差估计为2小时。
   为什么不把误差估计为0.9+1.7=2.6小时？因为0.9与1.7的误差估计值均指最大误差，出现最大误差的可能性已经很小了，二者同时出现最大值的可能性几乎为0。
   当然，平均误差仅为最大误差的六分之一左右。
   第(1)步与第(2)步计算无需考虑光行时、章动的修正计算，因为误差远超过这些修正量。
(3)a.余项计算：dW = W - Lon(to,21);
   b.时间计算：to = to + dW/v;
   c.无需再计算更精确的速度v
   周期项取21位的理论精度是4',对应的时间误差为8分钟。
   速度不精确引起的误差：2小时*0.004=0.008小时=0.48分钟。
   总误差不妨按8分钟估计。
(4)a.余项计算：dW = W - Lon(to,200);
   b.时间计算：to = to+dW/v;
计算200项的理论最大误差约2角秒,对应的时间误差为4秒钟
速度v不精确引起的误差：8分*0.004=0.3秒
总误差不妨取4秒钟。
to就是我们最终所需要计算的时间！
(5)如果相要再提高精度,继续计算：
   dW = W - Lon(to,全部);
   to = to+dW/v;
“全部”计算的精度为5毫角秒，相当于0.01秒
速度v不精确引起的误差：5秒*0.004=0.02秒
总误差不妨取0.03秒钟。
ELP/MPP02在J2000前后几十年内与DE406拟合得很好，差异为0.005角秒以内。
(6)说明：由于都是在邻近的时间内计算，可以认为速度不变(或者说加速度为0)。因为对ELP/MPP02求二次导数，我们发现加速度非常小。为什么会这样呢？其实，周期项的周期最小值为10天左右(个别为5天左右,但振幅非常小)，也就是说在10天范围内速度大小才小幅波动一次，在几小时或几分钟内，速度几乎不变。速度大小波动幅度最大的周期项的周期是27.5天左右，它对速度的影响可达10%以上，即便是这种周期项，10小时对速度的影响最多只有：10%*10小时/27.5天 = 0.15%，2小时误差对速度精度的影响为：0.03%
(7)说明2：如果要计算视黄经的,应注意光行差、光行时、岁差、章动修正等。月球光行时约为1.2秒(约修正0.6角秒)。
(8)计算量：以上针对月球的精确到5秒钟(指最大可能的误差为5秒,平均误差不到1秒，标准差约为3秒)的总计算量约为250个周期项的计算量。
三、地球经度为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秒。订正方法在前面的贴子中已经讲述。
十、使用真正的DE405/DE406方法的计算时间，方法同上，只是在计算Lon()时使用DE405/406方法
改进《农历天文算法完全解——第一部分》
算法本质：迭代技术
一、日月黄经速度的表达
(1)地球黄经速度的表达(Date黄道坐标中的速度,含岁差速度)
以下表达式中t是千年数,误差小于万分之3,数据由VSOP87地球黄经序列求导数得到:
v = 6283.32
    +210 *sin(1.5274 + 6283.07585*t) //误差小于万分3
    +4.4 *sin(1.4845 + 12566.15170*t)
    +12.9*sin( 2.6782 + 6283.07585*t)*t
    +0.55*sin( 1.0721 + 6283.07585*t)*t2;
(2)月球黄经速度的表达(Date黄道坐标中的速度,不含岁差速度)
岁差速度远小于月球黄经速度,所以下速度计算忽略了岁差速度。
平速度是(单位：单位弧度/世纪,式中的t是世经数)：v=8399.68473007193
考虑岁差速度后,应把式中的8399.68473007193换为：8399.68473007193弧度/世纪+5028.79角秒/世纪=8399.70911033384弧度/世纪
周期项速度(对ELP/MPP02黄经序列求导数得到)
dv=914*sin( 0.785 +8328.69142*t); //误差5%
dv=914*sin( 0.785 +8328.69142*t) //误差不超过0.8%
+179*sin( 2.54 +15542.7543*t )
+160*sin( 0.19 +7214.0629*t )
+ 62*sin( 3.14 +16657.383 *t )
+ 34*sin( 4.8 +16866.932 *t )
+ 22*sin( 4.9 +23871.45 *t )
+ 12*sin( 2.6 +14914.45 *t );
dv=914*sin( 0.7848 +8328.691425*t+ 1.523*t2 /10000) //误差小于0.3%
+179*sin( 2.543 +15542.7543*t )
+160*sin( 0.1874 + 7214.0629*t )
+62 *sin( 3.14   +16657.3828*t )
+34 *sin( 4.827 +16866.9323*t )
+22 *sin( 4.9    +23871.4457*t )
+12 *sin( 2.59   +14914.4523*t )
+7 *sin(0.23    + 6585.7609*t )
+5 *sin(0.9     +25195.624 *t )
+5 *sin(2.32    - 7700.3895*t )
+5 *sin(3.88    + 8956.9934*t )
+5 *sin(0.49    + 7771.3771*t )
+4 *sin(5.5     +24986.074 *t )
+4 *sin(4.3     +22756.817 *t )
+2 *sin(1       +32200.137 *t )
+2 *sin(1.9     +14428.126 *t )
+2 *sin(0.4     +31085.509 *t )
+2 *sin(1.528   + 628.302 *t );
dv的值可根据精度的需要从上面表达式中取一种。
月球真速度为：v=v-dv;
(3)月球平黄经与黄经的表达式(ELP/MPP02)
w1 = 3.81034409083
    +8399.68473007193*t
    -3.31895204255e-05*t2
    +3.11024944911e-08*t3
    -2.03282376489e-10*t4
注意:w1中不含岁差速度
(4)黄经表达式(相对于Date平分点): L=w1+周期项+泊松项，由于泊松项较小仅需考虑几项就可以了，除非有高精度计算的需要。
(5)周期项与泊松项的具体算法待叙。

二、求月球经度为W时所对应的时间
经度计算函数为 Lon(t,n); t为J2000.0起算的世纪数,t不超过+-30世纪;n为序列截取数量，计算步骤如下:
数据准备,表达式准备:
to=0;
v0 = 8399.70911033384; //平速度(含岁差速度)
v1(t) = 速度表达,0.3%精度(见上文)
(1)a.余项计算：dW = W - L(to,0)
   b.时间计算：to = to+ dW/v0
   c.速度计算：v=v1(to); //为下一次计算准备速度,速度误差0.3%左右
   以上两行相当于解方程 W=3.81034 +8399.70911033384*t,解为to = (W - 3.81034)/8399.70911033384
   ▲to是月球黄经为W是的“平”时间，to与真时间的差值不超过18小时
   ▲以上计算已经考虑了黄经的常数项(3.81034)和速度项(8399.70911033384)。因为to不超过30世纪，所以平黄经的高次项(加速度项)很小，不超过2度。而黄经周期项不超过8度，所以黄经的的高次项及周期项总和不超过10度。因此如果把to再次回代到完整的序列中计算，上述方程最多只有10度(约0.17弧度)的误差：
       W=3.81034 +8400*t+0.17，得to = (W -3.81034 -0.17)/8399.70911033384
   与原来的to比较不难发现，最大误差就是月球运动0.17(即10度)所用的时间,约18小时。这18个小时的误差就是截断高次项及周期项引起的。
   ▲其实最后得到的v的最大误差大于0.3%，下文还会讨论到,由于to误差18小时,也将造成v产生0.3%的误差。v的总误差估计为0.6%
   ▲请务必注意，“平速度”的精度要求比“速度”的精度要求要高得多。原因如下：
   t=(W-3.81034)/(v0+dv)=to+dt，式中dv是平速度的误差；利用二项式定理展开后得：
   t=[(W-3.81034)/vo] * (1-*dv/Vo) = to*(1-dv/vo)
   当W很大时(即to很大时)：dt = t - to = to*dv/vo，式中dt就是dv引起的时间误差
   设 to=20世纪，dv=1弧度/世纪，我们得到：dt=0.0024世纪=2087小时！而v的相对误差只有dv/v=1/8400=0.012%
   所以“平速度”的精度要求非常高。
   ▲不含岁差的平速度v=8399.68473007193，当然这指的是角速度，那么月亮平周期就是：
     T=2*3.1415926535897933/v=0.00074802632587923106世纪=27.32166155273891天(1世经=36525天)，这就是恒星月
    还应注意，月球黄经还有“加速度”项
(2)a.余项计算：dW = W - Lon(to,21);
   b.时间计算：to=to+dW/v;
   在步骤(1)中已经说到,误差为18小时,我们称之为余项(或称之为截断误差),而本次计算已经截取了21个主要周期项，可修正这18小时的余项。由于速度的精度只有0.6%，所以残留误差为18*0.6%小时=7分钟。
   不过，实际误差不止这些，因为本次计算也存在截断误差。
   ▲截断部分误差估计: 把被截断的最大10项的振幅取和作为最大误差估计,以后的项不用算(因为存在正负相消的作用),这10项取和为130角秒，对应的最大时间误差为5分钟。与刚才讨到的7分钟综合考虑，不妨把误差估计为12*0.7=10分钟。
   为什么不把误差估计为7+5=12分钟？因为7与5的均指最大误差，出现最大误差的可能性已经很小了，二者同时出现最大值的可能性几乎为0。
   当然，平均误差仅为最大误差的六分之一左右。
   ▲第(1)步与第(2)步计算无需考虑月球光行时、地球章动的修正计算，因为误差远超过这些修正量。
(3)a.余项计算：dW = W - Lon(to,159);
   b.时间计算：to = to + dW/v;
   周期项取159项的理论精度是3",对应的时间误差为6秒。
   速度不精确引起的误差：10分钟*0.6%=3.6秒，平均误差按1/6计算为0.6秒，不过这个估值比实测值要大许多。
   总误差不妨按 6+0.6=7秒钟估计。
(4)如果相要再提高精度,继续计算：
   dW = W - Lon(to,全部);
   to = to+dW/v;
“全部”计算的精度为5毫角秒，相当于0.01秒
速度v不精确引起的误差：5秒*0.006=0.03秒
总误差不妨取0.03秒钟。
ELP/MPP02在J2000前后几十年内与DE406拟合得很好，差异为0.005角秒以内。
(5)说明：由于都是在邻近的时间内计算，可以认为速度不变(或者说加速度为0)。因为对ELP/MPP02求二次导数，我们发现加速度非常小。为什么会这样呢？其实，周期项的周期最小值为10天左右(个别为5天左右,但振幅非常小)，也就是说在10天范围内速度大小才小幅波动一次，在几小时或几分钟内，速度几乎不变。速度大小波动幅度最大的周期项的周期是27.5天左右，它对速度的影响可达10%以上，即便是这种周期项，18小时对速度的影响最多只有：10%*18小时/27.5天 = 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，以便计算其它时刻的情形。

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

xjw01于十中，2008年5月27日晚编写

M_L1=new Array(
1.677,4.669,628.30196,-0.03,0,0,
0.516,3.372,6585.7609,-2.16,-19,0.1,
0.414,5.728,14914.4523,-0.64,6,0,
0.371,3.969,7700.3895,1.55,25,0,
0.276,0.74,8956.9934,1.5,25,0,
0.246,4.23,-2.3012,1.5,25,0,
0.071,0.14,7842.365,-2.2,-19,0,
0.061,2.50,16171.056,-0.7,6,0,
0.045,0.44,8399.679,-0.4,0,0,
0.040,5.77,14286.150,-0.6,6,0,
0.037,4.63,1256.604,-0.1,0,0,
0.037,3.42,5957.459,-2.1,-19,0,
0.036,1.80,23243.144,0.9,31,0,
0.024,0,16029.081,3,50,0,
0.022,1.0,-1742.931,-4,-44,0,
0.019,3.1,17285.685,3,50,0,
0.017,1.3,0.329,2,25,0,
0.014,0.3,8326.390,3,50,0,
0.013,4.0,7072.088,2,25,0,
0.013,4.4,8330.993,0,0,0,
0.013,0.1,8470.667,-2,-19,0,
0.011,1.2,22128.515,-3,0,0,
0.011,2.5,15542.754,-1,0,0,
0.008,0.2,7214.06,-2,0,0,
0.007,4.9,24499.75,1,0,0,
0.007,5.1,13799.82,-4,0,0,
0.006,0.9,-486.33,-4,0,0,
0.006,0.7,9585.30,1,0,0,
0.006,4.1,8328.34,2,0,0,
0.006,0.6,8329.04,2,0,0,
0.005,2.5,-1952.48,1,0,0,
0.005,2.9,-0.71,0,0,0,
0.005,3.6,30457.21,-1,0,0);

M_L2=new Array(
0.0049,4.67,628.3020,0,0,0,
0.0023,2.67,-2.301,1.5,25,0,
0.0015,3.37,6585.761,-2.2,-19,0,
0.0012,5.73,14914.452,-0.6,6,0,
0.0011,3.97,7700.389,2,25,0,
0.0008,0.7,8956.993,1,25,0);

function Mnn(F,t,t2,t3,t4,n){ //计算M_L0或M_L1或M_L2
n = Math.floor(n * F.length / M_L0.length+0.5)*6; //项数是以周期项为基准,n是周期项的计算项数
if(n>F.length) n=F.length;
var i,v=0; t2/=1e4,t3/=1e8,t4/=1e8;
for(i=0;i<n;i+=6)
    v+=F[i]*Math.cos(F[i+1] +t*F[i+2] +t2*F[i+3] +t3*F[i+4] +t4*F[i+5]);
return v;
}
function moonLon(t,n){ //月球经度计算,返回Date分点黄经,传入世纪数
if(n==-1) n = Math.floor(M_L0.length/6+0.1);
var t2=t*t, t3=t2*t, t4=t3*t, t5=t4*t;
var v = Mnn( M_L0, t,t2,t3,t4, n );
      v += Mnn( M_L1, t,t2,t3,t4, n )*t;
      v += Mnn( M_L2, t,t2,t3,t4, n )*t2;
var w = 3.81034409 + 8399.684730072*t -3.319e-05*t2 + 3.11e-08*t3 - 2.033e-10*t4; //月球平黄经
var p =5028.792262*t + 1.1124406*t2 + 0.00007699*t3 - 0.000023479*t4 -0.0000000178*t5; //岁差
return w + (v+p)/rad;
}

function moonV(t,jing){ //月球速度计算,传入世经数
var v=8399.70911033384;
if(jing==0) return v;
v-=914*Math.sin( 0.7848 +8328.691425*t+ 1.523*t*t/10000 );
if(jing==1) return v;                    //误差小于5%
v-=179*Math.sin( 2.543 +15542.7543*t ) //误差小于0.3%
    +160*Math.sin( 0.1874 + 7214.0629*t )
    +62 *Math.sin( 3.14   +16657.3828*t )
    +34 *Math.sin( 4.827 +16866.9323*t )
    +22 *Math.sin( 4.9    +23871.4457*t )
    +12 *Math.sin( 2.59   +14914.4523*t )
    +7 *Math.sin(0.23    + 6585.7609*t )
    +5 *Math.sin(0.9     +25195.624 *t )
    +5 *Math.sin(2.32    - 7700.3895*t )
    +5 *Math.sin(3.88    + 8956.9934*t )
    +5 *Math.sin(0.49    + 7771.3771*t )
    +4 *Math.sin(5.5     +24986.074 *t )
    +4 *Math.sin(4.3     +22756.817 *t )
    +2 *Math.sin(1       +32200.137 *t )
    +2 *Math.sin(1.9     +14428.126 *t )
    +2 *Math.sin(0.4     +31085.509 *t )
    +2 *Math.sin(1.528   + 628.302 *t );
   return v;
}

function moon_t(W){ //已知黄经求时间
var t,dW,v= 8399.70911033384;
t = ( W - 3.81034 )/v; v=moonV(t,2); //v的精度0.6%，详见原文
t += ( W - moonLon(t,20) )/v;
t += ( W - moonLon(t,159))/v;
return t;
}

testT=30; //世纪数
L=moonLon(testT,-1); //正算
t=moon_t(L); //反算
dt=(testT-t)*35625*86400;

document.write("高速迭代法求指定Date平分点黄经的发生时刻。测试如下：<br>");
document.write("正算时间T=" + testT + "(世纪) 时的黄经是:" + L + "(弧度)<br>");
document.write("反算黄经L=" + L + "(弧度) 时的时间是:" + t + "(世纪)<br>");
document.write("迭代误差(不是实际误差):" +dt +"秒<br>");

</script>

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

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

月球光行差问题
我们已讨论过精密的光行差计算，但是，在农历计算中根本无需精密计算，以下考虑简化计算。
从ELP的中序列截取5%精度的距离及速速表达式：
    日月距离：r=385001 +20905*sin( 0.7848 +8328.691425*t+ 1.523*t*t/10000 ); //最大误差2%
    黄经速度：v=8400   - 914*sin( 0.7848 +8328.691425*t+ 1.523*t*t/10000 ); //最大误差5%
为了求证方便，令：
    r0=385001，v0=8400
    φ= 0.7848 +8328.691425*t+ 1.523*t*t/10000
    dr= 20905*sin(φ)
    dv=   914*sin(φ)
    光速：c=300000千米/秒=9.467E14千米/世纪
那么有：r = r0+dr；v = v0+dv；τ= r/c； (式中τ光行时)
光行期间黄经的变化量：
    dL=v*τ=v*r/c，代入得：
    dL = (v0+dv)*(r0+dr)/c = (v0*r0/c)*(1+dv/v+dr/r) (略去二阶小量)
    dL = (3.416E-6弧度)*(1+(20905/385001-914/8400)*sin(φ))
       = (0.7角秒)*(1-0.055*sin(φ))
结论：月球的黄经光行差约为0.7角秒(最大误差16%)，如果考虑上式中的周期，最大误差5%

以下关于章动，取自IAU1980，最大误差不超过0.1角秒
function nutation(t){ //章动计算,t为世纪数,返回值为角秒单位
var c, dL=0,dE=0, t2=t*t, y0=-17.20-0.01742*t;
c= 2.1824    -33.75705*t +36e-6*t2; dL = y0 *Math.sin(c); dE = 9.20*Math.cos(c);
c= 3.5069 +1256.66393*t +11e-6*t2; dL+= -1.32*Math.sin(c); dE+= 0.57*Math.cos(c);
c= 1.3375 +16799.4182*t -51e-6*t2; dL+= -0.23*Math.sin(c); dE+= 0.10*Math.cos(c);
c= 4.3649    -67.5141*t +72e-6*t2; dL+= 0.21*Math.sin(c); dE+=-0.09*Math.cos(c);
c= 0.04   -628.302*t; dL+= -0.14*Math.sin(c);
c= 2.36 +8328.691*t; dL+= 0.07*Math.sin(c);
c= 3.46 +1884.966*t; dL+= -0.05*Math.sin(c); dE+= 0.02*Math.cos(c);
c= 5.44 +16833.175*t; dL+= -0.04*Math.sin(c); dE+= 0.02*Math.cos(c);
c= 3.69 +25128.110*t; dL+= -0.03*Math.sin(c);
c= 3.55   +628.362*t; dL+= 0.02*Math.sin(c);
this.dL=dL; //黄经章动
this.dE=dE; //交角章动
}

调用方法：nu = new nutation(testT);

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

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

rad = 180*3600/Math.PI;
function rad2mrad(v){   //对超过0-2PI的角度转为0-2PI
v=v % (2*Math.PI);
if(v<0) return v+2*Math.PI;
return v;
}
function rad2str(d,tim){ //将弧度转为字串
//tim=0输出格式示例: -23°59' 48.23"
//tim=1输出格式示例: 18h 29m 44.52s
var s="+";
var w1="°",w2="’",w3="”";
if(d<0) d=-d,s='-';
if(tim){ d*=12/Math.PI; w1="h ",w2="m ",w3="s "; }
else     d*=180/Math.PI;
var a=Math.floor(d); d=(d-a)*60;
var b=Math.floor(d); d=(d-b)*60;
var c=Math.floor(d); d=(d-c)*100;
d=Math.floor(d+0.5);
if(d>=100) d-=100, c++;
if(c>=60) c-=60, b++;
if(b>=60) b-=60, a++;
a="   "+a, b="0"+b, c="0"+c, d="0"+d;
s+=a.substr(a.length-3,3)+w1;
s+=b.substr(b.length-2,2)+w2;
s+=c.substr(c.length-2,2)+".";
s+=d.substr(d.length-2,2)+w3;
return s;
}

testT=0.08+(1-1.5)/36525; //世纪数
L=moonLon(testT,-1); //正算
t=moon_t(L); //反算
dt=(testT-t)*35625*86400;

document.write("高速迭代法求指定Date平分点黄经的发生时刻。测试如下：<br>");
document.write("正算时间T=" + testT + "(世纪) 时的黄经是:" + L + "(度)<br>");
document.write("反算黄经L=" + L + "(度) 时的时间是:" + t + "(世纪)<br>");
document.write("迭代误差(不是实际误差):" +dt +"秒<br>");

//以下是视位置计算
function nutation(t){ //章动计算,t为世纪数,返回值为角秒单位
var c, dL=0,dE=0, t2=t*t, y0=-17.20-0.01742*t;
c= 2.1824    -33.75705*t +36e-6*t2; dL = y0 *Math.sin(c); dE = 9.20*Math.cos(c);
c= 3.5069 +1256.66393*t +11e-6*t2; dL+= -1.32*Math.sin(c); dE+= 0.57*Math.cos(c);
c= 1.3375 +16799.4182*t -51e-6*t2; dL+= -0.23*Math.sin(c); dE+= 0.10*Math.cos(c);
c= 4.3649    -67.5141*t +72e-6*t2; dL+= 0.21*Math.sin(c); dE+=-0.09*Math.cos(c);
c= 0.04   -628.302*t; dL+= -0.14*Math.sin(c);
c= 2.36 +8328.691*t; dL+= 0.07*Math.sin(c);
c= 3.46 +1884.966*t; dL+= -0.05*Math.sin(c); dE+= 0.02*Math.cos(c);
c= 5.44 +16833.175*t; dL+= -0.04*Math.sin(c); dE+= 0.02*Math.cos(c);
c= 3.69 +25128.110*t; dL+= -0.03*Math.sin(c);
c= 3.55   +628.362*t; dL+= 0.02*Math.sin(c);
this.dL=dL; //黄经章动
this.dE=dE; //交角章动
}

nu=new nutation(testT);
L2=rad2mrad(L+nu.dL/rad-0.7/rad); //补章动与光行差
Ls=rad2str(L2,0);
document.write("视黄经：" + Ls+"<br>");

测试结果如下：

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

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" 本程序

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

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

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

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

十、地球黄经为W时对应的时间to的求解
v(t) = 6283.32 //误差小于万分3
+210 *Math.sin(1.527+f)
+4.4 *Math.sin(1.48+f*2)
+12.9*Math.sin(5.82+f)*t
+0.55*Math.sin(4.21+f)*t*t;
式中f=6283.07585*t;
(1)解方程 W = Lon(t) = a+v*t-dW
方程右边是地球黄经表达式，a是它的常数项，v是黄经速度，dW是黄经的高次项(t的1次方以上的)和周期项，称dW为方程的余项。
因dW较小，所以先忽略dW，解得“平”时间：
to = (W - a)/v = (W-1.753469512)/6283.319653318
v = v(to)
▲to是地球黄经为W是的“平”时间，to与真时间的差值不超过48小时
▲余项引起的误差：或余项已知，to回代到原方程，得t = (W-a+dW)/v = to+dW/v
显然dW/v就是to的误差值。
把to代入余项的表达式就可得到余项的值。
此外，也可以Lon()函数觉察到余项：Lon(to,10) = a+v*to - dW = W + dW，即dW = W-Lon(to,10)
▲因为to不超过6儒略千年，所以平黄经的高次项引起的余项不超过1度，黄经周期项不超过2度，合计余项最多3度。地球每天运动1度左右，所以to最大误差3天。
v(t)求导得最大可能的加速度为(3.6弧度/千年)/天，如果to误差3天，那么v(t)算得的v值误差为3.6*3=10弧度/千年。
这样，v的相对误差为 10/6283=0.17%
(2)在步骤(1)留下的主要余项：dW = W-Lon(to,10);
   那么修正to为：
   to = to+ (W-Lon(to,10)) / v;
   ▲to的最大误差由v及Lon()函数的截断误差引起的。
   v引起的最大可能误差为 3天*0.17% = 7分钟
   Lon(to,10)截断误差为17角秒，地球每分钟运动2.5角秒，所以to误差7分钟左右。
   合误差可以估计为10分钟(以上两个误差同时发生最大值的可能性非常小，所以误差不必取7+7=14分钟)。
(3)执行to = to +(W-Lon(to,120)) /v;
   v引起的误差为(这就是最后的迭代误差): 10分钟*0.17%=1秒，平均误差应小于0.15秒
   Lon(to,120)的截断误差是0.25角秒，对应时间误差为6秒。

这是精简的超高速的算法

//以下是地球黄经数据,最大误差0.25"
var E_L0=new Array(
33416565,4.6692568,6283.07584999,
348943,4.626102,12566.1517,34971,2.74412,5753.38488,34176,2.82887,3.52312,31359,3.62767,77713.77147,
26762,4.41808,7860.41939,23427,6.13516,3930.2097,13243,0.74246,11506.76977,12732,2.0371,529.69097,
11992,1.10963,1577.34354,9903,5.2327,5884.9268,9019,2.0451,26.2983,8572,3.5085,398.149,
7798,1.1788,5223.6939,7531,2.5334,5507.5532,5053,4.5829,18849.2275,4924,4.2051,775.5226,
3567,2.9195,0.0673,3171,5.849,11790.6291,2841,1.8987,796.298,2710,0.3149,10977.0788,
2428,0.3448,5486.7778,2062,4.8065,2544.3144,2054,1.8695,5573.1428,2023,2.4577,6069.7768,
1555,0.8331,213.2991,1322,3.4112,2942.4634,1262,1.083,20.7754,1151,0.6454,0.9803,
1029,0.636,4694.003,1019,0.9757,15720.8388,1017,4.2668,7.1135,992,6.21,2146.165,
976,0.681,155.42,858,5.983,161000.686,851,1.299,6275.962,847,3.671,71430.696,
796,1.808,17260.155,788,3.037,12036.461,747,1.755,5088.629,739,3.503,3154.687,
735,4.679,801.821,696,0.833,9437.763,624,3.978,8827.39,611,1.818,7084.897,
570,2.784,6286.599,561,4.387,14143.495,556,3.47,6279.553,520,0.189,12139.554,
516,1.333,1748.016,511,0.283,5856.478,490,0.487,1194.447,410,5.368,8429.241,
409,2.399,19651.048,392,6.168,10447.388,368,6.041,10213.286,366,2.57,1059.382,
360,1.709,2352.866,356,1.776,6812.767,333,0.593,17789.846,304,0.443,83996.847,
300,2.74,1349.867,254,3.165,4690.48,247,0.215,3.59,237,0.485,8031.092,
236,2.065,3340.612,228,5.222,4705.732,219,5.556,553.569,214,1.426,16730.464,
211,4.148,951.718,203,0.371,283.859,199,5.222,12168.003,199,5.775,6309.374,
191,3.822,23581.258,189,5.386,149854.4,179,2.215,13367.973,175,4.561,135.065,
162,5.988,11769.854,151,4.196,6256.778,144,4.193,242.729,143,3.724,38.028,
140,4.401,6681.225,136,1.889,7632.943,125,1.131,5.523,121,2.622,955.6,
120,1.004,632.784,113,0.177,4164.312,108,0.327,103.093,105,0.939,11926.254,
105,5.359,1592.596,103,6.2,6438.496,100,6.029,5746.271,98,1,11371.7,
98,5.24,27511.47,94,2.62,5760.5,92,0.48,522.58,92,4.57,4292.33,
90,5.34,6386.17,86,4.17,7058.6,84,3.3,7234.79,84,4.54,25132.3,
81,6.11,4732.03,81,6.27,426.6,80,5.82,28.45,79,1,5643.18,
78,2.96,23013.54,77,3.12,7238.68,76,3.97,11499.66,73,4.39,316.39,
73,0.61,11513.88,72,4,74.78,71,0.32,263.08,68,5.91,90955.55,
66,3.66,17298.18,65,5.79,18073.7,63,4.72,6836.65,62,1.46,233141.31,
61,1.07,19804.83,60,3.32,6283.01,60,2.88,6283.14,55,2.45,12352.85);

var E_L1=new Array(
2060589,2.6782346,6283.07585,43034,2.63513,12566.1517,4253,1.5905,3.5231,1193,5.7956,26.2983,
1090,2.9662,1577.3435,935,2.592,18849.228,721,1.138,529.691,678,1.875,398.149,
673,4.409,5507.553,590,2.888,5223.694,560,2.175,155.42,454,0.398,796.298,
364,0.466,775.523,290,2.647,7.114,208,5.341,0.98,191,1.846,5486.778,
185,4.969,213.299,173,2.991,6275.962,162,0.032,2544.314,158,1.43,2146.165,
146,1.205,10977.079,125,2.834,1748.016,119,3.258,5088.629,118,5.274,1194.447,
115,2.075,4694.003,106,0.766,553.569,100,1.303,6286.599,97,4.24,1349.87,
95,2.7,242.73,86,5.64,951.72,76,5.3,2352.87,64,2.65,9437.76,
61,4.67,4690.48,58,1.77,1059.38,53,0.91,3154.69,52,5.66,71430.7,
52,1.85,801.82,50,1.42,6438.5,43,0.24,6812.77,43,0.77,10447.39,
41,5.24,7084.9,37,2,8031.09,36,2.43,14143.5,35,4.8,6279.55,
34,0.89,12036.46,34,3.86,1592.6,33,3.4,7632.94,32,0.62,8429.24,
32,3.19,4705.73,30,6.07,4292.33,30,1.43,5746.27,29,2.32,20.36,
27,0.93,5760.5,27,4.8,7234.79,25,6.22,6836.65,23,5,17789.85,
23,5.67,11499.66,21,5.2,11513.88,21,3.96,10213.29,21,2.27,522.58,
21,2.22,5856.48,21,2.55,25132.3,20,0.91,6256.78,19,0.53,3340.61,
19,4.74,83996.85,18,1.47,4164.31,18,3.02,5.52,18,3.03,5753.38,
16,4.64,3.29,16,6.12,5216.58,16,3.08,6681.22,15,4.2,13367.97,
14,1.19,3894.18,14,3.09,135.07,14,4.25,426.6,13,5.77,6040.35,
13,3.09,5643.18,13,2.09,6290.19,13,3.08,11926.25,12,3.45,536.8);

var E_L2=new Array(
87198,1.0721,6283.07585,3091,0.8673,12566.1517,273,0.053,3.523,163,5.188,26.298,
158,3.685,155.42,95,0.76,18849.23,89,2.06,77713.77,70,0.83,775.52,
51,4.66,1577.34,41,1.03,7.11,38,3.44,5573.14,35,5.14,796.3,
32,6.05,5507.55,30,1.19,242.73,29,6.12,529.69,27,0.31,398.15,
25,2.28,553.57,24,4.38,5223.69,21,3.75,0.98,17,0.9,951.72,
15,5.76,1349.87,14,4.36,1748.02,13,3.72,1194.45,13,2.95,6438.5,
12,2.97,2146.17,11,1.27,161000.69,10,0.6,3154.69,10,5.99,6286.6,
9,4.8,5088.63,9,5.23,7084.9,8,3.31,213.3,8,3.42,5486.78,
7,6.19,4690.48,7,3.43,4694,6,1.6,2544.31,6,1.98,801.82,
6,2.48,10977.08,5,1.44,6836.65,5,2.34,1592.6,5,1.31,4292.33,
5,3.81,149854.4,4,0.04,7234.79,4,4.94,7632.94,4,1.57,71430.7,
4,3.17,6309.37,3,0.99,6040.35,3,0.67,1059.38,3,3.18,2352.87,
3,3.55,8031.09,3,1.92,10447.39,3,2.52,6127.66,3,4.42,9437.76,
3,2.71,3894.18,3,0.67,25132.3,3,5.27,6812.77,3,0.55,6279.55);

var E_L3=new Array(
2892,5.8438,6283.0758,168,5.488,12566.152,30,5.2,155.42,
13,4.72,3.52,7,5.3,18849.23,6,5.97,242.73,4,3.79,553.57,1,4.3,6286.6,1,0.91,6127.66);

var E_L4=new Array(77,4.13,6283.08,8,3.84,12566.15,4,0.42,155.42);
var E_L5=new Array(2,2.77,6283.08,1,2.01,155.42);

//地球黄纬数据,误差0.2"
var E_B0=new Array(
2796,3.1987,84334.6616,1016,5.4225,5507.5532,804,3.88,5223.694,438,3.704,2352.866,
319,4,1577.344,227,3.985,1047.747);

var E_B1=new Array(90,3.9,5507.55,62,1.73,5223.69);
var E_B2=new Array(17,1.63,84334.66);

//地球向径数据,误差0.00001AU
var E_R0=new Array(
1000139888,0,0,16706996,3.09846351,6283.07584999,
139560,3.055246,12566.1517,30837,5.19847,77713.77147,
16285,1.17388,5753.38488,15756,2.84685,7860.41939,
9248,5.4529,11506.7698,5424,4.5641,3930.2097,
4721,3.661,5884.9268,3460,0.9637,5507.5532,
3288,5.8998,5223.6939,3068,0.2987,5573.1428,
2432,4.2735,11790.6291,2118,5.8471,1577.3435,
1858,5.0219,10977.0788,1748,3.0119,18849.2275);

var E_R1=new Array(1030186,1.1074897,6283.07585,17212,1.06442,12566.1517,7022,3.1416,0);
var E_R2=new Array(43594,5.78455,6283.07585,1236,5.5793,12566.1517,123,3.142,0,88,3.63,77713.77);
var E_R3=new Array(1446,4.2732,6283.0758,67,3.92,12566.15);
var E_R4=new Array(39,2.56,6283.08,3,2.27,12566.15);
var E_R5=new Array(1,1.22,6283.08);

//迭代算法测试
t=1;
L=earthLon(t,-1);
t2=earth_t(L);
document.write("迭代算法测试<br>");
document.write("输入时间(千年):"+t+"<br>");
document.write("地球黄经(弧度):"+L+"<br>");
document.write("返算时间(千年):"+t2+"<br>");
document.write("迭代误差(秒):"+(t2-t)*365250*86400+"<br><br>");

//=======================
// 以下是视位置计算
//=======================
var nutB=new Array(
2.1824, -33.75705, 36e-6,-1720,920,
3.5069, 1256.66393, 11e-6,-132, 57,
1.3375,16799.4182, -51e-6, -23, 10,
4.3649, -67.5141, 72e-6, 21, -9,
0.04,   -628.302,   0,     -14, 0,
2.36,   8328.691,   0,       7, 0,
3.46,   1884.966,   0,      -5, 2,
5.44, 16833.175,   0,      -4, 2,
3.69, 25128.110,   0,      -3, 0,
3.55,    628.362,   0,       2, 0);
function nutation(t){ //章动计算
var i,c,a, t2=t*t,dL=0,dE=0;
for(i=0;i<nutB.length;i+=5){
c = nutB[i]+nutB[i+1]*t+nutB[i+2]*t2;
if(i==0) a=-1.742*t; else a=0;
dL+=(nutB[i+3]+a)*Math.sin(c);
dE+= nutB[i+4]   *Math.cos(c);
}
this.dL=dL/100; //黄经章动
this.dE=dE/100; //交角章动
}
function gxc_sun(t,L){ //太阳光行差(黄经),L为太阳黄经,t是千年数
var v = L-(282.93735+17.1946*t+ 0.046*t*t)/180*Math.PI; //真近点角
var e = 0.016708634-0.00042037*t-0.00001267*t*t;
return -20.49552 * (1+e*Math.cos(v));
}
function HCconv(JW,E){ //黄赤转换(黄赤坐标旋转)
//黄道赤道坐标变换,赤到黄E取负
var HJ=rad2mrad(JW[0]), HW=JW[1];
var sinE =Math.sin(E),cosE =Math.cos(E);
var sinW=cosE*Math.sin(HW)+sinE*Math.cos(HW)*Math.sin(HJ);
var J=Math.atan2( Math.sin(HJ)*cosE-Math.tan(HW)*sinE, Math.cos(HJ) );
JW[0]=rad2mrad(J);
JW[1]=Math.asin(sinW);
}

function hcjj (t){ //返回黄赤交角,t是千年数
var t2=t*t, t3=t2*t,t4=t3*t;
return (84381.4088 -468.36051*t -0.01667*t2 -1.99911*t3-0.00523*t4)/rad;
}

z=new Array();
da=1;
testT=0.008+(da-1.5)/365250;    //世纪数
earthCoor(testT,z,-1);         //地球黄经
z[0]+=Math.PI,z[1]=-z[1];      //太阳黄经黄纬

//L=z[0]-50287.9*(testT-0.008)/rad;
//L=rad2mrad(L);
//document.write(rad2str(L,0)+"<br>");

gx= gxc_sun(testT,z[0]);       //太阳光行差
nu= new nutation(testT*10);    //章动计算
z[0] += (nu.dL+gx)/rad;        //补章动与光行差
E = hcjj(testT)+nu.dE/rad;     //得到真黄赤交角
HCconv(z,E); //转到赤道坐标

document.write("坐标精度测试<br>");
document.write("testT：" + testT+"<br>");
document.write("视赤经：" + rad2str(z[0],1)+"<br>");
document.write("视赤纬：" + rad2str(z[1],0)+"<br>");
document.write("距离：" + z[2]+"<br>");
document.write("请与《2008中国天文年历》太阳视赤经赤纬比对。该年的第 "+da+" 天0h TD。注意：几乎分秒不差！<br>");
document.write("本程序中VSOP87地球数据序列以做了长期项拟合到DE405/406<br>");

</script>

songgz

关注

13
点赞
踩
26

收藏

觉得还不错? 一键收藏
0
评论
现代历法（整理）

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

扫一扫