为了便于天文计算,需要使用均匀流逝的时间标尺, 这种时间系统最早称为历书时, 后来,由于天文计算均以力学规律为基础, 又改称力学时。 人们发现原子内电子的能级跃迁振荡频率非常均匀, 如铯133的共振频率为每秒9192631770周。遵循这一原理, 科学家制作出了原子钟,作为均匀的时间标尺。因此, 原子时又成为力学时的等价时间。
通常, 人们用到的时间为世界时(UT)。 世界时是完全根据地球自转而确定的时间参考系统。 但是地球自转是不均匀的(章动), 同时也包括一些未知的因素, 造成了世界时与力学时之间的差异。
星历表、日月食、行星动态等, 都是先通过力学时为基准计算的, 再通过世界时与力学时之差转换为世界时。因此,
在天象计算中, Delta T是决定计算精度的重要因素之一。
但是Delta T只能通过实测来获得。 对远古时代的Delta T, 只能通过古代的可靠的天象记录进行推测, 例如, 2500年前, Delta T 可能达到4到5小时。
为此, 若需要某一特定时刻的Delta T, 则需要使用插值法。 插值法包括线性插值, 二次插值, 三次插值, 拉格朗日插值法等多种, 这里提供三次插值法的R函数, 以及相应的插值表, 以备参考。
详细参考 天文算法 http://www.fjptsz.com/xxjs/xjw/rj/117/03.htm
这里给出R函数, 计算Delta T
数据和函数改编自 许剑伟 老师的 寿星万年历 javascript 源代码。
int2 <- function (v) { floor**(v)**; }
dt_ext <- function**(y, jsd){** #二次曲线外推,用于数值外插
dy = **(y - 1820)/**100;
return (-20 + jsd ***** dy ***** dy);
}
dt_calc <- function**(y){** #传入年, 返回世界时UT与原子时(力学时 TD)之差, ΔT = TD - UT
y0 = dt_at**[length(dt_at)** - 1**]**; #表中最后一年
t0 = dt_at**[length(dt_at)]**; #表中最后一年的 deltatT
if**(y >= y0){**
jsd = 31; # sjd是y1年之后的加速度估计。
# 瑞士星历表jsd=31, NASA网站jsd=32, skmap的jsd=29
if**(y > y0 + 100){**
return**(dt_ext(y, jsd))**
};
v = dt_ext**(y, jsd)**; #二次曲线外推
dv = dt_ext**(y0, jsd)** - t0; # ye年的二次外推与te的差
return (v - dv(y0 + 100 - y)/100)*;
}
d = dt_at;
for**(i in seq(1, length(d), by = 5)** ) {
if**(y < d[i + 5])** break; # 判断年所在的区间
}
t1 = (y - d[i]) / (d[i + 5] - d**[i])** ***** 10 #### 三次插值, 保证精确性
t2 = t1 ***** t1
t3 = t2 ***** t1;
res <- d**[i + 1]** +d[i + 2] ***** t1 +d[i + 3] ***** t2 +d[i + 4] ***** t3
return**(** res );
}
# TD - UT1 插值表
dt_at = c**(**
**-**4000, 108371.7, **-**13036.80, 392.000, 0.0000,
**-**500, 17201.0, **-**627.82, 16.170, **-**0.3413,
**-**150, 12200.6, **-**346.41, 5.403, **-**0.1593,
150, 9113.8, **-**328.13, **-**1.647, 0.0377,
500, 5707.5, **-**391.41, 0.915, 0.3145,
900, 2203.4, **-**283.45, 13.034, **-**0.1778,
1300, 490.1, **-**57.35, 2.085, **-**0.0072,
1600, 120.0, **-**9.81, **-**1.532, 0.1403,
1700, 10.2, **-**0.91, 0.510, **-**0.0370,
1800, 13.4, **-**0.72, 0.202, **-**0.0193,
1830, 7.8, **-**1.81, 0.416, **-**0.0247,
1860, 8.3, **-**0.13, **-**0.406, 0.0292,
1880, **-**5.4, 0.32, **-**0.183, 0.0173,
1900, **-**2.3, 2.06, 0.169, **-**0.0135,
1920, 21.2, 1.69, **-**0.304, 0.0167,
1940, 24.2, 1.22, **-**0.064, 0.0031,
1960, 33.2, 0.51, 0.231, **-**0.0109,
1980, 51.0, 1.29, **-**0.026, 0.0032,
2000, 63.87, 0.1, 0, 0,
2005, 64.7, 0.4, 0, 0, #一次项记为x, 则 10x=0.4秒/年*(2015-2005), 解得x=0.4
2015, 69
);
## 输入年, 计算delta T
dt_calc**(1800)**
dt_calc**(1960)**
dt_calc**(1965)**
dt_calc**(2012)**