See the Pen グリニッジ恒星時 世界時 各日0時表 by Ikiuo (@ikiuo) on CodePen.
使用データ
計算用パラメータを求めるためのデータは、理科年表プレミアムの グリニッジ恒星時 1960 〜 2025 年(66年分)を修正して使用しています。
計算式
計算式は、一次式とsin関数4つの和です。sin関数の周期は、理科年表に記載されている
周期 | 日 |
---|---|
半年 | 365.24219 ÷ 2 |
1年 | 365.24219 |
月:対恒星交点順行周期 | 3232.605 |
月:対恒星交点逆行周期 | 6793.477 |
を用いています。パラメータは Python の numpy と scipy.optimize.curve_fit で求めました。
# import math が必要
# T は 2000/1/1 からの経過日数で、結果は秒[0,86400)が得られる
(sum([
(23992.216182023083 if T < -5478 else 23992.272751683977),
236.55536717968295 * T,
*(v * math.sin((T + p) * r) for v, p, r in (
(0.07913717139970983, 14619.815697014918, 0.03440558335924766),
(0.007770275860376258, 14601.259015974265, 0.01720279167962383),
(-0.01008260911585117, 14050.135958804629, 0.0019436910192181186),
(-1.0558999819068642, 14628.993680660877, 0.0009248850488755002)))
]) % 86400)
// T は 2000/1/1 からの経過日数で、結果は秒(-86400,86400)が得られる
((((T < -5478) ? 23992.216182023083 : 23992.272751683977)
+ 236.55536717968295 * T
+ [[0.07913717139970983, 14619.815697014918, 0.03440558335924766],
[0.007770275860376258, 14601.259015974265, 0.01720279167962383],
[-0.01008260911585117, 14050.135958804629, 0.0019436910192181186],
[-1.0558999819068642, 14628.993680660877, 0.0009248850488755002]
].reduce((a, d) => a + d[0] * Math.sin((T + d[1]) * d[2]), 0)) % 86400)
誤差は±0.1秒以内に収まりましたが、数値表現では理科年表と違うものが出てきます。
グリニッジ恒星時 世界時 0時の表
<!DOCTYPE html>
<html lang="ja">
<head>
<meta charset="utf-8">
<title>グリニッジ恒星時 世界時 各日0時</title>
<style>
table {
border: solid 1px lightgray;
border-collaspe: collaspe;
border-spacing: 0;
}
th, td {
border: solid 1px lightgray;
padding: 4px;
}
</style>
</head>
<body>
<script>
const start_year = 1950;
const end_year = 2050;
const num_years = end_year - start_year + 1;
const table = [...Array(num_years)].map(() => [...Array(12)].map(() => []));
function initialize() {
const edt = Date.UTC(end_year + 1);
var cdt = Date.UTC(start_year);
for (var p = 0; cdt < edt; p += 1, cdt += 86400_000) {
const dt = new Date(cdt);
const year = dt.getUTCFullYear();
const month = dt.getUTCMonth();
const ytab = table[year - start_year];
const mtab = ytab[month];
mtab.push(new Date(cdt - (9*3600_000)));
}
const year = Math.max(start_year, Math.min(new Date().getFullYear(), end_year));
document.body.innerHTML = [
'<table style="margin-bottom: 2px;">',
'<tr>',
'<td>グリニッジ恒星時 世界時 各日0時</td>',
'<td>',
'<input',
/**/ ' id="tagYear" name="tagYear" type="number"',
/**/ ` min="${start_year}" max="${end_year}" value="${year}"`,
/**/ ' onchange="setYear()">',
'<label for="tagYear">年</label>',
`<span>(${start_year}〜${end_year})</span>`,
'</td>',
'</table>',
'<table id="tagTable"></table>',
].join('');
setYear();
}
function getGST(T) {
return ((((T < -5478) ? 23992.216182023083 : 23992.272751683977)
+ 236.55536717968295 * T
+ [[0.07913717139970983, 14619.815697014918, 0.03440558335924766],
[0.007770275860376258, 14601.259015974265, 0.01720279167962383],
[-0.01008260911585117, 14050.135958804629, 0.0019436910192181186],
[-1.0558999819068642, 14628.993680660877, 0.0009248850488755002]
].reduce((a, d) => a + d[0] * Math.sin((T + d[1]) * d[2]), 0)) % 86400);
}
function setYear() {
const dty2000 = new Date(Date.UTC(2000, 1-1, 1) - (9*3600_000));
const month = [
"1", "2", "3", "4", "5", "6",
"7", "8", "9", "10", "11", "12",
];
const year = Math.max(start_year, Math.min(tagYear.value, end_year));
if (tagYear.value != year)
tagYear.value = year;
const ytab = table[year - start_year];
const fmt02 = ((n) => {
const s = `0${n}`;
return s.substring(s.length-2);
});
const fmt02_1f = ((n) => `${fmt02((n/10)|0)}.${n%10}`);
const fmt02_2f = ((n) => `${fmt02((n/100)|0)}.${fmt02(n%100)}`);
const hms = ((msec) => {
if (msec < 0) msec = (msec % 86400_000) + 86400_000;
const h = Math.trunc(msec / 3600_000); msec -= h * 3600_000;
const m = Math.trunc(msec / 60_000); msec -= m * 60_000;
return `${fmt02(h)}:${fmt02(m)}:${fmt02_1f((msec/100)|0)}`;
// return `${fmt02(h)}:${fmt02(m)}:${fmt02_2f((msec/10)|0)}`;
});
const getData = ((mon, day) => {
const mtab = ytab[mon];
return (day >= mtab.length) ? ' ' :
hms(getGST((mtab[day] - dty2000) / 86400_000) * 1000);
});
const html = [
'<tr>',
'<th>日\月</th>',
[...Array(12)].map((_, n) => `<th>${month[n]}</th>`),
[...Array(31)].map((_, day) => [
'<tr>',
`<th>${day+1}</th>`,
[...Array(12)].map((_, mon) => [
'<td>', getData(mon, day), '</td>',
].join('')),
'</tr>',
].flat().join('')),
].flat().join('');
tagTable.innerHTML = html;
}
window.onload = initialize;
</script>
</body>
</html>
誤差の除去手順
- データから前日との差分を累積 (1960/1/2〜2025/12/31)
前日との差分(秒)は236.5秒または236.6秒。一次式は係数は236.5〜236.6の間だから、さらにデータ列から236.5秒を引いた0.0秒と0.1秒を累積。地球の自転速度との関係なので全体像では直線に見える。
グラフは matplotlib.pyplot で作成し、縦は秒、横は 1960/1/2 からの日数です。
- 一次式との差分(1):先頭 30 日
元データが階段状なので一次式との差分では鋸波が出現する。
- 一次式との差分:全体
月の対恒星交点逆行周期(約18.6年)と一致。青の凸凹が鋸波。
- 月の対恒星交点逆行周期を除去したもの
1985年(パソコンが普及しはじめた頃)を境に段差(約0.05秒)が生じていた。
0.1秒未満の丸め処理(手計算と機械計算)の違いかな?
一次式+段差関数で一致させた。(計算式の三項演算子の原因)
- 段差を除去したもの
オレンジは半年周期(下記)。
先頭 4年を見ると半年周期と一致。
- 半年周期を除去したもの
オレンジは一年周期(下記)。
先頭 4年を見ると一年周期らしきものがあるので除去してみる。
- 一年周期を除去したもの
青の凸周期が月の対恒星交点順行周期(約8.85年)に思えるので除去してみる。
- 残り
残り:先頭 30 日
当たり前だけど、鋸波はしっかり残ってる...
周期関数部分の誤差は、地球の自転速度のふらつき?
計算結果と元データの差がプラスだと表示は一致しますが、マイナスだと0.1秒前で表示されます。誤差は±0.1秒の範囲に収まっていますが、±0.05秒未満なら四捨五入で表示も一致させられたでしょう。しかし、鋸波の高さが約0.05秒あります。ぐぅ...
パラメータが IEEE 754 厳密値の旧計算式
旧計算式のパラメータは下記の Python float, JavaScript Number 型である IEEE 754 倍精度浮動小数を正確に10進表現した厳密値を使用していました。16進表現にすると機械計算の演算精度に収まることが確認できます。
(sum([
(+23992.21618202308309264481067657470703125
if T < -5478 else
+23992.27275168397682136856019496917724609375),
+236.55536717968294624370173551142215728759765625 * T,
*(v * math.sin((T + p) * r) for v, p, r in (
(+0.079137171399709826946633484112680889666080474853515625,
+14619.81569701491753221489489078521728515625,
+0.03440558335924766286684217675428953953087329864501953125),
(+0.00777027586037625818404794841853799880482256412506103515625,
+14601.25901597426491207443177700042724609375,
+0.017202791679623831433421088377144769765436649322509765625),
(-0.01008260911585116999245048674538338673301041126251220703125,
+14050.13595880462889908812940120697021484375,
+0.00194369101921811862805633541739780412171967327594757080078125),
(-1.0558999819068641823349707919987849891185760498046875,
+14628.993680660876634647138416767120361328125,
+0.000924885048875500183178377522352775486069731414318084716796875)))
]) % 86400)
現計算式のパラメータは Python float, JavaScript Number では正確に表現できない近似値(2進/16進の厳密値は循環小数)ですが、実行時は厳密値と同じ float/Number 値に変換されることを期待しています。
(sum([
(23992.216182023083 if T < -5478 else 23992.272751683977),
236.55536717968295 * T,
*(v * math.sin((T + p) * r) for v, p, r in (
(0.07913717139970983, 14619.815697014918, 0.03440558335924766),
(0.007770275860376258, 14601.259015974265, 0.01720279167962383),
(-0.01008260911585117, 14050.135958804629, 0.0019436910192181186),
(-1.0558999819068642, 14628.993680660877, 0.0009248850488755002)))
]) % 86400)
各パラメータを、16進,2進,10進,10進近似の順で並べて確認します。16進のp±N
は2進数のシフト量です。
------------------
一次式のパラメータ
------------------
0x1.76e0dd5ed1fc8p+14
1. 7 6 e 0 d d 5 e d 1 f c 8
101110110111000.00110111010101111011010001111111001000
23992.21618202308309264481067657470703125
23992.216182023083
0x1.76e1174c37aa7p+14
1. 7 6 e 1 1 7 4 c 3 7 a a 7
101110110111000.01000101110100110000110111101010100111
23992.27275168397682136856019496917724609375
23992.272751683977
0x1.d91c591644052p+7
1. d 9 1 c 5 9 1 6 4 4 0 5 2
11101100.100011100010110010001011001000100000001010010
236.55536717968294624370173551142215728759765625
236.55536717968295
--------------------
半年周期のパラメータ
--------------------
0x1.442556b0f4884p-4
1. 4 4 2 5 5 6 b 0 f 4 8 8 4
0.00010100010000100101010101101011000011110100100010000100
0.079137171399709826946633484112680889666080474853515625
0.07913717139970983
0x1.c8de868c28142p+13
1. c 8 d e 8 6 8 c 2 8 1 4 2
11100100011011.110100001101000110000101000000101000010
14619.81569701491753221489489078521728515625
14619.815697014918
0x1.19d9bcea7d18ap-5
1. 1 9 d 9 b c e a 7 d 1 8 a
0.000010001100111011001101111001110101001111101000110001010
0.03440558335924766286684217675428953953087329864501953125
0.03440558335924766
--------------------
一年周期のパラメータ
--------------------
0x1.fd3b98b382c7ap-8
1. f d 3 b 9 8 b 3 8 2 c 7 a
0.000000011111110100111011100110001011001110000010110001111010
0.00777027586037625818404794841853799880482256412506103515625
0.007770275860376258
0x1.c84a1276f794ep+13
1. c 8 4 a 1 2 7 6 f 7 9 4 e
11100100001001.010000100100111011011110111100101001110
14601.25901597426491207443177700042724609375
14601.259015974265
0x1.19d9bcea7d18ap-6
1. 1 9 d 9 b c e a 7 d 1 8 a
0.0000010001100111011001101111001110101001111101000110001010
0.017202791679623831433421088377144769765436649322509765625
0.01720279167962383
---------------------------------
月:対恒星交点順行周期のパラメータ
---------------------------------
0x1.4a630e34996d9p-7
1. 4 a 6 3 0 e 3 4 9 9 6 d 9
0.00000010100101001100011000011100011010010011001011011011001
0.01008260911585116999245048674538338673301041126251220703125
0.01008260911585117
0x1.b711167191dbep+13
1. b 7 1 1 1 6 7 1 9 1 d b e
11011011100010.001000101100111000110010001110110111110
14050.13595880462889908812940120697021484375
14050.135958804629
0x1.fd86e571bc865p-10
1. f d 8 6 e 5 7 1 b c 8 6 5
0.00000000011111110110000110111001010111000110111100100001100101
0.00194369101921811862805633541739780412171967327594757080078125
0.0019436910192181186
---------------------------------
月:対恒星交点逆行周期のパラメータ
---------------------------------
0x1.0e4f76122310bp+0
1. 0 e 4 f 7 6 1 2 2 3 1 0 b
1.0000111001001111011101100001001000100011000100001011
1.0558999819068641823349707919987849891185760498046875
1.0558999819068642
0x1.c927f30ed8a91p+13
1. c 9 2 7 f 3 0 e d 8 a 9 1
11100100100100.111111100110000111011011000101010010001
14628.993680660876634647138416767120361328125
14628.993680660877
0x1.e4e7f84cc21f7p-11
1. e 4 e 7 f 8 4 c c 2 1 f 7
0.000000000011110010011100111111110000100110011000010000111110111
0.000924885048875500183178377522352775486069731414318084716796875
0.0009248850488755002
2進数で $1$ のある最下位桁に対応する10進数は必ず $5$ です。