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