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