1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
|
#include "calendar.h"
/* 需要注意的是日期与儒略日数的转换
* 儒略日/世纪数/千年数不是以日计算的,而是准确到具体时刻
* 因而不使用整数
* 一般地是采用公式,但不注意的话会有数据类型的问题
* 这里方便起见采用Unix时间戳,即计算机系统提供的从1970年1月1日0时到指定时刻的秒数
* Unix时间戳规定为32位整数,因而最大值为2038年1月19日3时14分7秒,最小值为1901年12月13日20时45分52秒
* 但是我目前提供的中time_t类型为64位整数,因而可以表示数千年的时间范围
*/
double Julian::d[23][5] = {
{-4000, 108371.700000, -13036.800000, 392.000000, 0.000000},
{-500, 17201.000000, -627.820000, 16.170000, -0.341300},
{-150, 12200.600000, -346.410000, 5.403000, -0.159300},
{150, 9113.800000, -328.130000, -1.647000, 0.037700},
{500, 5707.500000, -391.410000, 0.915000, 0.314500},
{900, 2203.400000, -283.450000, 13.034000, -0.177800},
{1300, 490.100000, -57.350000, 2.085000, -0.007200},
{1600, 120.000000, -9.810000, -1.532000, 0.140300},
{1700, 10.200000, -0.910000, 0.510000, -0.037000},
{1800, 13.400000, -0.720000, 0.202000, -0.019300},
{1830, 7.800000, -1.810000, 0.416000, -0.024700},
{1860, 8.300000, -0.130000, -0.406000, 0.029200},
{1880, -5.400000, 0.320000, -0.183000, 0.017300},
{1900, -2.300000, 2.060000, 0.169000, -0.013500},
{1920, 21.200000, 1.690000, -0.304000, 0.016700},
{1940, 24.200000, 1.220000, -0.064000, 0.003100},
{1960, 33.200000, 0.510000, 0.231000, -0.010900},
{1980, 51.000000, 1.290000, -0.026000, 0.003200},
{2000, 63.870000, 0.100000, 0.000000, 0.000000},
{2005, 64.700000, 0.210000, 0.000000, 0.000000},
{2012, 66.800000, 0.220000, 0.000000, 0.000000},
{2018, 69.000000, 0.360000, 0.000000, 0.000000},
{2028, 72.600000, 0.000000, 0.000000, 0.000000}};
double Julian::dt_ext(int y, double jsd) {
double dy = (double)(y - 1820) / 100.0;
return -20 + jsd * dy * dy;
}
//计算力学时与世界时之差,传入年份
double Julian::delta_t(int y) {
int y0 = d[22][0]; //表中最后一行的年份
double t0 = d[22][1]; //表中最后一行的delta_t
double jsd;
if (y > y0) {
jsd = 31;
if (y - y0 > 100) {
return dt_ext(y, jsd);
}
double v = dt_ext(y, jsd);
double dv = dt_ext(y0, jsd) - t0;
return v - dv * (y0 + 100 - y) / 100;
}
double res;
for (int i = 0; i < 22; i++) {
if (y < d[i + 1][0]) {
break;
}
double t1 = (y - d[i][0]) / (d[i + 1][0] - d[i][0]) * 10;
double t2 = t1 * t1;
double t3 = t2 * t1;
res = d[i][1] + d[i][2] * t1 + d[i][3] * t2 + d[i][4] * t3;
}
return res;
}
// 计算儒略日
double Julian::getJulianDay(time_t time) {
double t = (double)time;
// 2440587.5为1970年1月1日0时的儒略日数
return (double)t / 86400.0 + 2440587.5;
}
// 计算儒略千年数,自JD2000(即2000年1月1日12时)起算
double Julian::getJulianKiloYear(time_t time) {
double jd = getJulianDay(time);
// 2451545.0为2000年1月1日12时的儒略日数,365250.0为一千年的天数
return (double)(jd - 2451545.0) / 365250.0;
}
//儒略千年数转时间戳
time_t Julian::kiloYearToTime(double t, int year) {
double jd = t * 365250.0 + 2451545.0;
time_t time = (time_t)((jd - 2440587.5) * 86400.0 - delta_t(year));
return time;
}
|