1 : 21 春分間隔から春分年へ

← 前 1‒20 p ｜ ↑ もくじ i ｜ 次 1‒22 n →

「春夏秋冬」は「夏秋冬春」より長い

2017年11月26日
記事ID e71126

今の時代、地球の四季の循環周期（回帰年）は、365日5時間48分45秒。「春分→次の春分」の平均間隔はそれより長く、「夏至→夏至」はそれより短い。

例えば「春分→夏至」「夏至→秋分」「秋分→冬至」「冬至→春分」が、それぞれ平均93日・94日・90日・89日とすると、どこから数え始めても1周は 93 + 94 + 90 + 89 日で、ぴったり同じ所要時間になるように思える。なぜそうならないのだろう？

回帰年はだんだん短くなっているのに、春分間隔・春分年がだんだん長くなっているのは、なぜ？

PNG画像 (20 KiB): 春分年・夏至年・秋分年・冬至年の長さの変動は、位相がずれたサインカーブのような曲線を描く。
  1. §1　なぜ夏至→夏至は速いのか　（公転最速理論？）
  2. §2　恒星年・回帰年　（もしも地球の軌道が円だったら）
  3. §3　春分の間隔　（楕円軌道の素朴な計算）
  4. §4　解析的な解法　（瞬時の春分年）
  5. §5　考察
  6. 付録A　実装例

§1　なぜ夏至→夏至は速いのか　（公転最速理論？）

地球の軌道は、ほぼ楕円。太陽 S は（楕円の中心ではなく）楕円の焦点という場所にある。地球の公転（角速度）は、太陽のそばでは速い（太陽の重力が強く働くので）。太陽から離れるにつれて遅くなり、反対側では一番遅くなる。

とはいえ、軌道を完全に1周するのなら、どこからスタートしても、1周の所要時間にほとんど違いは生じない。高速区間と低速区間を全部1回ずつ通るからだ。

実際には、地球上で四季が一巡する期間（つまりは1年間）において、地球は、軌道をちゃんと 360° 公転しない！

太陽系の外から見ると、地球は1年につき 360° より少ない角度しか公転しないのである。…その理由は後回しにして、そこから生じる結果を考えてみよう。

四季の一巡に対応する公転角が (360 − x)° だったとしよう。つまり、360° より x° だけ公転量が「節約」された、とする。その場合:

A1 → A2 では高速区間をフルに通過するが、低速区間は一部省略される。B1 → B2 では低速区間を全部通過する羽目になるが、高速区間は一部省略される。A1 → A2 の方が、所要時間が短いことは明らかだろう。

現在の地球は、夏至のころ A1 側にいて、冬至のころ B1 側にいる。だから、夏至→次の夏至は時間が短いが、冬至→次の冬至は時間が長い。春分も、どちらかというと B1 寄りなので、春分→次の春分も、時間がやや長い。

表1: 春分・夏至などの平均間隔（365日5時間○○分××秒）
西暦年 春分 夏至 秋分 冬至
1000 48分51秒 47分58秒 48分51秒 49分41秒
2000 49分01秒 47分57秒 48分30秒 49分33秒
3000 49分09秒 47分59秒 48分12秒 49分20秒
4000 49分13秒 48分06秒 47分56秒 49分04秒

表1で、例えば「2000年」の春分「49分01秒」は「365日5時間49分01秒」の略。単位の「日・時・分・秒」は日常の単位と同じだが、うるう秒は ないものとする。

冬至→冬至は一番時間がかかる。春分→春分は少し時間がかかる。夏至→夏至は一番速い。けれど、この傾向は少しずつ変化している。冬至→冬至の所要時間はだんだん短縮され、春分→春分の所要時間はだんだん長くなっている。

これは「軌道上のどこで、どういう季節になるか」という関係が少しずつシフトするため。「そこを起点とする1年の長さ」を考えた場合、「春分」は、現在の「冬至」の立場（一番遅い）に近づきつつある。現在最速の「夏至」は、現在の「春分」の立場（少し時間がかかる）に近づきつつある。「秋分」は最速に近づき、現在いちばん遅い「冬至」は「少しだけ遅い」方向へ。

しかし、なぜ1年で (360 − x)° しか公転しないのだろうか。

次のように考えると、分かりやすい。

  1. 季節の変化は、地軸の傾きと太陽の位置関係によって決まる。例えば、南極側（地球の「下半分」）が太陽に向かってせり出している位置関係だと、南半球は日当たり良好で夏になり、北半球は冬になる。逆に北極側が太陽に向かってせり出していると、北半球が夏になる。こうした位置関係の変化の大部分は、もちろん地球の公転によって、1年周期で起きる。
  2. でも…。仮に地球が一切公転せず軌道上で止まっているとしても、もし地軸の向きが勝手に地球側で変わるなら（地球の「斜め具合」が変化して、南極側がせり出したり、北極側がせり出したりするなら）、やはり季節変化のようなことが起きるだろう。
  3. 現実には、季節の変化の大部分は公転によるのだが、地軸の向きも少し変化している。その少しの変化も、季節の進行に多少貢献している。そのため、完全に 360° 公転する少し手前で、季節は完全に一巡する。

正確に言うと、地軸の向きのゆるやかな変動（歳差現象）から説明されるのは x の8割くらいで、残りの約2割は「近日点（きんじつてん）移動」という別の現象による。

https://commons.wikimedia.org/wiki/File:Apsidendrehung.png

近日点移動は、地球の軌道の楕円自体がクルクル回る現象。公転する地球は、太陽系の外から見ると「エスカレーターの上を歩いている人」のような立場にある。エスカレーターの働きにより、外から見ると「自力で歩いている以上に進む」。言い換えると「その分、自力で歩いていない＝軌道上の公転量が少なくなる」。x の2割程度は、この現象の結果。

具体的な数値例（後述の方法で計算した「ケプラー的」平均日数）:

  1. 2000年の春分→夏至 a = 92.757822
  2. 2000年の夏至→秋分 b = 93.649445
  3. 2000年の秋分→冬至 c = 89.842381
  4. 2000年の冬至→春分 d = 88.992727
  5. 2001年の春分→夏至 e = 92.757074
  6. 2001年の夏至→秋分 f = 93.649836
  7. 2001年の秋分→冬至 g = 89.843104
  8. 2001年の冬至→春分 h = 88.992361

2000年春分→2001年春分は a + b + c + d = 365.242375 となり、2000年夏至→2001年夏至は b + c + d + e = 365.241627 となる。同じ「4個の季節1回ずつ」といっても、この場合「春」の長さ（a と e）が違うので、「足せば同じ」ではない。

§2　恒星年・回帰年　（もしも地球の軌道が円だったら）

数学的結論だけ知りたい方は、§4へ。

恒星年

地球軌道が円だったとして、中心から見て地球は1日1°動くとする。位置を「0°～360°」で表すとして、ある日の地球の位置が100°なら、翌日の位置は101°、翌々日の位置は102°、一般に d 日後の位置は 100 + d となる。1年＝365.25日とすれば、y 年後の位置は 100 + 365.25y となるだろう。

その単純なことを精密に表現したのが、次の式（出典 [1]）。

t は2000年1月1日12時からの経過時間（2000年1月1日12時より前なら負）で、単位は1000年。t = 1 は1000年後、t = −0.2 は200年前…などとなる。この式や、[1] の他の式は、t = −6 から t = +6 の範囲で有効（紀元前4000年から紀元後8000年くらい）。

右辺の単位は角度の「°」。3600で割っているので、1次の項の係数は約 359993.7° に当たる。1年で約360°動くから、1000年では1000周36万°というわけ。

単位の「1000年」は「365250日」と約束する（1年＝365.25日）。この日数は計算上の約束で、現実の暦の日数とは異なる。時間・時刻の基準は、とりあえず「日本時間マイナス9時間」と考えていい。

Ex2.1
2000年3月20日12時の地球位置を求めよう。この時刻は、2000年1月1日12時のちょうど 31 + 29 + 19 = 79 日後だから、時間変数は
  • t = 79/365250
これを公式にあてはめると:
  • earth(t) = 178°.3295…
ざっと「178°」と分かった。

この角度は、黄道の星座の中の特定の場所を基準に、黄道の星座に沿ってぐるっと測ったもの。太陽から見て地球が「178°」ということは、地球から見て太陽は「358°」の方角にあったことになる（「太陽から見た地球」と「地球から見た太陽」は180°反対方向）。

地球の軌道は本当は円ではないので、earth(t) は理論上の「架空の地球」の位置にすぎない。でも、この理論的位置は「公転周期が現実の地球と平均的に一致する」ように設定されていて、いろいろな計算の基礎となる。

earth(t) を t について微分【※】して365250で割ると:

  1. rate_earth(t) = (1295977422.834 29 − 2 × 2.044 11t − 3 × 0.005 23t2) / 3600 / 365250
  2. = (1295977422.834 29 − 2 × 2.044 11t − 3 × 0.005 23t2) / (3600 × 365250)

これは、時刻 t における「黄道の星座に対する、地球の公転の速さ」に当たる（単位: °/日）。365250で割ったのは、1000年当たりの角速度を1日当たりの角速度に換算するため。

【※】 「位置」を表す式を微分すると「速度」を表す式が得られる。微分の計算の仕方は次の通り。もともとの定数項は消滅し、もともとの1次の項の係数が、新たに定数項となる。もともとの2次の項 At2 は 2At に変わる。もともとの3次の項 Bt3 は 3Bt2 に変わる。要するに、それぞれの項について、係数を「もともとの指数」倍して、もともとの指数から1引けばいい。

t = 0 のとき:

これは、2000年1月1日12時における「地球の公転の、背景の星座に対する角速度」。地球から見れば、「黄道の星座の中を動いていく太陽の、位置の変化率」（1日当たりの角度）。360°をこの値で割ると、「太陽が黄道の星座の中を1周するのに必要な日数」が得られる:

これが、2000年1月1日12時における恒星年の長さ（星座から見て、地球が太陽の周りを1周する所要時間）。

r は瞬間速度で、D 日間における平均速度ではない。つまり、地球が D 日間公転したときの移動量は、正確に360°とは限らない。「今、時速100キロで走っている車が、1時間後にちょうど100キロ先にいるとは限らない」のと同じこと。

回帰年

上記の角度は、黄道の星座を基準にしたものだった。原点に当たる「0°」の場所は、大ざっぱに言うと、2000年の春分に、地球から見て太陽があった方向。正確に言うと「2000年1月1日12時の平均春分点の方向」だが…。

なぜ「2000年1月1日12時の」春分点…などと時刻を指定するのだろうか？

実は「春分点」自体は、星座の中に固定されたものではなく、黄道の星座の中を逆行している（逆行＝黄道上の太陽の動きとは逆向きに動く）。2000年には「うお座」にあって、今も「うお座」にあるのだが、「みずがめ座」の方向に少しずつ移動している。…「春分点が動く」というのは、ある意味、地球人の誇大妄想（？）で、「地球の自転軸の向きが少しずつ変わる」という地球ローカルの現象を、宇宙スケールの天球に投影して見ている。足元がフラフラしている酔っぱらいが「宇宙が揺れていて歩きにくいなあ」と主張するようなもの。

例えば「2001年の春分点を基準にする」「2002年の春分点を基準にする」…「今の春分点を基準にする」というのは、全部原点の位置が異なり、星座の同じ場所でも違う座標値になってしまう。そんないいかげんな座標系は使いたくないのだが、地球の四季の周期を決定しているのは、この動く「春分点」。仕方がないので、動いた量を補正しながら計算しよう。

地軸の向きが変わるこの現象は、歳差と呼ばれる（正確に言うと、地軸の向きの変動のメインの成分に当たる）。「歳差現象の結果として、春分点が黄道上を移動した角度」も歳差と呼ばれることがある。「2000年1月1日12時から通算して、春分点が黄道上を動いた合計量」を推定したのが、この式（出典 [2]）。

引き続き t の単位は1000年（＝365250日）、結果の単位は角度の「°」。

この記事の計算システムでは、「星座を基準」は、2000年1月1日12時の「春分点を基準」と同じ。それ以降については、時間が進めば進むほど、「星座を基準」と「春分点を基準」の意味がずれる。2000年1月1日12時より前についても同様で、時間が戻れば戻るほど、逆方向にずれる。prec(t) は、その「ずれ」の大きさを表す。「星座を基準とした位置（角度）」に prec(t) を足せば、「春分点を基準とした位置」が得られる。時刻 t において、春分点を基準とする地球の位置は:

前半同様、prec(t) を t について微分して365250で割ると:

これは時刻 t における「星座に対する、春分点の移動速度」（単位: °/日）を表している。春分点は「星座の中の太陽」とは逆向きに動くので、春分点を基準とすると、太陽は実際以上に大きく動いたように見える。原点がバックするので…。バックといっても円周上のことなので、逆方向から近づいていることになり、太陽と春分点が再び出会う日が早まる。春分点を基準とした太陽の位置の変化は、1日当たり:

これは (★) を時間で微分したもの。t = 0 のとき、計算済みの earth(0) と組み合わせると:

  1. f′(0) = rate_earth(0) + rate_prec(0)
  2. = 0.985 609 113 114 525… + 50287.961 95 / (3600 × 365250) = 0.985 647 357 819 028…

この場合、春分点を基準にして、太陽が 360° 動くのに必要な時間は:

これが、2000年1月1日12時における回帰年の長さ、つまり四季の一巡の平均周期。前半の D（星座に対して地球が1周する所要日数）と比べると、20分ほど短い。§1の x の大部分（ざっと8割）は、ここから生じる。

Ex2.2
t = 1 として同様の計算をすると、結果は365日5時間48分39.99…秒。これは、1000年後（西暦3000年）の回帰年の長さ。

現在、回帰年の長さは、100年につき約0.53秒の割合でだんだん短くなっている。これは春分点移動の（見かけの）加速に起因する変化で、地球の自転の遅れに起因する変化を含まない（それも含めると、回帰年はさらに短くなると見込まれる）。

§3　春分の間隔　（楕円軌道の素朴な計算）

春分の間隔を求める常識的な方法は、2回の春分の時刻を調べて、引き算することだろう。意外なことに、この常識的発想は、数学的には良くない！

でも、まずは常識的にやってみましょう…。

春分というのは、地球から見た太陽が、春分点を基準として 0° のポイントを通過する時刻。円軌道上で言えば、§2 の

が 180° になる瞬間に当たる（もちろん 360° の整数倍の差があってもいい）。この「180°」は太陽から見た地球の黄経【※】で、そのとき、地球から見た太陽黄経は 0°。同様に、地球から見た太陽黄経が 90°、180°、270° になる瞬間が、それぞれ夏至・秋分・冬至。

【※】 黄経（こうけい）は、§2で考えた「0°（原点）～360°」のこと。太陽から見た地球と、地球から見た太陽は向きが正反対なので、この値が180°ずれる。

しかし §2 で見たように、円軌道上の（理論上の）地球は「常に毎日1°」のような単純な動きをする。つまり、地球の公転角速度に季節差がないことになり、「春分の間隔」「夏至の間隔」などは、どれも回帰年とほとんど同じになってしまう。そこで以下では、楕円軌道上を公転する（リアルに近い）地球を考えてみよう。そうすれば、現実の春分・夏至などに近い値を求めることができる。

楕円軌道上の位置の表し方

PNG画像: 近日点W、楕円軌道上の天体P、円軌道上の天体Q

「地球が太陽の周りを楕円運動している」と考えた場合、太陽は（楕円の中心 O ではなく）楕円の焦点 S に位置する。軌道の楕円の長軸の両端の2点のうち、太陽に近い側を近日点（きんじつてん）または近点という（図の W）。画像は、3次元的な図解ではなく、軌道面の上の平面的なもの（楕円が大げさに描かれている。現実の軌道は、もっと円に近い）。

「近日点をスタートラインとして地球がどれだけ動いたか」を表す角度を「近点離角」という。楕円軌道上の位置 P に対応する「近点離角」、つまり ∠WSP を真近点離角と呼び、v（ブイ）で表す。これは「ケプラー運動する惑星が実際にその位置にある」という角度。

楕円軌道上の惑星が P まで動いたとき、それに対応する「円軌道上の（理論上の）惑星の位置」が Q だったとしよう（円軌道を考える場合、太陽の位置は、円の中心 O）。Q に対応する「近点離角」、つまり ∠WOQ を平均近点離角といい、M で表す。Q は円軌道上をほぼ等速運動する。

図の惑星 P は、近点から反時計回りにざっと 60° くらい動いている（角 v）。ところが、Q は、その半分くらいしか動いていない（角 M）。これは、近日点の近くでは太陽の重力が強く働き、現実の惑星は速く動くため。円軌道上の Q の動き方は、公転速度の季節差を平均化した理論上のものなので、実際の動きが速い場所では、相対的に遅く見える。角 v と角 M の関係を解明できれば、現実と理論を結び付けられる。

PNG画像: 円軌道の中心O、春分点の方向Y、近日点W、惑星の位置Q

説明のための例として、次の画像の Y が春分点の方向だとして、近日点 W が黄経 100° だとする（k = ∠YOW = 100°）。円軌道上の（理論上の）地球 Q の位置を表す M = ∠WOQ が 30° だとすると、地球 Q の黄経 λ は、λ = k + M = 100° + 30° = 130°。逆に、Q の黄経 λ が 130°、W の黄経 k が 100° とするなら、地球の平均近点離角は:

M を v に置き換えれば、楕円軌道上の真近点離角についても、同様の計算が成り立つ。

「地球軌道の近日点 W が Y から測って何度なのか」ということ（k = ∠YOW）を近日点黄経という（※注1）。これは「地球の軌道面を取り囲む宇宙の星々の中で、軌道の楕円の長軸が今どっちを向いているのか」を表す。星座を基準にした場合（つまり Y が2000年1月1日12時の春分点の方向に固定されている場合）、その具体的な値は（出典 [1]）:

近日点の黄経と近点離角は、上述のように、単純な足し算・引き算の関係にある。地球の平均近点離角 を M とすると:

Ex3.1
Ex2.1 を参考にすると、2000年3月20日12時おける地球の平均近点離角は:
  • M = earth(79/36525) − perih(79/36525) = 75°.3915…

※注1: 本当は、計算の基準になる平面と、惑星の軌道面は、同一とは限らない（＝軌道傾斜角が0°とは限らない）。それに関連して、k = perih(t) の幾何学的意味は、厳密には、上記説明と少し異なる（§4参照）。

真近点離角の求め方（近似式）

円と楕円の関係は、楕円の形（円に近いのか・細長いのか）によって異なる。楕円の形の違いを表すため、離心率 e という値が使われる。地球軌道の離心率は、時代によって少しずつ変わる。その値は（出典 [1]）:

v と M の関係は、ケプラー方程式というもので表される（§4）。ところが、M を与えて v を求める方向では、この方程式は（普通の意味では）数学的に解けない。仕方がないから、とりあえず v = M + C と置いて、ある方法で C の近似値を求めると…

ここで π は円周率。B は:

C は中心差と呼ばれる。

【※】 sin（サイン）は、角度を長さに変換する関数の一つで、電卓などを使って簡単に計算できる。B と C は、単位が違うだけで、同じ角度を表している。B の単位は「ラジアン」、C の単位は「°」。1ラジアンは (180/π)° に当たる。

Ex3.2
Ex3.1 に引き続き、2000年3月20日12時の、地球の真近点離角 v を求めよう。
  1. e = ecc(79/36525) = 0.016 708…
  2. M = 75°.3915…
  3. sin M = 0.967 671…
  4. 2M = 150°.7830…
  5. sin 2M = 0.488 117…
これらの値を使うと、中心差の近似値は:
  1. B = 2e sin M + (5e2/4) sin 2M = 0.032 507…
  2. C = B × (180/π) = 1°.8625…
従って、真近点離角の近似値は:
  • v = M + C = 77°.2540…
太陽位置の簡易計算

引き続き t = 79/365250 の例を考えよう。近日点黄経と真近点離角を足せば、楕円軌道上の地球の黄経が得られるのだった。星座を基準にすると:

そこに春分点の移動量（座標系のずれ）を足せば、春分点を基準とする地球の黄経が得られる:

代わりに、円軌道上の地球の黄経（星座が基準）に「春分点の移動量」と「中心差」を直接足し算してもいい:

太陽から見て、楕円軌道上の地球の黄経（春分点が基準）が 180°.1951… ということは、その地球から見た太陽の黄経は 0°.1951…。「春分点を基準とする太陽黄経」が 0° になる瞬間が春分なので、上記の時刻（2000年3月20日12時）は、春分の少し後に当たる（☆）。

この簡易計算法を要約すると:

  1. §2の公式を使って、与えられた時刻 t における earth(t) と prec(t) を求める。
  2. その時刻における離心率 e = ecc(t) と平均近点離角 M = earth(t) − perih(t) を求め、それらを使って中心差 C の近似値を計算する。
  3. その時刻において、ケプラー運動する地球から見た太陽の黄経（基準は春分点）を ☉ とすると、☉ = (earth(t) + 180°) + prec(t) + C である。
Ex3.3
Ex3.2 は、2000年3月20日12時の話だった。その半日前、2000年3月20日0時の太陽位置を簡易計算してみよう。半日前なので、
  • t = 78.5 / 365250
であり、上記の手順に従うと:
  1. earth(t) = 177°.8367…
  2. prec(t) = 0°.0030…
  3. e = ecc(t) = 0.016 708…
  4. M = earth(t) − perih(t) = 74°.8987…
  5. C = 1°.8585…
  6. ☉ = (earth(t) + 180°) + prec(t) + C = 359°.6983…

2000年3月20日0時は、太陽黄経が 0° に戻る少し前であることが分かった（☆☆）。

逆の問題を考えてみよう。

Ex3.4
地球がケプラー運動すると仮定した場合、2000年の春分は いつだろうか。（☆☆）と（☆）から、3月20日の午前中に太陽黄経が 0° になることは容易に想像できる。いつ 0° になるのか突き止めるため、1/24 日ごと（要するに1時間ごと）の太陽黄経を計算してみよう。
  • d = 78.5 + n/24 [n = 0, 1, 2, …, 12]
  • t = d / 365250
と置いて、3月20日0時～12時の1時間ごとの太陽位置を簡易計算すると、7時と8時の間で 0° になることが分かる。そこで今度は
  • d = 78.5 + (7 + n/60) / 24 [n = 0, 1, 2, …, 59]
  • t = d / 365250
と置いて、7時台の太陽位置を1分ごとに簡易計算すると、7時17分と7時18分の間で 0° になることが分かる。…このような手順を繰り返せば「2000年3月20日7時17分10秒くらいが春分」という答えが出る。手順は原始的だが、発想は単純明快。ただし、この計算結果は、ケプラー方程式による厳密解（同7時17分34秒）と比べて24秒もずれている！

t = 0.1 が2100年1月1日12時、t = −0.1 が1899年12月31日12時であることを参考に、上記同様、1900年と2100年の春分の時刻を推定すると、それぞれ「1900年3月21日01時35分54秒」「2100年3月20日12時59分55秒」。両者の時間差「73048日11時間24分01秒」（2100年には2月29日がないことに注意）は、春分間隔200回分なので、200で割ると…。「1900～2100年の平均春分間隔は、365日5時間49分01.205秒」となる。この結果は、理論上の真の値（365日5時間49分01.186秒）にかなり近く、誤差0.1秒未満。秒の端数を四捨五入して、とりあえず「西暦2000年ごろの平均春分間隔は365日5時間49分01秒」と結論できる。

いいかげんな計算だったくせに、どうしてこんなに正確なのだろう？

タイミングの誤差が最悪50秒だとしても、2個のタイミングの間隔の誤差は最悪で100秒。2個のタイミングが200年の間隔なら、それを平均したとき、真の平均値と比べて、1年当たりの誤差は最悪でも0.5秒で済む。この方法なら、誤差の大きいデータからでも、意外と精度の高い答えを出せる【※】。これに加えて、中心差の誤差は周期的なので、「2回の春分」のような場合、2個のデータが同様の誤差を持ち、引き算すると誤差が消える傾向がある。

【※】 「そういう理屈なら、200年ではなく、1000年や1万年の平均を求めたら、さらに正確になるのでは？」とも思える。残念ながら、長期間のうちには春分間隔そのものが、複雑に大きく変動する。従って、あまりに長い期間の平均だと「考えている時点での春分間隔」とずれる可能性が高い。「200年の平均」にしても、「その期間なら良い」という保証があるわけではないが、とりあえず実用になる。

代わりに「1900年と2100年において、太陽黄経90°になる時刻」を求めると、夏至の平均間隔が分かる。秋分間隔・冬至間隔も同様。

Ex3.4 の「1時間ごとの位置・1分ごとの位置…を求める」というのは、計算量的に無駄が多いが、それは枝葉の問題。この計算法には、実はもっと根本的な「欠陥」がある――結果の数値は実用になるものの、概念的な不整合をはらんでいる。次の §4 で、より良い解法を紹介したい。「不整合」の性質も、そこで明らかになるだろう。

計算上の時刻 vs. リアルな時刻

現実の地球の位置は、月や惑星の作用によって多少乱されるし、地球上で現実の春分を観測する場合、その他にも微妙な「ずれ」が生じる。上記で計算したのは、現実の春分の時刻ではなく、純粋なケプラー運動を仮定した場合の（乱れ・ずれが一切なかったら…という理想の状況における）理論上の春分の近似値。そのため、理科年表のようなものに記載される「現実の春分」の時刻と、完全に一致するとは限らない。

でも「理論上の計算で、現実と合わないから駄目」とは限らない。月・惑星などの影響による「乱れ」は、長期の平均を考えればほぼプラスマイナス・ゼロになるため、純粋なケプラー運動を仮定した計算結果は「現実の春分の平均間隔」と、非常によく一致する。夏至・秋分・冬至についても事情は同じ。現代を中心にする約1万年間の範囲に関して、この節の方法で春分などの平均間隔を計算した場合、誤差は1秒未満。となれば、最初から乱れがゼロとして計算した方がスマートだろう。同じ答えが、簡単な計算で得られるのだから…。

注意点として、計算の基礎となっているいろいろな式の時間変数 t は、「太陽系力学時」（TDB）という時刻系による。このシステムは日本時間マイナス9時間とだいたい同じだが、うるう秒がない。つまり、1日は常に 24 × 3600 = 86400 秒。日常の時刻は、ときどき「うるう秒」が入って1日が 86401 秒になるので（地球の自転の遅れを反映）、1日が平均的に少し長い。

残念ながら、今後うるう秒がどのくらい入るかは誰にも分からないため、日常の時刻を基準とする正確な計算は、誰にもできない。

この記事では、過去や未来の地球の自転について当てずっぽうの予想をせず、1日の長さは一定とする。この仮定のため、計算結果は、リアルな時刻と微妙に異なる。「数千年後の回帰年・春分間隔・春分年などの日数の誤差が1秒未満」といっても、それは「基準の時刻系における」誤差。リアルな時計（協定世界時や日本標準時）で測った時間とは、最大10秒くらいずれるかもしれない。

半面、この方法なら、「地球の自転の乱れ」に左右されない「正味の時間」を知ることができる。「回帰年の長さと比べて、現在、春分間隔は長く夏至間隔は短い」といった相対的な比較については、上記の問題にかかわらず、リアルな結論が得られる。

メモ

1. 真近点離角・平均近点離角は、それぞれ真近点角・平均近点角と呼ばれることも多い。

2. 楕円の離心率は e で表されることが多い。「自然対数の底」も e で表されるが、それとは全く意味が異なる。近日点黄経は、本来 ϖ（ギリシャ文字の「パイ」の筆記体）で表される。分かりやすさのため、ここでは代わりに k とした。

3. この節の計算では、特定の年の（ケプラー的な）春分の時刻を正確に決定できなかった（誤差が、ざっと30秒もあった）。時刻計算を正確にしたい場合、最も確実な方法は、ケプラー方程式を数値的に解くこと。簡易計算としては、中心差の近似式を次のものに置き換えるだけで、誤差0.1秒未満になる。

精度が向上するのは、あくまで「ケプラー的な春分・夏至」などの時刻。現実の春分・夏至などの正確な時刻を知るには、章動・摂動などを補正する必要がある（平均間隔の計算では、現実の春分・夏至などにこだわるのは、かえって良くない）。

4. 歳差の式 prec(t) と rate_prec(t) は、他の式同様 t = −6（紀元前4000年）から実用になる。未来については有効範囲が狭く、t = +3（紀元後5000年）を過ぎると急激に精度が悪化する。1秒の精度で十分なら西暦2000年±6000年の1万2000年間で使えるが、0.1秒の精度が欲しい場合、紀元前3500年から紀元後5000年までの8500年間の計算が可能。

§4　解析的な解法　（瞬時の春分年）

ここからが本題。平均速度に基づく近似、瞬間速度に基づく厳密解の順で話を進める。

平均速度に基づく春分年の近似

§3 では2回の春分のタイミングを求め、引き算により春分間隔を決定した。アプローチは単純だが、実際の計算は面倒だった。よく考えれば、これはケプラー運動する地球が「軌道上のそこで春分になる点」を通過する周期の問題で、地球と「その点」の相対角速度を求めて、360° をその角速度で割った方が、見通しが良いだろう。そのように定義された周期を春分年と呼ぶことにする。夏至年・秋分年・冬至年も同様に定義される。

1. 円軌道上を反時計回りに平均運動する地球を考える。その地球の平均近点離角は（§3）:

直観的には、これを時間で微分すると、近日点に対する地球の角速度 ρ が得られる。単位を「°/日」として:

ここで rate_earth(t) は earth(t) の変化率。§2 で計算済み。rate_perih(t) は、恒星を基準とした近日点の位置 perih(t) （§3参照）を時間で微分して 365250 で割ったもの。同じ単位で測った近日点移動の速度を表す。

近日点黄経 perih(t) は本当は角度ではないが、一定条件下で、角度のように振る舞う（12.）。

以下では、一貫して（瞬時の平均）春分点を基準にするが、計算上 ρ をそのまま使うことができる。なぜなら、春分点を基準にする場合、地球黄経 earth(t) と近日点黄経 perih(t) の両方に（春分点移動に関する）座標系の補正 prec(t) が加わるのだから、それらの差である離角は最初と変わらず、従って、その離角の変化率も最初と変わらない。

PNG画像 (12 KiB)

2. 楕円軌道上において、「もし今が春分だったら地球はそこにある」という点 P を考える。P から見た太陽 S は、春分点を基準に黄経 φ = 0° であり、従って点 P の真近点離角は:

第1項は、（地球軌道上の）点 P から見た太陽の黄経 φ を、太陽から見た P の黄経 Φ に直したもの。Φ = φ + 180° である（または、同じことだが Φ = φ − 180°）。第2項は、春分点を基準にした近日点黄経（つまり perih(t) プラス「春分点の移動量」）。

近日点黄経 perih(t) も、春分点の移動量 prec(t) も、時間とともに増大する。だから v は時間とともに減少する。つまり点 P は地球とは逆回りに動き、負の「角速度」を持つ。このセクションの P は単なる軌道上の点で、地球の位置ではない。例えば2000年1月1日12時は春分ではないので、地球はその点にはないが、「その点」自体は、軌道上のどこかに、常に存在している。

この点は、太陽から見ると（春分点ではなく）秋分点方向に当たる。

3. P に対応する円軌道上の点を Q とする。P は天体ではなく、ケプラー運動するわけではない。典型的なケプラー方程式の問題では、（円軌道上の天体の）平均近点離角 M の変化率がほぼ一定で、v の変化率が場所によって異なる。ここでは逆に、（楕円軌道上の点 P の）真近点離角 v の変化率がほぼ一定で、M の変化率が場所によって異なる。M の瞬間変化率 μ は、直観的には「近日点 W に対する Q の角速度」（ここでは単位を「°/日」とする）。この値が分かれば、春分年の長さは 360° / (ρ − μ) として容易に決定される。

（解説）

  1. 地球については円軌道上で考える。地球の瞬間角速度を楕円軌道上で考えるのは面倒だし、そんなことをすれば、地球の角速度（ひいては春分年などの長さ）に1年周期の変動が生じてしまう。円軌道上の地球の（W に対する）瞬間角速度 ρ なら簡単に計算可能で、それを使えば、（本質と関係ない）1年周期の変動は自然に除去される。地球の具体的位置は不要。円軌道上のどこかをこの角速度で動いている、ということだけ分かればいい。
  2. 「その点」については、まずは楕円軌道上で考える。「春分年・夏至年などの違い」はここから生じる。実際の計算では P を Q に変換、円軌道上の問題に帰着させる。Q は、円軌道上において「もし今が春分だったら地球はそこにある」という点。P の離角 v は既知。それを Q の離角 M に変換することは容易（4. 参照）。
  3. M の瞬間変化率（つまり Q の W に対する瞬間角速度）を μ と置く。地球と逆回りに動く Q は、円周の反対側から地球に接近している。従って、地球と Q の相対角速度は ρ − μ。それに基づくと、地球と Q が出会う周期は 360° / (ρ − μ) 日。

4. P を通り楕円の長軸と直交する直線 ℓ を考える。楕円に外接する円と ℓ の交点のうち、長軸から見て P と同じ側にあるものを R とする。E = ∠WOR を離心近点離角と呼ぶ（O は円の中心）。ケプラー方程式に関する基本事項として、角度の単位をラジアンとするとき、平均近点離角 M と真近点離角 v について、E を仲立ちとする次の関係が成り立つ（詳細は別記事）。

ここで e = ecc(t) は楕円の離心率。普通ケプラー方程式というと、M から v を求める問題を指し、それは難しい。しかし、逆に v から M を求めるのは易しい。(4.1) から:

(4.2) から:

  1. tan (E/2) = tan (v/2) / √(1 + e) / (1 − e)
  2. tan (E/2) = tan (v/2) × √(1 − e) / (1 + e)

u = tan (v/2), w = (1 − e) / (1 + e) と置くと:

  1. tan (E/2) = u√w ……… (4.4)
  2. E/2 = arctan (u√w)
  3. E = 2 arctan (u√w) ……… (4.5)

−π ≤ v ≤ π とすると、v → ±π (= ±180°) のときには u = tan (v/2) → ±∞ だが、その場合、数値的には E → ±180°。幾何学的にも、極限では点 P が遠日点 A に一致し、点 R も A に一致するのだから E = ±180°。そのとき M = E − e sin E = ±180° となる（sin E = 0 なので）。

5. 以上によって、実用上、問題は既に解決している。

例えば、t = t1 = 0（2000年1月1日12時）とすると:

  1. t1 における P の真近点離角（2. より）: v1 = (0° + 180°) − (prec(t1) + perih(t1)) = 77°.062 651 92 = 1.344 997 006 322 329 860… [radians]
  2. それに対応する Q の平均近点離角: M1 = 75°.201 891 032 205 203…

…離角の変換の方法は、4. の通り: その瞬間の離心率 e1 = ecc(t1) = 0.016 708 634 2 から

を求めて、

と組み合わせて (4.5), (4.3) から

  1. E1 = 2 arctan (u1√w1) = 1.328 742 140 323 818 993… [radians]
  2. M1 = E1 − e1 sin E1 = 1.312 520 602 237 977 889… [radians] = 75°.201 891 032 205 203…

とすればいい…。

ちょうど1日後の t2 = 1/365250 において、同様の計算を行うと:

  1. t2 における真近点離角: v2 = 180° − (prec(t2) + perih(t2)) = 77°.062 604 843 936 954…
  2. 対応する平均近点離角: M2 = 75°.201 844 453 748 219…

この1日での Q の平均近点離角の変化は:

とりあえず ΔM を1日当たりの角速度だとしよう（それは瞬間速度ではなく、24時間の平均速度だが…）。ところで、t1 における、地球の平均近点離角の変化率（瞬間角速度）は:

ρ − ΔM を、2000年1月1日12時における地球と点 Q の相対角速度の近似値とすることができる（単位は「°/日」）。このとき、春分年の長さの近似値は:

この結果は（計算として）小数第9位まで正しく、計算誤差は8600ナノ秒程度（ナノ秒＝10億分の1秒）。

「φ = 0°, Φ = φ + 180°」の φ = 0° を φ = 90°, 180°, 270° に置き換えれば、それぞれ夏至年・秋分年・冬至年の計算となる。

6. 上記は単純化した例。実際には、与えられた時刻 t を中心に、挟み込むように、ちょうど半日前からちょうど半日後までの「1日」の ΔM を考えた方が計算精度が高くなる。「1日」の代わりに「1秒」や「1ミリ秒」（1000分の1秒）、いっそのこと「1ナノ秒」当たりの「角度の変化」を調べて、単位が「°/日」になるように定数倍すれば、事実上、瞬間速度を求めたのと変わらない。このような「1日」とか「1秒」とかの計算区間設定を刻み幅と呼ぶ。

刻み幅をミリ秒やナノ秒にすることで瞬間速度に近づくことができるのは、十分な計算精度（例えば有効数字50桁）がある場合に限られる。倍精度演算（コンピューター上の一般的な計算。有効数字15～16桁）の場合、40日前から40日後までの「80日」くらいの刻み幅が良いようだ。下記 7. の方法を使うと、さらに良い。

瞬時の春分年

7. 上記の ΔM の代わりに M の瞬間変化率 μ を使えば、正確な瞬時の解析解が得られる。ラジアンを単位とした角度 M = E − e sin E を、時間で微分すると:

  1. μ = M′ = E′ − [e′ sin E + e(cos E)E′]
  2. μ = E′(1 − e cos E) − e′ sin E ……… (1)

ここで

は、§3 の ecc(t) を時間で微分して 365250 で割ったもの（1日当たりの変化量）。後は E′ の値さえ分かれば、式 (1) の右辺を計算できる。E = 2 arctan (u√w) なので:

  1. E′ = 2 × {1 / [1 + (u√w)2]} × (u√w)′
  2. = 2(u√w)′ / (1 + u2w) ……… (2)

(u√w)′ の値は次の通り。u = tan (v/2) なので

  1. u′ = [1 / cos2 (v/2)] × (v/2)′
  2. = [1 + tan2 (v/2)] × v′/2
  3. = (1 + u)2v′/2 ……… (3)

であり、(√w)′ = w′/(2√w) なので:

  1. (u√w)′ = (1 + u)2v′/2 × √w + u × w′/(2√w)
  2. = (1 + u)2v′√w / 2 + (uw′/√w) / 2
  3. = [(1 + u2)v′√w + uw′/√w] / 2 ……… (4)

(4) を (2) に代入すると:

  1. E′ = {2[(1 + u2)v′√w + uw′/√w] / 2} / (1 + u2w)
  2. = [(1 + u2)v′√w + uw′/√w] / (1 + u2w) ……… (5)

実用上、ほとんどの場合 (5) から E′ を計算できる。けれど、例えば v = π (= 180°) の場合 u は定義されず、(5) は計算不能。(5) の右辺に w = (1 − e) / (1 + e) と w′ = −2e′/(1 + e)2 を代入し、分子・分母を (1 + e) / (1 + u2) 倍すると、次の簡潔な形が導かれ（詳細は付録B）、任意の v に対して E′ を計算できるようになる。

ここで j = √1 − e2 である。v の定義から、「°/日」を単位とする v′ は次の通り。

式 (1)～(6) の μ, E′ の単位は「ラジアン/日」。v′ の単位を「ラジアン/日」に変換して (6) から E′ を求め、それを使って (1) を計算し、結果を「°/日」単位に再変換すればいい。

8. 時刻 t = 0 に対する春分年（φ = 0°）の計算例:

  1. v = (φ + 180°) − (prec(t) + perih(t)) = 77°.062 651 92
  2. e = ecc(t) = 0.016 708 634 2
  3. e′ = rate_ecc(t) = −1.150 897 741 273 100 616 016… × 10−9
  4. w = (1 − e) / (1 + e) = 0.967 131 912 451 698 150 435…
  5. u = tan (v/2) = 0.796 328 979 416 363 861 888…
  6. E = 2 arctan (u√w) = 1.328 742 140 323 818 993 591… [radians]

式 (6) と式 (1) を使うと:

  1. j = √1 − e2 = 0.999 860 401 027 648 754 418…
  2. v′ = −rate_prec(t) − rate_perih(t) = −4.707 606 270 438 816 640… × 10−5 [°/日] = −8.216 322 930 668 795 438 911… × 10−7 [radians/日]
  3. E′ = −8.173 382 392 119 454 813 360… × 10−7 [radians/日]
  4. μ = −8.129 474 395 384 374 034 353… × 10−7 [radians/日] = −4°.657 845 725 151 913 050… × 10−5 [°/日]

一方:

ゆえに、この瞬間の春分年の日数は:

この結果は [1] と [2] を前提とするもので、計算上は全部の桁が正しい。「その前提が、現実とどこまで一致するか」「実際の有効数字は何桁か」は別問題。

検証用として、過剰な精度で数値を掲載した。コンピューター上の倍精度演算では、通常どの数値も左から15～16個目の数字（整数部分が3桁なら小数12～13桁目）に誤差が入り、上記と完全には一致しない。付録Aとして、JavaScript による実装例を収録。

9.（別解）　ρ と v′ は、どちらも rate_perih(t) を含んでいる。これを「くくり出す」ことで、次の表現が得られる。

  1. C0 = 1 − e cos E
  2. C1 = jC0
  3. C2 = 1 + e cos v
  4. C3 = [(C0/j) sin v + C2 sin E] × 180/π

と置き、

とすると、求める日数は:

詳細は付録C。実装例は付録A3。

10. 表2は、グレゴリオ暦・紀元前1000年～紀元後5000年の計6000年について、毎年1月1日0時（太陽系力学時）における春分年・夏至年・秋分年・冬至年の日数を計算し、計算誤差の絶対値を調べたもの。倍精度Aは 5. の計算法（刻み幅80日）、倍精度Bは 7. の計算法、倍精度Cは 9. の計算法。誤差の単位を U = 2−44 日（およそ5ナノ秒）とした。倍精度の浮動小数において、U は 365 付近の数の ULP に当たる（＝原理的に、0.5U までの誤差は避けられない）。

表2: 春分年・夏至年・秋分年・冬至年の日数の計算精度
各6000サンプルの誤差分布
誤差の絶対値 倍精度A 倍精度B 倍精度C 4倍精度
0.5 U 以内 3755 12565 13034 24000
0.5–1 U 3719 7920 7976 0
1–2 U 7016 3476 2983 0
2–4 U 7414 39 7 0
4–8 U 2031 0 0 0
8–16 U 65 0 0 0
16–32 U 0 0 0 0

倍精度A・倍精度B/Cは、いずれも誤差100ナノ秒未満。それぞれ日数の小数11桁程度・12桁程度の精度を持つ。これは理論と現実の違いの評価ではなく、浮動小数点演算の誤差（丸め誤差など）の評価。

追記: 倍精度演算の場合、rate_earth(t) を次のように実装すると、計算誤差が減る（他の関数についても、同様の工夫が可能。実装例 v2 参照）。

表2.1は、この1カ所の改善の効果を示す。

表2.1: 演算の順序などを工夫した場合の例
誤差の絶対値 倍精度C（旧） 倍精度C（新）
0.5 U 以内 13034 16749
0.5–1 U 7976 6518
1–2 U 2983 733
2–4 U 7 0
春分間隔と春分年

11. 「春分年」は「瞬間変化率に基づく値」だが、「春分間隔」は 6. の「平均変化率に基づく近似値」の一種。正確に言うと、春分間隔というのは、

に当たる。

その気になれば、この刻み幅を使って「春分とは限らない任意の時刻から測った仮想春分間隔」を定義することもできる。もちろんそれは、刻み幅 → 0 の「瞬時の春分年」とは異なる。

春分間隔を「春分年」と定義することは不可能ではないが、実際上、いろいろとまずい事態を招く…。

（ア） 表面的な問題として、瞬時の年の近似としては、約365日という刻み幅は長過ぎる。

（イ） もっと本質的な問題として、「春分間隔を刻み幅にする」ということは、計算で求めようとしている春分間隔を計算の入り口で使うことであり、堂々巡りを発生させる。つまり「春分間隔」を決定するためには、その期間における歳差・近日点移動の量が必要だが、その量を知るには「その期間」の長さが必要。ところが「その期間」というのは、今から求めようとしている「春分間隔」…。逐次近似（回帰年の長さを第0近似とする）によってこの「もつれ」は克服可能だが、用語の定義の中にその用語が含まれているようなもので、良い状況ではない。

（ウ） さらに深刻な問題として、現代天文学における「回帰年」は瞬間変化率に基づいて定義されている。「瞬間変化率に基づく春分年・夏至年・秋分年・冬至年など」はそれと整合的だが、「春分年とは現実の春分の平均間隔」などと定義してしまうと、「回帰年」の概念と不整合になる。

例えば、「回帰年は、春分年・夏至年・秋分年・冬至年の平均に等しいか？」という考察（あるいは数値実験）には、精密な概念構成（そしてミリ秒より細かい計算精度）が必要。「春分から測った春分間隔」「夏至から測った夏至間隔」「秋分から測った秋分間隔」「冬至から測った冬至間隔」という4つの日数を足して4で割り、「年初の回帰年」と比較するのは、意味がない（5回の測定はどれもタイミングが異なる上、間隔と瞬時の年は直接比較できない）。ちなみに、「回帰年」と「春分年・夏至年・秋分年・冬至年の平均」は、一般には等しくない。

瞬時の年も仮想間隔も、任意の瞬間に対して定義される。例えば、ある年の「1月23日4時56分における春分年」や「冬至における夏至年」を考えることもできる。一方、「春分年（または春分間隔）は、春分を基準に測定するもの」という常識にとらわれると「いつが春分なのか」を決定する手間が発生し、その場合、ケプラー方程式の複雑さから逃れられない。「春分・夏至などの時刻と無関係に、春分年・夏至年などを考える」という発想の転換によって、概念が整理され、問題が（「平均→真」ではなく逆向きの）「真→平均」離角の変換に帰着される。この逆転は「解析的に解けるか、数値解しかないか」の大きな違いとなる。これが第1の分かれ道。

「瞬時の春分年」「仮想春分間隔」のどちらについても、次の事実が成り立つ。

ところが「春分間隔そのものを刻み幅にする」という設定だと、（イ）の循環論法が絡むため、別の場所で逐次近似が必要になる。これが第2の分かれ道。

「逐次近似が必要になる」ということは、数学的な間違いではなく、計算量的な問題にすぎない。どちらの分かれ道をどちらに曲がったとしても、しっかり数値計算をすれば、結果は（選択された定義に照らして）いくらでも正確になる。とはいえ、どの道が良いかは明らかだろう。

付記: 「2000年春分から2001年春分までの日数」（春分間隔）は、その期間のちょうど中央で測った春分年の日数とほぼ等しいが、「2000年春分（または2001年春分）における春分年の日数」とは等しくない。これと同様の現象は、恒星年や回帰年や近点年においても発生する。ある瞬間に、地球の黄経（恒星を基準にする）が λ だったとして、その瞬間における恒星年の長さが D 日だったとしよう。この場合、その瞬間のちょうど D 日後の地球の黄経（恒星を基準にする）が、λ + 360° になるとは限らない（平均運動の永年加速のため）。

技術的な補足

12. 近日点黄経 perih(t) は、軌道傾斜がゼロの場合には「軌道面の上で、（瞬時の、または特定時点の）春分点と近日点の成す角度」だが、軌道傾斜がゼロではない場合には「昇交点黄経 Ω と近日点引数 ω の和」と定義される。ここで、昇交点黄経は、黄道面（太陽系全体の基準面）の上で春分点と昇交点が成す角であり、近日点引数は、軌道面の上で昇交点と近日点が成す角。軌道傾斜がゼロではない場合、この二つは別の平面なので、（Ω と ω は、それぞれ普通の角度だが）Ω + ω は、普通の意味での角度ではない。

perih(t) は2000年1月1日12時の春分点を基準とする数値で、この座標系では、t ≠ 0 のとき、わずかながら軌道傾斜が存在する。従って、perih(t) を角度と呼んだり rate_perih(t) を角速度と呼んだりするのは、厳密には正しくない。その場合、天体や軌道上の点の黄経というのは、「折れ曲がって（二つの平面にまたがって）定義される近日点黄経」を経由して「春分点→昇交点→近日点→その点」と測った「角度もどき」。近日点から先の離角は普通の角度だが、全体としては「角度」ではない。

にもかかわらず、われわれの計算の背景において、Ω と ω はいずれも時間の多項式で表現され、時間で微分可能なのだから、perih(t) = Ω + ω を微分してその変化率を得ることは、数学的に可能。その結果は、「角度もどき」である黄経の変動の様子を、正確に表す。

13. t = 0 において、定義上、地球の軌道面は、基準面（黄道面）と一致する。幾何学的には、その瞬間に降交点と昇交点が180°ジャンプして入れ替わる…。しかし、そうだとすると、Ω と ω は軌道上を滑らかに動かないことになり、t = 0 における rate_perih(t) の意味について、特別な考察が必要になる。この記事の基になっている [1] では、t < 0 の地球の軌道傾斜角を負としているため、この問題は発生しない。つまり、その期間において、地球が基準面を上から下に横切る場所を「負の傾斜角で下から上に上昇する場所」と見なしている（負の傾斜角の上昇とは、要するに下降だが…）。

§5　考察

現実との比較

§4 の方法で求めた春分年・夏至年などの日数（理論値）は、現実の平均（視黄経に基づく現実の春分・夏至などについての値）と、どの程度一致するのだろうか。

「現実の平均」という場合、「平均」の定義が難しい。数年間の平均では章動・摂動の影響を除去できないが、かといって、100年間や1000年間の平均を考えると、その期間における回帰年・春分年の長さの（永年的）変動に影響されてしまうだろう。

差し当たっては、J. Meeus (1991), TABLE 26.A の近似式の、1次の係数との比較を行う。Meeus の式は、現実の春分・夏至などのタイミングになるべく合うように調整された近似多項式。その意味で「現実の春分年・春分間隔」などと結び付いている。結果は、表3-1の通り。計算の瞬間のユリウス日 JD は Meeus が採用したもの（力学時）。同一の瞬間における春分年・夏至年などの比較ではない。

表3-1: 紀元前1年ごろの春分年・夏至年・秋分年・冬至年の日数
歳差モデル 日数 計算の瞬間（JD）
春分年 La86? 365.242 137 40 J.M. 1721 139.29189
La86365.242 137 403 964
IAU06365.242 137 525 401
LT11365.242 137 656 107
夏至年 La86? 365.241 725 62 J.M. 1721 233.25401
La86365.241 725 616 450
IAU06365.241 725 727 522
LT11365.241 725 859 620
秋分年 La86? 365.242 495 58 J.M. 1721 325.70455
La86365.242 495 588 366
IAU06365.242 495 688 235
LT11365.242 495 813 821
冬至年 La86? 365.242 882 57 J.M. 1721 414.39987
La86365.242 882 573 737
IAU06365.242 882 683 641
LT11365.242 882 807 895

J.M. は Meeus の係数、それ以外は t = (JD − 2451545)/365250 として、§4 の方法で求めた日数（精度20桁以上の計算を行い、小数第13位を四捨五入）。§2 の prec() と rate_prec() は、IAU 2006 歳差モデルによる（IAU06: 文献 [2]）。1990年代の Meeus の計算では、背景にある歳差理論が異なる。その系統誤差を緩和するため、ここでは Laskar (1986) による一般歳差の10次式を使い、別計算を行った（La86: 文献 [3]）。Meeus の係数は、La86 を使った場合の値の、小数第9位を四捨五入（秋分年に関しては切り捨て）したものと一致する。参考として、Vondrák et al. (2011) の長期モデル（LT11: 文献 [4]）による計算値も、表に含めた（恐らく現実に近い）。

Meeus (1991), TABLE 26.B は、西暦2000年を中心とする同様の近似式。それとの比較は次の通り。西暦2000年前後では、IAU06 と LT11 はほとんど同じなので、この表には後者を含めていない。

表3-2: 西暦2000年ごろの春分年・夏至年・秋分年・冬至年の日数
歳差モデル 日数 計算の瞬間（JD）
春分年 La86? 365.242 374 04 J.M. 2451 623.80984
La86365.242 374 044 692
IAU06365.242 374 884 709
夏至年 La86? 365.241 626 03 J.M. 2451 716.56767
La86365.241 626 023 616
IAU06365.241 626 898 275
秋分年 La86? 365.242 017 67 J.M. 2451 810.21715
La86365.242 017 666 350
IAU06365.242 018 519 235
冬至年 La86? 365.242 740 49 J.M. 2451 900.05952
La86365.242 740 492 026
IAU06365.242 741 311 714

この場合も、Meeus の採用値は、La86 による計算値の小数第9位を四捨五入（夏至年に関しては切り上げ）したものに等しい。これを単純に解釈すると、現実と理論は高い精度で一致している（誤差1ミリ秒未満）。別の解釈として、現実に合わせた Meeus の式自体、実は本稿と似た方法で生成されているのかもしれない。「Meeus の計算と、この記事の計算が同じ系統誤差を持っていて、そのため相対的に一致している」という可能性もある（実際、どちらも VSOP87 ベース）。

表から明らかなように、歳差モデル次第で、日数の小数第6位が変わる場合がある。常識的に考えると「20年前の La86 より現代の IAU06 の方が正確（※注2）。理論の違いによる結果の違いが小数6桁目なら、20年前の理論を使った場合の有効数字は小数6桁。今の有効数字は小数7桁以上」と期待される。悲観的に見れば「まだ6桁目は怪しい」が、20年間の技術の進歩で1桁くらい精度が上がっていても、おかしくない。さらに「歳差理論の不確実性と比べ、地球の公転理論に関する不確実性は無視できるくらい小さい」と仮定すると、われわれの計算の（現実に対する）不確実性は、歳差モデルの不確実性程度。IAU06 モデルを使った場合、−3500年から+5000年までの8500年間において、春分年などの日数の誤差は0.1秒未満と予想される。その内訳は、IAU06 の多項式表現から生じる誤差（20万年間有効な LT11 との比較による分析）が0.05秒未満、そこから類推して「LT11 自身から生じる不確実性」が0.05秒未満。

※注2　実際には、時間変数の範囲によっては、必ずしもそうとは言えない。やや遠い未来になると、La86 の方がむしろ良いようだ。しかし、現代の値に関しては、現代の理論の方が正確だろう。

それとは別に、地球の軌道要素について多項式近似を使った場合、遠い過去・未来において信頼度が低下する。1990年代の [1] と IAU 2006歳差モデル [2] には、計算の枠組みとなる座標系の違いもある。今の枠組みから派生的に再定義された「J2000座標系」と、前世紀に「J2000座標系」とされていたものは、完全に同じでは ないのかもしれない。そのため、本稿のように [1] と [2] を単純に組み合わせて計算するのは、厳密には正しくないのかもしれない。もっとも、黄経の座標値が一律に50ミリ秒角くらい変わるとしても、それだけでは計算結果に有意な影響は生じないようだ。

LT11 の時刻系はTDBではなくTTで、計算の基準時刻が最大2ミリ秒ほどずれるが、その影響はモデル自体の不確実性に比べて十分に小さいだろう。

暫定的な結論として（大ざっぱな目安だが）、現代の理論水準においては、春分年などの日数は小数第6位（0.1秒）まではほぼ正確で、順当にいけば7桁目（0.01秒）もまあまあ信用できる。8桁目（1ミリ秒）以下は、多少現実とずれるだろう。この理論精度を生かすには小数9桁以上の計算精度が必要だが、倍精度B/Cは、そのニーズを手軽に満たしてくれる。公転理論の古さ（精度）は問題だが、この記事のテーマは「計算アルゴリズムの提示」。rate_earth() その他の関数の中身は、後から自由に置換できる。

以上の議論は太陽系力学時に基づくもので、ΔT の補正は考慮されていない。

他の資料との比較

外部サイトのグラフ「Länge des tropischen Jahres nach verschiedenen Definitionen」（異なる定義による回帰年）は興味深い。中央の黒線は（通常の意味での）回帰年の長さ、それを取り囲む4種類のカーブは「春分間隔・夏至間隔・秋分間隔・冬至間隔」らしい（2006年の資料）。下記のグラフ1は、§4 の春分年・夏至年・秋分年・冬至年を同様にプロットしたもの。

PNG画像 (24 KiB)
グラフ1: 瞬時の春分年・夏至年などの日数の変動（太陽系力学時）
「2000」は2000年1月1日12時、1000年＝365250日とする。

文献 [1] による地球の軌道要素（多項式表現）の有効範囲は J2000 ±6000年だが、比較に便利なように、こちらも±8000年分を出力した。§2 の prec() と rate_prec() では有効期間が足りないため、グラフの作成では LT11 [4] を利用した（§2 の式を使った場合、「10000」年における誤差は約3秒と推定される）。

二つのグラフは、かなりよく一致しているが、微妙な違いもある。例えば、春分年の長さを表す曲線の左端の値が 365.2420 のグリッドより上か下かという違い。これは、主に歳差モデルの違いだろう。La86, IAU06, LT11 のどれかを使った場合、計算上、「−6000」年における春分年の日数は 365.24199 台だが、[1] の歳差モデルを使うと、この値が 365.24200 台になる。このグラフのスケールでは「瞬時の年」と「間隔」の違いは問題にならないが、「Länge des tropischen Jahres」の作者は、それらの概念的不整合に一定の注意を払っている。

別のサイトには、「The Lengths of the Seasons」というコンテンツがある。数値的に「季節の長さ」が考察され、そこから発展して「春分間隔」のような概念も扱われている。怪しい点もあるが、参考になる点も多い。同サイトのチャート10は、本稿の結果とほぼ一致している。数値計算なので解析的方法より扱える期間が長く、その点は うらやましい。

「「冬至年」について」（2009年の資料）では、「平均→真近点離角」という変換（本稿の §4 の逆方向）を使った簡易計算が行われている（記述には不正確な部分もある）。その計算モデルは回帰年の長さを一定とするもので、100年につき0.53秒のオーダーで、現実とずれる。これはモデル選択の問題で、アルゴリズムの本質的欠陥ではない。「平均→真」でも、計算量を気にしなければ、望むだけの精度で変換を行うことができる。その場合、中心差の近似級数またはケプラー方程式の逐次近似が必要だが、「「冬至年」について」では、何と1次の中心差が使われている。たった1次でも意外と計算誤差が出ないところが、注目に値する。

「「冬至年」について」と本稿の §4 は、アルゴリズムは異なるものの、どちらも次のアイデアに基づいている。

  1. 春分年（あるいは平均春分間隔）などの計算においては、「視黄経に基づく複雑な計算をしてから、平均処理により、章動などの影響を排除する」より、初めから章動などを無視して、単純なケプラー運動を仮定する方がいい。
  2. ということは、実際に春分のタイミングを決定して引き算で春分間隔を求める必要もなく、春分年・春分間隔そのものを解析的計算の対象とすることができる。

「回帰年と春分年」というブログ記事（2011年）には、「厳密に計算すると西暦1872年時点での春分年は365.24236日」と記されている。§4 の方法で計算すると、グレゴリオ暦1872年1月1日0時と1873年1月1日0時の春分年の日数は、それぞれ 365.242 3613 台、365.242 3614 台（古い La86 モデルでは、どちらも 365.242 360 台）。「回帰年と春分年 補足2」の計算式は、一見簡単そうだが…。よく見ると、左辺の「年」を定義する右辺に「1年分の」という言葉があり、閉形式になっていない。その根底には、「刻み幅のある計算」と「瞬時の年」の概念的不整合がある。一般向けのブログ記事なので、分かりやすさのために厳密性を犠牲にした、という面もあるのだろう。

展望

春分年などの変動に関して、回帰年との差を考えると印象深い。

表4: 回帰年との長さの差（単位:秒）。時刻系は太陽系力学時。
「2000」は2000年1月1日12時。1000年＝365250日。
年 春分年 (M) 夏至年 (J) 秋分年 (S) 冬至年 (D)
−1000 −29.88 −43.93 +30.30 +43.51
0 −14.94 −50.52 +16.02 +49.44
1000 +0.84 −52.14 +0.47 +50.82
2000 +15.93 −48.69 −14.84 +47.60
3000 +28.90 −40.64 −28.40 +40.14
4000 +38.61 −28.96 −38.82 +29.17
PNG画像 (17 KiB)
グラフ2: 表4の数値を図示したもの

西暦1000年ごろに、興味深いイベントが続発したことが見て取れる。回帰年との比較で、春分年がマイナスからプラスに、秋分年がプラスからマイナスになったこと。夏至年が極小、冬至年が極大になったこと。これらのイベント（その正確なタイミングは、モデルの選択によって異なる）の状況を検討すると、面白いだろう。

§4 の点 P が比較的太陽に近い場合、「公転の節約」が高速区間で起きるので、「年」の所要時間は平均より長くなる（§1）。「春分年＝秋分年」のとき、P の近点離角は、おおざっぱに±90°で（正確に±90°ではない）、比較的太陽に近い。従って「春分年＝秋分年」のとき、それは「回帰年」より長い。この結果、2個の曲線の交点は、回帰年を表す曲線（上記においては水平線）より高い位置になり、3種の曲線に囲まれた三角形状の領域が発生する。

PNG画像 (13 KiB)
グラフ3: 「謎三角」の拡大図

この繊細な現象は、0.1秒や0.01秒オーダーの計算精度が確立されて、初めて数値的に検出される。±1秒の精度では「春分年＝秋分年＝回帰年」が計算誤差の範囲内で成立するのか、しないのか、判定が難しい。

そういう目で見ると、グラフ1においても、「謎三角」を視認できる。「Länge des tropischen Jahres」のグラフでも、あるいは「The Lengths of the Seasons」のグラフ（チャート10、PDF版）でさえ、心なしか「謎三角」が発生しているようだ。特に「−4000年」ごろの「冬至年＝夏至年」交点が回帰年より少し上にあることは、比較的分かりやすい。

掘り下げると面白い問題が、いろいろありそうだ！

結び

数学的な意味で、「平均→真」近点離角という中心差の計算は難しい。その方向ではケプラー方程式が解けないため、強引な数値計算や、級数による近似計算が必要になってしまう。「瞬時の年」を定義とすれば、春分年・夏至年などの日数計算は、逆向きの「真→平均」近点離角の変換に帰着される。その方向なら、三角関数・逆三角関数を使って、ケプラー方程式を容易に解くことができる。この三角関数は、近似級数の成分ではなく、厳密解を表す。

われわれは この「逆方向の抜け道」を利用して、春分年・夏至年などの日数について新しい計算法を提示した。結果の数値に新規性はないが、このアルゴリズムは簡潔・高速・高精度で、単なる倍精度演算でも計算自体の誤差が100ナノ秒未満となる（モデルと現実の誤差は、現代±3000年の範囲で0.1秒未満と予想される）。prec() や rate_prec() などはモジュール化されていて、自由に内容を置き換えることができる。「どのような軌道要素・歳差理論を選択するか」とは別の問題として、アルゴリズムの骨組みが示されている。

「春分年」などの概念は便宜上のものであり、現代の天文学において、これらの用語が正式に定義されている わけではない。それどころか、春分点や黄道座標系は、過去の遺物となりつつある のかもしれない…。

付録A　実装例

JavaScript の例。sun_phase = 0, 90, 180, 270 とすると、それぞれ t における春分年・夏至年・春分年・冬至年の日数が返る。sun_phase として任意の太陽黄経を渡せば、その黄経についての「季節年」の日数が返る。どのアルゴリズムも「ケプラー方程式に関する級数を適当に打ち切った近似解」ではなく、正確な解析解を表す。十分に正確な「地球の平均軌道要素と歳差」の式が与えられていることを前提とする。

ソースコードとデモ。

  1. この記事の初版（2017年11月26日）に添付されたバージョンは years-v1-20171126.js と years-v1-20171126.html。
  2. 2017年12月3日: years-v2-20171203.js を公開（変更点）。
  3. 2017年12月10日: years-v3-20171210.js を公開。新表記に合わせた。
  4. 2017年12月31日: years-v4-20171231.js を公開（変更点）。
A1　基本形

最初に（少し冗長だが）基本形を提示する。JavaScript の三角関数・逆三角関数では角度がラジアン単位なので、角度 v は、コード上ではラジアン単位。それとは別に、ケプラー方程式との関連で、vdot（v̇ のこと。見やすいように、本文では v′ としている）もラジアン単位で表す必要がある。

function radians( deg ) { return deg * Math.PI / 180; }
function degrees( rad ) { return rad * 180 / Math.PI; }

function get_year_length_B( t, sun_phase )
{
    var v_deg = ( 180 + sun_phase - (prec(t) + perih(t)) ) % 360;
    var v = radians( v_deg );

    var z = ecc(t);
    var zdot = rate_ecc(t);
    var w = (1 - z) / (1 + z);
    var j = Math.sqrt( 1 - z*z );

    // *1
    var u, E;
    if( v_deg === +180 || v_deg === -180 ) E = Math.PI;
    else {
        u = Math.tan( v/2 );
        E = 2 * Math.atan( u * Math.sqrt(w) );
    }

    // Eq. (6)
    var vdot = radians( -rate_prec(t) - rate_perih(t) );
    var N6 = j * vdot - zdot/j * Math.sin(v);
    var D6 = 1 + z * Math.cos(v);
    var Edot = N6 / D6;

    var mu = degrees( Edot * (1 - z * Math.cos(E)) - zdot * Math.sin(E) );
    var rho = rate_earth(t) - rate_perih(t);

    return 360 / ( rho - mu );
}

v = ±π, ±3π, ±5π, … は関数 tan (v/2) の定義域に含まれないのだから、（2π の整数倍の違いを無視した上で）v = ±π (= ±180°) とそれ以外について場合分けする必要がある。E は後から sin, cos の引数として使われるだけなので、E = ±π の場合、符号の区別は必要ない。

現実の JavaScript では、*1 の7行の分岐の代わりに、単に次のように書いても動作する。というのは、通常の倍精度演算でラジアンを単位とする場合、v/2 = ±π/2 になることはない（π は無理数で、普通の有限小数では ±π/2 という状態を正確に表現できないため）。

    // *1 手抜きバージョン（悪い見本）
    var u = Math.tan( v/2 );
    var E = 2 * Math.atan( u * Math.sqrt(w) );

他の言語（特に三角関数が「°」モードの環境）でこんな手抜きをすると、tan の引数が90°のとき、定義域エラーになる可能性がある。

論理的には、次のように書くこともできる（cos の戻り値が 0 にならない言語では意味がない）。

    var v = radians( 180 + sun_phase - (prec(t) + perih(t)) );
    var u, E;
    if( Math.cos( v/2 ) === 0 ) E = Math.PI;
    else {
        u = Math.tan( v/2 );
        E = 2 * Math.atan( u * Math.sqrt(w) );
    }
A2　前半の別の実装

離心近点離角 E と真近点離角 v について、次が成り立つ［別記事の式 (4.4)］。

atan2 関数がサポートされている環境では、これを利用して次のように書けば、関数の前半を5行にできる（計算精度は、最初の方法と同じ）。

    var v = radians( 180 + sun_phase - (prec(t) + perih(t)) );
    var z = ecc(t);
    var zdot = rate_ecc(t);
    var j = Math.sqrt( 1 - z*z );
    var E = Math.atan2( j * Math.sin(v), z + Math.cos(v) );

この厳密なコードは「手抜きバージョン」よりも1行短い。

普通の atan を使って E を表現する場合には、(♥) に基づくと、角度の象限を決定する手間が発生するばかりか、引数の分母に（比較的分かりにくい）零点があって、話がややこしくなる（E = ±90°）。「基本形」の方法なら、角度の象限が自動決定され、定義域の「穴」（v = ±180°）もよく見えるので、紛れが少ない。

付記: (♥) の両辺を微分して §4 の式 (6) を導くことも可能。

A3　付録Cによる実装

別解の方法により、次のように書くことができる。

function get_year_length_C( t, sun_phase )
{
    var v = radians( 180 + sun_phase - (prec(t) + perih(t)) );
    var z = ecc(t);
    var j = Math.sqrt( 1 - z*z );
    var E = Math.atan2( j * Math.sin(v), z + Math.cos(v) );

    var C0 = 1 - z * Math.cos(E);
    var C1 = j * C0;
    var C2 = 1 + z * Math.cos(v);
    var C3 = degrees( C0/j * Math.sin(v) + C2 * Math.sin(E) );
    var X1 = C1 * rate_prec(t) + (C1-C2) * rate_perih(t) + C3 * rate_ecc(t);

    return 360 / ( rate_earth(t) + X1/C2 );
}

rate_perih() の呼び出しが2回から1回になり、ラジアンへの変換も1回節約された。環境にもよるだろうが、計算精度も少し向上する（倍精度C）。

このコードの（ささいな）問題点は、sin v と cos v がそれぞれ2回ずつ計算されていること。外形的にそれを直すことは難しくない。

A4　長期歳差モデルの実装・他

例えば、prec() と rate_prec() の代わりに prec_LT() と rate_prec_LT() を使えば、長期歳差モデルによる計算になる。

付録B　§4 の式 (6) の導出

§4 の 7. では、次の (5) から (6) への途中計算が省略されている。

  1. E′ = [(1 + u2)v′√w + uw′/√w] / (1 + u2w) ……… (5)
  2. E′ = [jv′ − (e′/j) sin v] / (1 + e cos v) ……… (6)

この変形の主目的は、定義域の「穴」の解消。v = ±π, ±3π, ±5π, … のとき u = tan (v/2) は定義されず、従って (5) の右辺は定義されない。つまり (5) の定義域には、所々に「穴」がある。他方、(5) や (6) の分母が 0 になる心配はない。なぜなら e は地球軌道の離心率で、地球のケプラー軌道が円または楕円であることは既知なので 0 ≤ e < 1。この制約のため、w = (1 − e) / (1 + e) が 0 や負になることはなく、1 + e cos v や j = √1 − e2 が 0 になることもない。

計算の準備

w = (1 − e) / (1 + e) と定義しているので:

一方、(1 + e) > 0 なので (1 + e) = √(1 + e)2 であり:

  1. √w × (1 + e) = √(1 − e) / (1 + e) × √(1 + e)2
  2. = √(1 − e)(1 + e) = √1 − e2

j = √1 − e2 と定義しているので、要するに:

ところで u = tan (v/2) なので、次の関係が成り立つ（倍角の公式、または幾何学的考察による）。

(5) から (6)

式 (5) の分子を N5 としよう。そこに (#1) を代入すると:

  1. N5 = (1 + u2)v′√w + u[−2e′ / (1 + e)2] / √w
  2. = (1 + u2)v′√w − 2ue′ / [√w(1 + e)2]

両辺を (1 + e) 倍して (#2) を使うと:

  1. N5 × (1 + e) = (1 + u2)v′√w(1 + e) − 2ue′ / [√w(1 + e)]
  2. = (1 + u2)v′j − 2ue′ / j

両辺を (1 + u2) で割って (#3) を使うと:

  1. N5 × (1 + e) / (1 + u2) = v′j − 2ue′ / [j(1 + u2)]
  2. = jv′ − (e′/j) sin v ……… (*)

一方、(5) の分母を D5 とすると:

  1. D5 = 1 + u2w = 1 + u2[(1 − e) / (1 + e)]
  2. D5 × (1 + e) = (1 + e) + u2(1 − e) = (1 + u2) + e(1 − u2)

両辺を (1 + u2) で割って (#4) を使うと:

形式的には E′ = N5/D5 = (*) / (**) であり、確かに (6) が成り立つ。

厳密性の確認

式 (6) = (*) / (**) は、式 (5) = N5/D5 と比べ、分子・分母がそれぞれ

倍になっている。例えば v = π のとき u は定義されないのだから、「両辺を (1 + u2) で割って」という処理は実行不可能。その場合でも (6) は有効だろうか？

この疑問を解決するため、(5) の右辺を v の関数 f5(v) と見なし、とりあえずその定義域を 0 ≤ v < π または π < v < 2π としよう。(6) の右辺を v の関数 f6(v) と見なし、その定義域を 0 ≤ v < 2π としよう。v = π 以外の（定義域の）任意の点において、この2関数は等しい値を持つ。なぜなら後者は、前者の分子・分母に一定の数（式 (7) の値）を掛け算したもの。

f(v) = f5(v) − f6(v) と置くと、その右辺の「マイナス」の前と後ろは（v = π 以外では）等しいのだから、v = π 以外の任意の点において f(v) = 0。そればかりか、v → π のときにも次が成り立つ:

これを示すためには、「任意の正の数 ε に対して正の数 δ が存在して、v の値が π から距離 δ 未満（ただし 0 ではないとする）なら、f(v) と 0 の誤差（絶対値）は ε 未満」ということを立証すればいい。ところが、定義域のいたるところで f(v) = 0 であり、従ってこの誤差はどこでも 0。ゆえに、条件を満たす δ は確かに存在する。実際、どんな ε が与えられても、気にせず δ = 1 とでもすればいい。われわれの誤差は 0 なので、どんなに小さい ε が来ても、こっちの方がもっと小さい。

(***) により、v → π における f5 の極限値と f6 の極限値は等しい。ところが f6 は v = π において連続だから、v → π における f6 の極限値は、関数値 f6(π) に等しい。

同様に、f6 の定義域を実数全体として、f5 の定義域を「実数全体から π の奇数倍を除いたもの」とすれば、f6 によって f5 の定義域の「穴」は全てカバーされる。要するに、式 (5) が不定形になる場合も含めて、(5) の代わりに (6) を計算すれば、正しい結果が得られる。

上記の議論では f5, f6 を v の関数と見なした。正確には v は独立変数ではなく、時刻 t の関数 v = g(t)。といっても g は連続関数なので、本質的には何も変わらない。単に g(t) = π となる時刻 t を t1 と置いて、t → t1 における F5(t) = f5(g(t)) と F6(t) = f6(g(t)) の極限値を考えれば、全く同じ議論が成り立つ。

付録C　§4 の別解

§4 の 9. について。

求めたい日数は 360° / (ρ − μ) だが、一方において:

他方において: μ の式は E′ を含み、E′ の式は

を含む。このように ρ, μ の両方が rate_perih(t) を含むのだから、「ρ − μ を簡約できないか」と考えるのは自然だろう。

と置いて、360° / (ρ − μ) の代わりに

  1. 360° / (rate_earth(t) − rate_perih(t) − μ)
  2. = 360° / (rate_earth(t) − ξ) ……… (9)

を計算する形にしたい。

「ラジアン/日」単位で表した rate_prec(t), rate_perih(t) をそれぞれ α, β とする。「ラジアン/日」を単位とすると v′ = −α − β だから、式 (6) をこう書くことができる:

(10) を (1) に代入すると、「ラジアン/日」を単位とした μ は、次のように表現される。

ここで、こう置く。

  1. C0 = 1 − e cos E
  2. C1 = jC0
  3. C2 = 1 + e cos v

すると (11) は、こうなる。

β と (12) を使って (8) を書き換えると、「ラジアン/日」を単位とした ξ は次の通り。

全部の項を C2 倍すると:

  1. C2ξ = C2β + [j(−α − β) − (e′/j) sin v]C0 − C2e′ sin E
  2. = C2β − jC0α − jC0β − (e′C0/j) sin v − C2e′ sin E
  3. = C2β − C1α − C1β − [(C0/j) sin v + C2 sin E]e′
  4. = −C1α + (C2 − C1)β − [(C0/j) sin v + C2 sin E]e′ ……… (13)

(13) の単位は「ラジアン/日」で、後から「°/日」に戻す必要がある。けれど右辺第1項・第2項については、α, β の代わりに「°/日」単位の rate_prec(t), rate_perih(t) をそのまま使えば、「°/日」→「ラジアン/日」→「°/日」の往復変換を回避できる。第3項も「°/日」単位にするため、[ ] 内を 180/π 倍した値を考えよう。それを

と置く。(13) の右辺を「°/日」単位に書き換えたものを X とすると:

  1. X = −C1 × rate_prec(t) + (C2 − C1) × rate_perih(t) − C3 × e′
  2. = −C1 × rate_prec(t) + (C2 − C1) × rate_perih(t) − C3 × rate_ecc(t) ……… (14)

何やら面白い式になったぞ！

μ の式と E′ の式には e′ = rate_ecc(t) が別々に含まれていたが、この散漫さも解消された。

X は、(13) の左辺 C2ξ の値を「°/日」単位に変換したものに等しい。だから、「°/日」を単位とした ξ は ξ = X/C2。これを (9) に入れれば、計算が完了する。(14) の代わりに、符号を反転させた

を考え、(9) の代わりに次の計算をしてもいい。

コードとしては合計10行程度で、本文の方法より少し簡潔。実装例は、付録A3。この場合 v′, E′, ρ, μ を求める必要はない。

文献

数式などの出典

地球の平均軌道要素（対恒星）は、次の 5.8.3 による。
[1] Simon, J. L.; Bretagnon, P.; Chapront, J.; Chapront-Touze, M.; Francou, G.; Laskar, J.: Numerical expressions for precession formulae and mean elements for the Moon and the planets (1994)

一般歳差の通算量は、次の式 (39) による。
[2] N. Capitaine; P. T. Wallace; J. Chapront (2003): Expressions for IAU 2000 precession quantities

これは、下記の Table 4 の pA に当たる。
[2a] N. Capitaine; P. T. Wallace; J. Chapront (2005): Improvement of the IAU 2000 precession model

§5 では、比較検討のための「古い」歳差理論として、下記の Table 8 を利用した。
[3] Laskar, J. (1986): Secular terms of classical planetary theories using the results of general theory

実際のグラフのプロットでは、[2] の代わりに下記の式 (10) と Table 3 を使った。
[4] J. Vondrák; N. Capitaine; P. Wallace (2011): New precession expressions, valid for long time intervals

§5 の検討で利用したのは、Meeus (1991): Astronomical Algorithms にある近似式。

参考文献

Meeus & Savoie (1992): The history of the tropical year は「回帰年」の概念の歴史的変遷をたどり、「春分間隔」などとの違いを論じた軽妙な読み物。

下記では [2] の誤差が評価されている。
A. Fienga; H. Manche; J. Laskar; M. Gastineau (2008): INPOP06: a new numerical planetary ephemeris

ケプラー方程式、中心差、離心近点離角などについて:

追記: 「楕円軌道とケプラー方程式」の式 (2) の数行下の一番長い根号内、第2項の正しい係数は −1。次の点は名前が不統一: C=O, D=G, K=M, H=J。式 (3) の正しい左辺は cos θ。Q’J’ は、添付画像の KL に当たる。P は楕円上の点なので、| PT | は | RT | の b/a。△OKL と △ORT は相似で、各辺の長さが | OK |:| OR | なので、| KL | も | RT | の b/a。従って | PT | = | KL |。準線を用いた証明では「e ≠ 0 なら」という条件が抜けている。e = 0 のケースは別に証明する必要がある。

「春夏秋冬」は「夏秋冬春」より長い > 作成メモ・更新履歴

この記事は、画像・ソースコードも含め、全てパブリックドメインです（再利用自由）。

  1. 2017年8月24日: 「13日は金曜になりやすく31日は水曜になりにくい」の付録「回帰年 vs. 春分の間隔」として書き始めた（グレゴリオ暦の誤差の評価と関連して）。
  2. 2017年9月3日: 「13日は～」を公開。上記の付録は公開せず。
  3. 2017年9月15日: ふと「純粋なケプラー運動から平均春分間隔を計算できる」ことに気付いた。
  4. 2017年9月18日: 独立した記事「「春夏秋冬」は「夏秋冬春」より長い！　地球軌道の幽玄」として再構成開始。歳差だけでなく近日点移動が絡むのが、幽玄だと感じた。
  5. 2017年9月24日: §5 まで一通り下書き。
  6. 2017年9月30日: 春分年の長さを得た（現在の §4）。
  7. 2017年10月9日: §5 まで2回目の下書き。天文計算的な議論を減らし、解析的な解法を追加。この時点では、解析解とは別に、§4 でケプラー方程式を数値的に解いていた。
  8. 2017年10月10日: 旧 §4 を削除。題名も単純に「「春夏秋冬」は「夏秋冬春」より長い」に変更。「幽玄」な部分（近日点移動の影響）については、長くなるのでカット。ケプラー方程式の数値解も、この記事には必要ないのでカット。
  9. 2017年10月15日: アルファ版（未完成の下書き）を暫定公開。
  10. 2017年10月22日: アルファ版（その2）を暫定公開。「古い歳差理論」として Laskar (1986) の式を使い、§5 の検討をやり直した。以前の検討では「Meeus (1991) の計算の背後にある歳差理論は、Simon et al. (1994) の歳差理論と同じ」と勘違いしていて、そのため約70ミリ秒の系統誤差を消せず、ずっとモヤモヤしていた。
  11. 2017年10月29日: ベータ版（作成中の記事）を暫定公開。
  12. 2017年11月1日: ベータ版（その2）を暫定公開。
  13. 2017年11月5日: ベータ版（その3）を暫定公開。グラフを追加。
  14. 2017年11月12日: ベータ版（その4）を暫定公開。§4 式 (5) が不定形になるケースをきちんと扱った。§4 に別解（9.）を追加。グラフの歳差理論を IAU (2006) から Vondrák et al. (2011) に置き換えた。本文は一部未完成。
  15. 2017年11月19日: ベータ版（その5）を暫定公開。本文と付録A～Cを書き終えた。
  16. 2017年11月26日: 初版公開。付録Bの極限の議論を少し簡単化。ソースコードv1を添付。
  17. 2017年12月3日: バージョン2
    1. 計算精度の改善について §4 10. に追記。ソースコードv2を添付。
    2. 表現の細部の改善。「地球軌道上の P から見た太陽の黄経」→「（地球軌道上の）点 P から見た太陽の黄経」。その他、約15カ所を調整、冗長な文を削除（2カ所）。
  18. 2017年12月5日: §4 8. の v′ の値の、指数の誤りを修正。「−4.707（中略） × 10−4」→「−4.707（中略） × 10−5」
  19. 2017年12月10日: バージョン3
    1. 式 (5), (6): 旧版から見ると、それぞれ分子・分母を √b で割った形、分子・分母を j で割った形に変更。この方が (5) から (6) の導出などが簡単になる。ソースコードv3を添付。
    2. 付録B・付録Cの流れを一部変更、C0, C1 などの符号を以前と逆にした（本質的な違いはない）。式 (11) から (13) の導出を簡単化。
    3. 数式の変形の各ステップの説明を title テキストとして付加。中間ステップ（途中計算）の表示も少し増やした。
    4. 表現の細部の改善（約15カ所）。
  20. 2017年12月11日: §3 の誤字修正。「春分間隔、冬至間隔も同様。」→「秋分間隔・冬至間隔も同様。」
  21. 2017年12月15日: §4 の不正確な記述を改めた。「v = π (= 180°) の場合 u = ∞ で…」→「v = π (= 180°) の場合 u は定義されず…」
  22. 2017年12月17日: バージョン3.1
    1. §4 3. に解説を追加。
    2. Ex3.4 のすぐ次のパラグラフ: 「73048.47501…日」→「73048日11時間24分01秒」のように単位を「時・分・秒」に。
    3. 「参考文献」に追記。
    4. 表現の細部の改善（約15カ所）。
    5. 根号の中身にCSSでオーバーラインを引くようにした。
  23. 2017年12月24日: 付録B (#3)(#4) の導入を別の記事にした。
  24. 2017年12月31日: バージョン4
    1. §4: b = (1+e)/(1−e) と置いて √b で割る代わりに、w = (1−e)/(1+e) と置いて √w 倍するようにした。微分が少し簡潔になり、倍精度Aではわずかに精度も向上。ソースコードv4を添付。
    2. 式 (2) を (4.4) からではなく (4.5) から導く形に変更。
    3. 表現の細部の調整（約5カ所）。一部の画像を描き直した。
  25. 2018年1月10日: 根号のオーバーラインの一部を border-top に変更。
  26. 2018年1月14日: 関連記事「ケプラー方程式（微積・三角公式を使わないアプローチ）」を公開したのでリンク設定。
  27. 2018年1月21日: 表現の改善（1カ所）。「長軸が楕円の周と交わる2点」→「長軸の両端の2点」
  28. 2018年12月30日: 「太陽位置の簡易計算」: 「1000年や1万年の平均を求めたら、さらに正確になる？」という点について、注釈を付けた。
  29. 2020年11月24日: 離心率を表す文字を z から e に変更。
  30. 2022年7月19日: 修正。本文中、離心率を表す文字が一部 z のままで不統一だった。ソースコードの変更はバグのもとなので、コード上は引き続き z のままにしておく。
  31. 2022年9月1日: 表現の細部の調整（3カ所）。

この記事のURL

パブリックドメイン

← 前 1‒20 p ｜ ↑ もくじ i ｜ 次 1‒22 n →
<メールアドレス>
▲
［妖精現実 フェアリアル］ Home h > Articles z > §1 i > 1‒21