天文計算 — 星はどこにあるか

「実際の夜空を描く」というとき、計算しているのは一つのことです。ある地点・ある時刻に、その天体がどちらの方向に見えるか。この節では、その方向をどんな量から導いているか、恒星と太陽系天体でなぜ計算の経路を分けたか、どこまでを既存のライブラリに委ねたかを説明します。

その計算はすべて純粋な計算層に閉じています。分けている理由はテストです。外部の状態に触れず、ブラウザもWebGLも要らない関数の集まりにしておけば、Node上のテストランナーでそのまま走り、天文年鑑（その年の天体の位置を表にした刊行物。暦書とも呼びます）に載る値や、別実装の計算結果と直接突き合わせられます。

天球と、二つの座標系

天体は、観測者を中心とする半径一定の仮想的な球面、すなわち天球の上の一点として扱います。Asterariumのシーンでも天球は実在していて、恒星は半径1000の球面上に置かれます。シーンの長さは1単位が1メートル相当と決めてあるので（この取り決めは「座標系とシーンの姿勢」の節にあります）、これは半径約1kmの球面ということになります。この数字は見え方には影響しません。半径は星の向きさえ保てばよいからです。条件は一つだけ、カメラが描く範囲の内側に収まることです。その範囲は、手前がニアプレーンの2単位、つまり2メートル。奥がファープレーンで、1e7単位、つまり1千万単位、およそ1万kmです。

天球上の一点を表す角度の組は、基準の取り方で何通りもあります。使うのは二つです。赤道座標は地球の自転軸を基準に取り、赤経と赤緯で位置を表します。地平座標は観測者の足元を基準に取り、方位（北が0°、東が90°）と高度（地平線が0°、天頂が+90°）で表します。

二つを分けるのは、時間に対するふるまいです。恒星の赤経・赤緯はほとんど変わりません。地球の自転が変えるのは観測者の向きであって、星そのものの位置ではないからです。逆に地平座標は数分で目に見えて変わります。星のカタログが赤道座標で配布されるのはこのためで、Asterariumはその値をそのまま使います。

二つを結ぶのが恒星時、春分点を基準にした地球の自転角です。地球は太陽に対して24時間で一周しますが、遠い星に対しては約23時間56分で一周します。その差の約4分ぶん、同じ星が南中する時刻は毎日約4分ずつ早まり、1年でほぼ一周ぶんずれます。純粋な計算層はグリニッジの恒星時を天文ライブラリから取り、観測者の経度を足して地方恒星時にします。ただし、恒星を画面に出す経路はこの値を通りません。天球を回す行列は、大気散乱ライブラリ@takram/three-atmosphereが日付から作るからです。大気の見え方を計算するには地球の姿勢そのものが要るので、そのライブラリが地球の姿勢を表す行列をすでに持っている、という事情です。恒星時を明示的に使うのは、情報パネルの表示・天体の検索・照準のための簡易な変換と、投影機のページ（/projector/）の計算です。

見かけの位置を歪める三つの効果

幾何学的な位置と、実際に目に見える位置は少しずれます。主な原因は三つで、Asterariumは対象ごとに扱いを変えています。ずれの大きさは秒角と分角で表します。見えるかどうかの目安には二つの水準があります。一つは、点が1つずれたと分かる限界で、およそ1分角、つまり60秒角。もう一つは、星座の形が変わったと分かる限界です。こちらは数十度に広がる星の並びの中での話で、星座線で結ばれる隣り合う星の間隔は目安として数度あります。その中では数分角のずれにはまだ気づけず、数十分角に達してようやく形の違いとして見えてきます。

• 光行差 — 観測者自身が動いているため、光の来る向きが進行方向へわずかに傾く現象。地球の公転による分は最大で約20秒角です。太陽系天体では計算に含めています。恒星ではカタログのJ2000位置をそのまま使い、補正しません。
• 大気差 — 大気の屈折で天体が実際より高く見える現象。地平線付近で約34分角に達し、太陽や月の見かけの直径、およそ0.5度ほどもあります。太陽系天体では標準的な補正を入れて地平座標へ変換します。恒星の描画では入れていません。
• 視差 — 観測地が地球の中心ではなく地表にあることによる方向のずれ。月では最大約1°に達するので無視できません。観測者の緯度・経度・標高を渡し、トポセントリックな位置として計算します。恒星は遠すぎて完全に無視できます。

恒星に光行差と大気差を入れない理由は、それぞれ違います。光行差の約20秒角は、点1つのずれとして見える限界の3分の1しかなく、そもそも見えません。大気差はそうではありません。高度45°でも約1分角、地平線の近くではその30倍を超える約34分角に達します。それでも入れていないのは、大気差の量が高度によって変わるからです。高度ごとに量の違うずれは、天球をまるごと一つの回転で回す現在の経路では表せません。数万個の恒星は、まさにその一つの回転行列で回すだけで画面に出ます。行列の中身は、2000年の基準から当日の基準への歳差に地球の自転を掛けたものです（その組み立て方は「座標系とシーンの姿勢」の節にあります）。

とはいえ、道が閉じているわけではありません。星を描く頂点シェーダーは、天球の回転行列を通した後の各星の向きをすでに計算していて、その上下成分は高度の正弦そのものです。実際そこでは、きらめきの強さを高度によって変えるのにこの値を使っています。同じ値から高度ごとの持ち上げ量を求めて星を上へずらすことは、原理的にはできます。実装していないだけです。

一方、太陽・月・惑星は1天体ずつ計算するので、大気差を入れています。この選び分けには、良い面と悪い面が両方あります。太陽そのものの高度は実際の空に近づきます — 沈む太陽が地平線の上に浮き上がって見えるのと同じ効果です。一方で背後の星との相対位置は、星に大気差を入れていないぶん実際の空とずれ、地平線の近くでは太陽や月が最大約34分角だけ星の並びより高く描かれます。後者を承知のうえで残しています。

二つの経路 — 恒星はビルド時に、太陽系天体は毎フレームに

恒星と太陽系天体では、位置の決まり方がまったく違います。恒星は数万個あるけれども、この用途では動かないとみなせる。このアプリが扱う太陽系天体は、太陽と月に、地球を除く7つの惑星、つまり水星・金星・火星・木星・土星・天王星・海王星を足した9個しかないけれども、速く動き、しかも恒星よりも高い精度 — 肉眼の限界にあたる分角級 — が要る。精度が要るのは、点で描く恒星と違って見かけの大きさを持つからです。月は直径約0.5度の円盤で、その欠け方も、惑星や恒星への近づき方も、数分角の差で目に見えて変わります。この非対称が、そのまま二つの経路になっています。

恒星の位置は配信前に確定します。決まるのはビルド時、利用者がページを開くよりも前です。ビルドスクリプトが恒星カタログAT-HYG（古くからあるHYG星表を、Tychoの測光やGaiaの視差といった観測で拡張した公開のカタログ）を読み、各星のJ2000赤道座標を単位ベクトルへ変換し、二進のファイルへ焼き込みます。ランタイムでは、1星ずつ位置を計算しません。天球全体を一つの行列で回すだけです。だから星が何万個あっても、位置計算の仕事は行列を一つ組み立てることに変わりません。描画そのものの負荷は星の数に比例するので、そちらは端末の性能に応じて描く数を絞ります。絞るのは品質段階の役目で、明るい順に並べておいて先頭から必要な数だけ描くというデータの作り方は「データパイプライン」の節にあります。

太陽系天体は逆に、毎フレーム計算します。太陽・月と、地球を除く水星から海王星までの7惑星、あわせて9天体について、bodyState()という関数が観測地と時刻から地平座標・当日基準の座標での赤経赤緯・等級・視半径を返します。太陽以外の8天体には、位相角と輝面比も付きます。9天体分をまとめて呼んで約0.12ミリ秒。開発時に手元の機械で計った、描画を含まない位置計算だけの値で、端末が変われば変わります。比べる相手は、描画まで含めたフレーム全体に許された時間 — 毎秒60回なら16.7ミリ秒 — です。位置計算はその1%に満たないので、フレームを間引くような小細工はしていません。

const time = MakeTime(new Date(t))
const observer = new Observer(loc.latDeg, loc.lonDeg, loc.elevationM)
const aeBody = BODY_MAP[body]

// Topocentric equatorial coordinates, of-date, with aberration.
const eq = Equator(aeBody, time, observer, true, true)
const hor = Horizon(time, observer, eq.ra, eq.dec, 'normal')
const illum = Illumination(aeBody, time)

// eq.dist is topocentric distance in AU. Angular radius = asin(R / d).
const distKm = eq.dist * KM_PER_AU
const angularRadiusDeg = (Math.asin(RADIUS_KM[body] / distKm) * 180) / Math.PI

const state: BodyState = {
  horizontal: { azDeg: hor.azimuth, altDeg: hor.altitude },
  equatorialOfDate: { raHours: eq.ra, decDeg: eq.dec },
  mag: illum.mag,
  angularRadiusDeg,
}
// Then: phase angle and illuminated fraction, for every body but the Sun.

短いですが、この節で挙げた三つの補正がそのまま並んでいます。観測者を緯度・経度・標高から作るので位置は地表基準になり、これが視差の補正にあたります。赤道座標を求める呼び出しの末尾の二つのtrue引数は、2000年時点ではなく当日の赤道と春分点を基準にすること、つまり当日基準の座標で返すことと、光行差を含めることを指示します。地平座標へ移す呼び出しのnormalは標準的な大気差モデルの指定です。視半径だけは自前で、天体の半径と距離から求めています。

基準の取り方が恒星と食い違わないのか、という疑問が残ります。恒星は2000年時点の赤道と春分点を基準にした座標のまま、太陽系天体は当日基準の座標です。答えは、二つが別々の道を通って同じ当日の空に着く、です。恒星は天球全体を回す行列を通り、その行列が2000年の基準から当日の基準への歳差ぶんも一緒に動かします。太陽系天体はこの行列をまったく通りません。bodyState()が返す地平座標を、そのままワールド座標の方向ベクトルに直して置くだけです。地平座標はその時点の空を基準にした量なので、最初から当日の基準になっています。

その行列は大気散乱ライブラリが日付から作ります。大気の計算には地球の姿勢が要るので、ライブラリはこの行列をもともと持っていて、外へ公開しています。中身は前に触れたとおり、2000年の基準から当日の基準への歳差に、地球の自転を掛けたものです。歳差の量は2026年時点で約0.4度あり、星座の位置が肉眼で分かる程度に動く大きさなので、これが行列に入っているかどうかは効いてきます。

恒星に固有運動を入れていないのも、その量が見えないからです。ほとんどの星の固有運動は年0.1秒角未満で、カタログの基準である2000年から2026年時点で26年たっても、移動は数秒角にとどまります。AT-HYGから抜き出しているのは8等より明るい約38,000星です。8等まで持つのは、実際に描く下限の6.8等に1.2等ぶんの余裕を足した数だからです（カタログをどう二段に分け、どの設定がどこまで描くかは「データパイプライン」の節にあります）。そのなかで最も速いのがグルームブリッジ1830の年約7.1秒角、明るい星ではアルクトゥルスの年約2.3秒角です。26年ぶんに直すと約3分角と約1分角。数十度の広がりを持つ星座の形が変わって見える大きさではありません。

同じ判断が、目的の違うページでは逆転します。恒星間フライトのページ（/starflight/）は固有運動そのものが主題なので、実測の固有運動と視線速度を持つ別のカタログを作り、前後10万年まで外挿します。無視してよいかどうかは星の性質ではなく、そのページが何を見せたいかで決まる、ということです。

どこまでをライブラリに任せたか

天体位置のエンジンはastronomy-engine一つだけです。決め手は三つありました。MITライセンスなので、サイトに含めて配布できること。級数展開だけで自己完結していて、あらかじめ計算した位置を並べた表を読み込まずに済むこと。そして公表されている精度が太陽系天体で1分角級、つまり肉眼の限界とほぼ同じであること。サーバー側で動くコードを持たないという前提に噛み合うのは、このうち前の二つです。ブラウザの中だけで計算が完結し、追加のダウンロードも増えません。精度のほうはサーバーの有無とは無関係で、肉眼で見える空を描くのに足りるかどうかという別の基準で測っています。ただし、要る精度も公表精度も1分角級だというのは、余裕が大きくないということでもあります。だから、分角の差がそのまま絵に出る日食・月食のページ（/eclipse/）は、この汎用の位置計算を使いません。同じライブラリが持つ日食専用の探索を呼んで、観測地ごとの接触時刻を直接得ます。毎フレームではなく、観測地を選んだときに一度だけ走る計算です。

ライブラリに委ねているのは、当日基準の座標での、地表を基準にした赤経赤緯、大気差込みの地平座標への変換、等級と位相角と輝面比、恒星時、出没時刻と薄明時刻の探索、月相、月の秤動です。自前で書いたのはその上に載る薄い層で、角度から単位ベクトルへの変換、月の輝縁の位置角、1日分のイベントのまとめ方、そして低精度な地平座標変換です。

低精度版をわざわざ別に持っているのは、用途が違うからです。情報パネルの表示、天体の検索、そしてカメラをその天体へ向けるための照準は、恒星や星団のようにほぼ動かない対象について、赤道座標から地平座標への標準的な回転で方位・高度を出せば足ります。この回転を決めるのは観測地の緯度と地方恒星時だけです。大気差も光行差も入れません。太陽系天体だけはこの経路を使わず、必ずbodyState()を通す決まりです。恒星向けの簡易な地平座標変換と、太陽系天体向けのbodyState()を並べて持つ以上、どちらを何に使うかを決めておかないと取り違えが起きます。

等級、月、そして薄明

星の明るさは等級で表します。数字が小さいほど明るく、5等級の差がちょうど光量100倍と定義されています。1等級あたりの比はその5乗根、約2.512で、ポグソン比と呼ばれます。カタログが持つ等級を描画に使える線形の強度へ直すのは、次に挙げる関数のreturnの行、式にして一行です。やっていることは言葉にすれば一つで、等級が1つ増えるごとに強度を約2.512分の1にし、基準の等級でちょうど1にする、というだけです。その基準の等級は引数で選べます。星を描くシェーダーは同じ式を6.5等基準で使い、暗い空で肉眼が見える限界あたりの星を強度1に置いています。

const POGSON = Math.pow(100, 1 / 5) // 2.511886431509…

export function magToIntensity(mag: number, refMag = 0): number {
  return Math.pow(POGSON, refMag - mag)
}

月については、位相角と輝面比のほかに、輝縁の位置角と秤動も計算しています。欠けた月の明るい側が画面のどちら向きになるかは、この輝縁の位置角で決まり、太陽と月の赤道座標から古典的な式で求めます。ただし満ち欠けの向きを決めるときに位相角は使いません。位相角は満月付近で0、新月付近で180になり、上弦でも下弦でも同じ90前後の値を取るので、行きと帰りを区別できないからです。代わりに月と太陽の黄経差 — 新月で0°、満月で180°、次の新月で360°と一方向に増える量 — を使います。

薄明は太陽の高度で三段階に分けます。市民薄明が−6°、航海薄明が−12°、天文薄明が−18°です。その日のイベントをまとめる関数は、日出没・月出没とあわせて、六つの薄明時刻 — 夕方に太陽の高度が各閾値を下向きに横切る時刻と、朝に上向きに横切る時刻 — を探索します。空の暗さそのものは大気散乱の計算が決めますが、いつ何が起きるかという時刻の側はここが持ちます。

高緯度では、太陽が一日中沈まない日も、一日中昇らない日もあります。そのときイベントの時刻は存在しません。純粋な計算層はそれを0や当日の終わりで埋めることをせず、フィールドを未定義のまま返します。呼び出し側に「その日は起きなかった」ことが伝わるので、白夜や極夜の空でも表示が破綻しません。存在しない値を無理に作らない、という規律です。

この節の判断は、どれも一つの問いに集約されます。どれだけの精度が、どの頻度で必要か。肉眼で見分けられない差は計算しない。見える差でも、描画の経路に載らないなら入れない。この線引きが、サーバー側で動くコードを持たない静的なページのまま、ブラウザだけで毎秒60回の空を描く余裕を作っています。