天文算法讨论 [中华农历论坛]

《天文算法》一书中给出了求各种天文现象发生时刻的方法，不过作者并没有深入讨论算法的效率、简洁性及具体原理。事实书中算法的精度是固定的，不一定满足现代天文计算的要求。
比如试图把农历的精度计算到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理论本身

第(8)贴《农历天文算法完全解——第一部分(分两部分讲述)》
算法本质：迭代技术
一、日月黄经速度的表达
(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=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方法

[此贴子已经被作者于2008-5-27 8:59:26编辑过]

第(8)贴的改进《农历天文算法完全解——第一部分》
算法本质：迭代技术
一、日月黄经速度的表达
(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秒。订正方法在前面的贴子中已经讲述。

[此贴子已经被作者于2008-5-28 7:23:55编辑过]

第(8)贴的算法实现(测试程)

这是一个超高速的算法，基于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" 本程序

[此贴子已经被作者于2008-5-28 15:29:01编辑过]

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

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

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

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

通过年的天干地支可以推导月的天干地支，知道日的天干地支就可以知道时的天干地支了，请问：怎样知道日的天干地支？

谢谢！

错误订正：

前面的贴子中的“地球速度”计算有误，请更正如下：

(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;

第(8)贴继续

十、地球黄经为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秒。

回复标题：
上传附件：

	签名：不显示显示

主题：天文算法讨论

天文算法讨论