世界时与力学时之差:Delta T的计算

为了便于天文计算,需要使用均匀流逝的时间标尺, 这种时间系统最早称为历书时, 后来,由于天文计算均以力学规律为基础, 又改称力学时。 人们发现原子内电子的能级跃迁振荡频率非常均匀, 如铯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)**