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


<p class="codepen" data-height="1200" data-default-tab="result" data-slug-hash="ByBOYgM" data-pen-title="グリニッジ恒星時 世界時 各日0時表" data-user="ikiuo" style="height: 1200px; box-sizing: border-box; display: flex; align-items: center; justify-content: center; border: 2px solid; margin: 1em 0; padding: 1em;">
  <span>See the Pen <a href="https://codepen.io/ikiuo/pen/ByBOYgM">
  グリニッジ恒星時 世界時 各日0時表</a> by Ikiuo (<a href="https://codepen.io/ikiuo">@ikiuo</a>)
  on <a href="https://codepen.io">CodePen</a>.</span>
</p>
<script async src="https://cpwebassets.codepen.io/assets/embed/ei.js"></script>

# 使用データ

計算用パラメータを求めるためのデータは、[理科年表](https://official.rikanenpyo.jp/)プレミアムの グリニッジ恒星時 1960 〜 2025 年(66年分)を[修正](https://qiita.com/ikiuo/items/fcd5624dc8585b82f64a)して使用しています。


# 計算式

計算式は、一次式とsin関数4つの和です。sin関数の周期は、理科年表に記載されている

| 周期 | 日 |
|:-:|:--|
| 半年 | 365.24219 ÷ 2 |
| １年 | 365.24219 |
| 月:対恒星交点順行周期 | 3232.605 |
| 月:対恒星交点逆行周期 | 6793.477 |

を用いています。パラメータは Python の numpy と scipy.optimize.curve_fit で求めました。

```python: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: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時の表

```html: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 = [
             "１", "２", "３", "４", "５", "６",
             "７", "８", "９", "１０", "１１", "１２",
         ];

         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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/226677/f841c563-25d4-422f-0d3c-07f122276e63.png)
1. 一次式との差分(1)：先頭 30 日
   元データが階段状なので一次式との差分では鋸波が出現する。
   ![Figure_2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/226677/c1fc8836-e95c-81ed-ca8e-8310a28ff32a.png)
1. 一次式との差分：全体
   月の対恒星交点逆行周期(約18.6年)と一致。青の凸凹が鋸波。
   ![Figure_3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/226677/10a92132-8291-0794-a96f-3817e95cd838.png)
1. 月の対恒星交点逆行周期を除去したもの
   1985年(パソコンが普及しはじめた頃)を境に段差(約0.05秒)が生じていた。
   0.1秒未満の丸め処理(手計算と機械計算)の違いかな？
   一次式+段差関数で一致させた。(計算式の三項演算子の原因)
   ![Figure_4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/226677/c9f9f375-9e99-b249-e478-8636cf2c78f3.png)
1. 段差を除去したもの
   オレンジは半年周期(下記)。
   ![Figure_5.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/226677/fa56e3ac-1805-397e-1613-5647c3b817eb.png)
   先頭 4年を見ると半年周期と一致。
   ![Figure_5a.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/226677/5d29d7cb-a2d0-7f71-105b-6b5311f6c290.png)
1. 半年周期を除去したもの
   オレンジは一年周期(下記)。
   ![Figure_6.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/226677/06baa991-992e-5c17-9439-9ec3d0245c62.png)
   先頭 4年を見ると一年周期らしきものがあるので除去してみる。
   ![Figure_6a.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/226677/1eb3672f-5940-ff48-4627-af8aa8d610cd.png)
1. 一年周期を除去したもの
   青の凸周期が月の対恒星交点順行周期(約8.85年)に思えるので除去してみる。
   ![Figure_7.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/226677/815cf5af-22dc-fafd-1c20-ea18ebdf62a3.png)
1. 残り
   ![Figure_8.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/226677/d8085479-3dc6-2cbd-d4e9-4cad161f1848.png)
   残り：先頭 30 日
   当たり前だけど、鋸波はしっかり残ってる...
   ![Figure_8a.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/226677/627f73c1-df02-4338-8731-67f9e606f7ad.png)

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

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


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

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

```Python: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: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進の<code>p±N</code>は2進数のシフト量です。

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