1
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

グリニッジ恒星時 世界時 0時表

Last updated at Posted at 2025-01-21

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 で求めました。

Python
# 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)
JavaScript
// 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時の表

sample.html
<!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 = [
             "", "", "", "", "", "",
             "", "", "", "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) ? '&nbsp;' :
                    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>

誤差の除去手順

  1. データから前日との差分を累積 (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 からの日数です。
    Figure_1.png
  2. 一次式との差分(1):先頭 30 日
    元データが階段状なので一次式との差分では鋸波が出現する。
    Figure_2.png
  3. 一次式との差分:全体
    月の対恒星交点逆行周期(約18.6年)と一致。青の凸凹が鋸波。
    Figure_3.png
  4. 月の対恒星交点逆行周期を除去したもの
    1985年(パソコンが普及しはじめた頃)を境に段差(約0.05秒)が生じていた。
    0.1秒未満の丸め処理(手計算と機械計算)の違いかな?
    一次式+段差関数で一致させた。(計算式の三項演算子の原因)
    Figure_4.png
  5. 段差を除去したもの
    オレンジは半年周期(下記)。
    Figure_5.png
    先頭 4年を見ると半年周期と一致。
    Figure_5a.png
  6. 半年周期を除去したもの
    オレンジは一年周期(下記)。
    Figure_6.png
    先頭 4年を見ると一年周期らしきものがあるので除去してみる。
    Figure_6a.png
  7. 一年周期を除去したもの
    青の凸周期が月の対恒星交点順行周期(約8.85年)に思えるので除去してみる。
    Figure_7.png
  8. 残り
    Figure_8.png
    残り:先頭 30 日
    当たり前だけど、鋸波はしっかり残ってる...
    Figure_8a.png

周期関数部分の誤差は、地球の自転速度のふらつき?

計算結果と元データの差がプラスだと表示は一致しますが、マイナスだと0.1秒前で表示されます。誤差は±0.1秒の範囲に収まっていますが、±0.05秒未満なら四捨五入で表示も一致させられたでしょう。しかし、鋸波の高さが約0.05秒あります。ぐぅ...

パラメータが IEEE 754 厳密値の旧計算式

旧計算式のパラメータは下記の Python float, JavaScript Number 型である IEEE 754 倍精度浮動小数を正確に10進表現した厳密値を使用していました。16進表現にすると機械計算の演算精度に収まることが確認できます。

Python(旧)
(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 値に変換されることを期待しています。

Python(現)
(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進数のシフト量です。

16進,2進,10進,10進近似
------------------
一次式のパラメータ
------------------

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$ です。

1
0
3

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
1
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?