金星・水星の最大光度のタイミングの違いを解く

この記事は、数ヶ月前にニューラルネットワークをScalaでゼロから書いてたときに勾配降下法の応用としてやっていた内容です。

金星・水星の光度変化を数式で表現し、その最大値のタイミングを調べました。最大値は勾配降下法を使わなくても解けそうですが、せっかく勾配降下法を作ったので試したものです。

金星・水星の最大光度とは

水星と金星は地球よりも内側を周回している惑星です。内惑星といいます。内惑星は地球から見ると月と同様に満ち欠けをします。

ただし、内惑星の満ち欠けが月と大きく違うのは、地球との距離が大きく変わるため見た目の大きさが変わることです。内惑星が満月状のときが地球から最も遠く、見た目も最も小さくなります。新月状のときが地球に最も近く、見た目も最も大きくなります。

満月状のときを「外合」、新月状のときを「内合」といいます。

f:id:suzuki-navi:20200906221643p:plain

月が最も明るくなるのは満月のときですが、内惑星はいつ最も明るくなるのでしょうか。遠いけど満月状のとき（外合）か、細い姿だけど地球から近いとき（内合前後）か。

実際は、2つの内惑星で光度変化の挙動が大きく異なります。水星は遠いにもかかわらず満月状のとき（外合）が最も明るくなり、金星は細いけど地球に近い三日月状のとき（内合の前後）が最も明るくなります。

水星は満月状のとき（外合）が最も明るくなります。外合は太陽と同じ方向で、太陽が眩しいため、最も明るいにも関わらず水星を人が観測することはできません。

金星は三日月状のとき（内合の前後）が最も明るくなります。金星が最大光度のころは、宵の明星または明けの明星として、夕方または日の出前にとても明るく輝いて見えます。最大光度の金星は、太陽、月についで明るく見える天体です。

この記事では、内惑星の光度を適当な式で表現し、その式から勾配降下法で光度が最大になるタイミングを計算しました。

式で表現したならばその微分が0になるところを直接計算してもいいのですが、勾配降下法の例題として考えたので、この記事ではその計算していません。

計算式を考える

位置

まず、位置関係についての計算式です。

太陽を座標の原点として、視点となる地球の位置を \displaystyle x=1, y=0 に固定して、平面上を内惑星が円軌道で周回するとします。

\displaystyle r_{vs} は内惑星ごとの定数です。 \displaystyle \theta_s は内惑星の位置を表すパラメータです。他の変数の値はすべてこの2つから計算できます。 \displaystyle \theta_s=0 が内合、 \displaystyle \theta_s=\pi が外合です。

  • \displaystyle \theta_s : 太陽から見て地球と内惑星がなす角
  • \displaystyle \theta_e : 地球から見て太陽と内惑星がなす角
  • \displaystyle \theta_v : 内惑星から見て太陽と地球がなす角
  • \displaystyle x_v, y_v : 内惑星の位置の座標
  • \displaystyle r_{vs} : 内惑星と太陽の距離（日心距離）（一定値）
  • \displaystyle r_{ve} : 内惑星と地球の距離
 \displaystyle
\begin{align}
x_v &= r_{vs}\cos \theta_s \\
y_v &= r_{vs}\sin \theta_s =  r_{ve} \sin \theta_e \\
{r_{ve}}^2 &= (1-x_v)^2+y_v^2 \\
y_v &= r_{ve} \sin \theta_e \\
\pi &= \theta_s + \theta_e + \theta_v \\
\end{align}

実際には地球も内惑星も楕円軌道を描き、もっとずっと複雑です。単純化のため円軌道を仮定しています。

明るさ

次に、明るさに関する計算式です。明るさを計算する一般的なモデルはわかりません。この記事では次の計算式に従うものと仮定しました。実際には惑星の表面や大気の構成によっても変わりそうです。

  • \displaystyle a : 距離を一定にし、内惑星の見た目の形( \displaystyle \theta_v )から計算される明るさで、満月状のときを基準にした相対値
  • \displaystyle b : 地球からの距離と \displaystyle a から計算される明るさ（距離1で満月状の場合に明るさ1とした相対値）
 \displaystyle
\begin{align}
a &= (1 - \frac{\theta_v}{\pi}) \cos \theta_v + \frac{1}{\pi} \sin \theta_v \\
b &= \frac{a}{{r_{ve}}^2} \\
\end{align}

\displaystyle a の式はちょっと複雑な形をしていますが、次の積分を計算して式を作りました。この積分も私が適当にモデルを考えたもので、合っているのかどうかはよくわかりません。

 \displaystyle
\begin{align}
a &= \frac{2}{\pi} \int_0^{\pi-\theta_v} \sin \theta \sin(\theta_v + \theta) d \theta \\
&=  \frac{2}{\pi} \int_0^{\pi-\theta_v} \sin \theta ( \sin \theta_v \cos \theta + \cos \theta_v \sin \theta ) d \theta \\
&=  \frac{2}{\pi} \int_0^{\pi-\theta_v} \sin \theta_v \sin \theta \cos \theta + \cos \theta_v (\sin \theta)^2 d \theta \\
&=  \frac{2}{\pi} \int_0^{\pi-\theta_v} \frac{1}{2} \sin \theta_v \sin 2\theta + \frac{1}{2} \cos \theta_v (1 - \cos 2\theta) d \theta \\
&=  \frac{2}{\pi} \left[ -\frac{1}{4} \sin \theta_v \cos 2\theta + \frac{\theta}{2} \cos \theta_v - \frac{1}{4} \cos \theta_v \sin 2\theta \right]_{\theta=0}^{\theta=\pi-\theta_v} \\
&=  \frac{2}{\pi} \left[ \frac{\theta}{2} \cos \theta_v - \frac{1}{4} \sin (\theta_v + 2\theta) \right]_{\theta=0}^{\theta=\pi-\theta_v} \\
&= \frac{2}{\pi} \left( \frac{1}{2} (\pi - \theta_v) \cos \theta_v - \frac{1}{4} \sin (2\pi - \theta_v) + \frac{1}{4} \sin \theta_v \right) \\
&= \frac{2}{\pi} \left( \frac{1}{2} (\pi - \theta_v) \cos \theta_v + \frac{1}{2} \sin \theta_v \right) \\
&= (1 - \frac{\theta_v}{\pi}) \cos \theta_v + \frac{1}{\pi} \sin \theta_v \\
\end{align}

ちなみに \displaystyle a は \displaystyle \theta_v の関数ですが、グラフにすると次のような形をしています。1つ目は横軸 \displaystyle \theta_v の単位をラジアンにしたもの、2つ目は横軸 \displaystyle \theta_v の単位を度にしたものです。

f:id:suzuki-navi:20200906221715p:plain

f:id:suzuki-navi:20200906221727p:plain

ここまでの式をもとに、明るさ \displaystyle b は \displaystyle \theta_s から以下の順番に計算できます。 \displaystyle r_{vs} は定数とします。

 \displaystyle
\begin{align}
x_v &= r_{vs}\cos \theta_s \\
y_v &= r_{vs}\sin \theta_s \\
{r_{ve}}^2 &= (1-x_v)^2+y_v^2 \\
r_{ve} &= \sqrt{r_{ve}^2} \\
\theta_e &= \sin^{-1} \frac{y_v}{r_{ve}} \\
\theta_v &= \pi - (\theta_s + \theta_e) \\
a &= (1 - \frac{\theta_v}{\pi}) \cos \theta_v + \frac{1}{\pi} \sin \theta_v \\
b &= \frac{a}{{r_{ve}}^2} \\
\end{align}

光度の変化をグラフで見る

\displaystyle r_{vs} に実際の水星と金星の値を入れてみて、内惑星の位置 \displaystyle \theta_s と明るさ \displaystyle b の関係をグラフにしました。1つ目は水星、2つ目は金星です。横軸は \displaystyle \theta_s の単位を度にしたものです。

f:id:suzuki-navi:20200906221749p:plain

f:id:suzuki-navi:20200906221803p:plain

水星と金星とでグラフの形が大きく違います。

勾配降下法を使う

最も明るくなる箇所を勾配降下法で探します。

勾配降下法は、与えられた関数の勾配（微分）を数値的に計算して、関数が最大または最小となる箇所を探す方法です。

微分の値を計算するには、式全体を微分してから数値を当てはめればいいのですが、複雑な式を微分するのは面倒なので、ニューラルネットワークのバックプロパゲーション（誤差逆伝播法）と同様に、またTensorFlowの自動微分と同じように、計算ノードごとに微分を計算することにします。

いくつかの微分式をここに書いておきます。

 \displaystyle
\begin{align}
\frac{d}{dx} \sin x &= \cos x \\
\frac{d}{dx} \cos x &= -\sin x \\
\frac{d}{dx} \sin^{-1} x &= \frac{1}{\sqrt{1-x^2}} \\
\frac{d}{dx} \frac{1}{x} &= -\frac{1}{x^2} \\
\frac{d}{dx} \sqrt x &= \frac{1}{2\sqrt{x}} \\
\end{align}

これで、定数 \displaystyle r_{vs} に対して \displaystyle b が最大となる \displaystyle \theta_s を勾配降下法で計算できます。

内惑星の日心距離と最大光度のタイミングの関係を見る

内惑星の日心距離（太陽と内惑星の距離 \displaystyle r_{vs} ）から、内惑星の光度（ \displaystyle b ）が最大となる \displaystyle \theta_s を勾配降下法により計算します。

横軸に \displaystyle r_{vs} 、縦軸に最大光度の \displaystyle \theta_s を置いたものが以下のグラフの青い線です。ついでにそのときの \displaystyle \theta_e 、 \displaystyle \theta_v の線もひいています。

f:id:suzuki-navi:20200906221822p:plain

青: 太陽から見て地球と内惑星がなす角 \displaystyle \theta_s 橙: 地球から見て太陽と内惑星がなす角 \displaystyle \theta_e 緑: 内惑星から見て太陽と地球がなす角 \displaystyle \theta_v

日心距離0.45付近に大きな崖があります。

日心距離が0.45付近より小さい内惑星では、 \displaystyle \theta_s=180° つまり外合のときに最大光度となることがわかります。日心距離がそれより大きくなると、最大光度となる \displaystyle \theta_s が急激に内合のほうに近くなることがわかります。

日心距離が閾値より小さい内惑星は、地球から見ると太陽のすぐ近くを小さく周回しており、地球との距離の変化が小さくなります。そのため、月と同じように満月状のときが最大光度になります。

逆に日心距離が閾値より大きい内惑星は、地球から見ると太陽の周りを大きく周回しており、地球に近いときは極端に近くなります。地球との距離の変化が激しくなるため、距離のほうが光度の変化に大きく寄与することになり、内合に近いタイミングで最大光度になります。

水星は日心距離0.3871で0.45付近のしきい値より太陽に近いため、外合のときが最大光度となります。金星は日心距離が0.7233ですので、内合に近いときが最大光度となります。

惑星の表面の違いによる光度変化の違いを考察する

今回明るさを計算するモデルとして以下の式を使いました。

 \displaystyle
a = (1 - \frac{\theta_v}{\pi}) \cos \theta_v + \frac{1}{\pi} \sin \theta_v \\

実は最初はもっとシンプルな以下の式を使っていました。

 \displaystyle
a = \frac{1 + \cos \theta_v}{2}

\displaystyle a はグラフにすると以下のような形をしています。

f:id:suzuki-navi:20200906221843p:plain

f:id:suzuki-navi:20200906221854p:plain

もっと複雑な最初に紹介した \displaystyle a のグラフを再掲します。

f:id:suzuki-navi:20200906221749p:plain

f:id:suzuki-navi:20200906221803p:plain

見た目ほとんど同じですが、シンプルな \displaystyle a のほうが、真ん中あたりの値が少し大きいです。

シンプルな \displaystyle a を使って計算したグラフを並べておきます。

水星の光度変化。

f:id:suzuki-navi:20200906221947p:plain

金星の光度変化。

f:id:suzuki-navi:20200906222002p:plain

日心距離と最大光度の位置の関係。

f:id:suzuki-navi:20200906222015p:plain

シンプルな \displaystyle a のほうが閾値が小さいです。

金星はどっちであっても結果に変わりありませんが、水星の光度変化は \displaystyle a をシンプルにすることで大きく変わってしまいました。事実とはだいぶ違ってしまっています。

\displaystyle a の計算式のモデルの説明は省略しますが、シンプルな \displaystyle a は濃い大気に覆われている天体を想定したモデル、複雑な \displaystyle a は大気が少なく地表がクレーターなどごつごつしている天体を想定したモデルです。あくまで素人が勝手に考えたモデルなので、違っている可能性高いのですが、水星の光度変化が水星の表面を推測する根拠になりうることには驚きました。

リンク

過去の勾配降下法に関する私の記事