太陽の位置の季節変化「アナレンマ」を計算

2018年3月21日は[春分の日](https://ja.wikipedia.org/wiki/春分の日)。昼と夜がちょうど同じ長さ…ではないのだが、太陽の動きを観察したとき1年間の節目のひとつになっている。

毎日同じ時刻（例えば昼の12時）に太陽の位置を記録したら、1年間で空にどんな軌跡を描くだろうか。太陽高度は冬に低く夏に高いから、上下方向の直線になりそうな気がする。

実は、太陽は上下（南北）だけでなく左右（東西）にもずれる。例えば冬至は昼が最も短いが、日の出は冬至を過ぎてもしばらく更に遅くなる。このずれが加わることで、太陽の軌跡は大きな**8の字**になることが分かっている。**「アナレンマ」**（analemma）と呼ばれるこの軌跡を求めてみたい。

[Google画像検索：「analemma photo」](https://www.google.com/search?tbm=isch&q=analemma+photo)

## 目的
アナレンマを計算してグラフで表す。

* 計算の天文的意味を知る
	* 用いるパラメーターを最小限に絞る
	* 歳差・章動・近点移動などの、長期的または微小振幅の変化は無視する
	* 軌跡を描くことに集中し、時刻の関わる計算は実行しない
* 「地球」における「太陽」のアナレンマ*以外*も描く
	* 必要なパラメーターの収集方法をまとめる

## 8の字ができる理由
アナレンマの南北方向の変化は普通に**太陽高度の季節変化**であり、東西方向の変化は**均時差**と呼ばれるものである。均時差と言うと専門的で難しそうだが、機械時計に対する日時計のずれとして簡単に確かめられる[^eot-corr]。ただ、これらの変化は独立ではないので、8の字の要因でなく結果と考えたほうが理解しやすい。

では8の字は何なのかというと、太陽が地球を1年かけて回っているという見方において[^sun-revolution]、天の赤道上を等速円運動する仮想的な太陽（平均太陽：下図の赤色）を基準としたときの実際の太陽（視太陽：下図の橙色）の**相対位置の軌跡**として得られる。平均太陽は季節変化などない完璧に一定の動きをしていて、地上から見れば毎日同じ時刻に空の同じ位置にいる。（仮想なので見えないが）

<!-- ![analemma_sun.png](https://qiita-image-store.s3.amazonaws.com/0/127130/efec7ca7-1e63-cdaf-fafb-aaa02417a7c2.png) -->
<!-- ![analemma_sun.png](https://qiita-image-store.s3.amazonaws.com/0/127130/eed22450-47a0-52fe-ac2a-28f4caecf010.png) -->
![analemma_sun.png](https://qiita-image-store.s3.amazonaws.com/0/127130/95285bb5-6529-29cf-60c5-8be1d6817bdb.png)

この理屈は、準天頂衛星（みちびきなど）の[8の字軌道](http://www.jaxa.jp/countdown/f18/overview/orbit_j.html)と同じことである[^geosync]。人工衛星ならパラメーターを変える想像がしやすいので、**3つの軌道要素**が8の字の形を決める様子を見ていく。

[^eot-corr]: 逆に言えば、日時計で正確に時刻を知るには均時差の補正が必要である。
[^sun-revolution]: 1年間の地球の公転運動を地球中心に見た状態で、地球からは太陽が天球上の星々の合間をゆっくり進むように見える（→[黄道](https://ja.wikipedia.org/wiki/黄道)）。混同してはいけないこととして、毎日の太陽の出没は、太陽の周回でなく地球自身の自転によるものと考え続ける。
[^geosync]: 地表に対して8の字と説明されるが、静止衛星を基準に8の字とも考えることができる。

<!-- ![analemma_satellite_0.png](https://qiita-image-store.s3.amazonaws.com/0/127130/3bfa6d7e-9142-01af-8612-e30b379b5fef.png) -->

### 軌道傾斜角
地球の自転と同じ周期で等速円運動する2機の人工衛星を用意し、ひとつは赤道上空（→静止衛星）、もうひとつは**赤道面に対して斜め**に配置する（→準天頂衛星の候補）。どちらも速さは常に同じなのだが、同時刻における相対位置は記録していくと時間変化が現れる。

* 赤道面を横切るときは、斜めに進むため静止衛星より**遅れる**＝相対的に**西へ移動**していく
* 北端・南端付近のときは、自転軸に近づいているので静止衛星より**進む**＝相対的に**東へ移動**していく

これと南北移動を合わせると、右図のように静止衛星を基準とした8の字軌道に見えることになる。

<!-- ![analemma_satellite_1.png](https://qiita-image-store.s3.amazonaws.com/0/127130/c08bd7bb-7ece-24eb-ab08-48cd4ee62bf6.png) -->
![analemma_satellite_1.png](https://qiita-image-store.s3.amazonaws.com/0/127130/9223c97d-8807-50c9-376d-7996d0cef200.png)

### 軌道離心率
傾斜軌道の衛星をコピーし、軽い**楕円軌道**に変える（南端で地球に最接近し、北端で最遠になるように）。楕円軌道の場合は地球に近づくほど速く動くので、角速度が一定でなくなり、

* 南端付近では速くなり、静止衛星に対する追い越しが促進される＝8の字のループが大きくなる
* 北端付近では遅くなり、静止衛星に対する追い越しが緩和される＝8の字のループが小さくなる

これによって8の字軌道は南北非対称となる。なお、楕円軌道が極端になると北端で追い越せなくなり、ループがひとつ消えて8の字でなくなる。

<!-- ![analemma_satellite_2.png](https://qiita-image-store.s3.amazonaws.com/0/127130/81770281-8d93-4c97-ddcf-a9b16519cff8.png) -->
![analemma_satellite_2.png](https://qiita-image-store.s3.amazonaws.com/0/127130/3bb3def8-46d7-d1f9-fa74-37fd4a276390.png)

※楕円軌道なので本当は球の中心からの距離が変化するが、図では天球上に投影しているため距離の変化は示していなく、左図では前節の円軌道と重なっている。

### 近点引数[^periapsis]
準天頂衛星なら以上で終わりだが、**楕円軌道が地球に最接近する位置**は南端から変えることができる。その場合も8の字軌道が歪められ、東西にも非対称な形に変化する。前節では8の字でなくなる場合があることを述べたが、実は近点を天の赤道に近づければ8の字が復活するので、アナレンマの形状に意外と重要なパラメーターである。

<!-- ![analemma_satellite_3.png](https://qiita-image-store.s3.amazonaws.com/0/127130/58a6b72c-7682-084b-b8de-b36701660954.png) -->
![analemma_satellite_3.png](https://qiita-image-store.s3.amazonaws.com/0/127130/3bfcf63d-db52-e287-9ae4-435363db9ecc.png)

[^periapsis]: 近点は特に、主星が地球なら近地点、主星が太陽なら近日点という。

### 周期
ところで、ここまでの話は地球の自転周期と同じ静止衛星や準天頂衛星についてだったが、周期が関わる話は一切出てこなかった。なので、他の公転周期をもつ人工衛星であっても、同じ周期で赤道上空を等速円運動する衛星を基準に考えれば8の字軌道に見える。**アナレンマの軌跡を考えるだけなら周期は不要**である。

ある時刻に衛星がどこに位置するかを求めるにはその衛星の公転周期が必要になる。その際はさらに起点とする時刻と位置が1組いる。

また、地上から同じ方角を定期観察することを考えると、アナレンマが昇り沈みする周期も必要になる。これは衛星の公転周期と地球の自転周期から求められる。太陽の場合はこの周期がちょうど1日なので、毎日同じ時刻に観察すればいい。

```math
\frac{1}{T_{\rm{culm}}} = \frac{1}{T_{\rm{rot}}} -\frac{1}{T_{\rm{rev}}}
```

今回は軌跡だけ求められれば満足なので、これらのパラメーターや計算方法にはあまり踏み込まないことにする。

## 計算
衛星のアナレンマについて考えておき、太陽についてはパラメーターを工夫して衛星と見なすことにする。

### 数式
以下のチャートに沿って計算する。最終行の赤い文字 $\varDelta h, \delta$ がアナレンマに必要な変数で、青い文字 $e, i, \omega$ （必要なら $M\_0, T, t\_0$ も）は準備するパラメーター。両矢印の部分は逆算できるので、独立変数はそれらの式の中から1つ選べばいい。

<!-- ![analemma_equations.png](https://qiita-image-store.s3.amazonaws.com/0/127130/d4e2ec9b-b52e-f3a9-827d-4b3ec8aa1216.png) -->
<!-- ![analemma_equations.png](https://qiita-image-store.s3.amazonaws.com/0/127130/f45d6fe5-2b39-f8e0-256e-0ed1ce85473d.png) -->
<!-- ![analemma_equations.png](https://qiita-image-store.s3.amazonaws.com/0/127130/bd7e1bd7-fe1d-9928-491f-3a4e7122511f.png) -->
![analemma_equations.png](https://qiita-image-store.s3.amazonaws.com/0/127130/1aadede5-6e86-ae2b-3331-dc95af70c1c5.png)

* $E$ : 離心近点角
* $e$ : 軌道離心率
* $h$ : 時角
	* $\varDelta h$ : 均時差
* $i$ : 軌道傾斜角
* $M$ : 平均近点角
	* $M_0$ : 元期平均近点角
* $T$ : 公転周期
* $t$ : 時刻
	* $t_0$ : 元期
* $u$ : 緯度引数
* $\alpha$ : 赤経
* $\delta$ : 赤緯
* $\nu$ : 真近点角
* $\pi$ : 円周率 ※数学定数
* $\omega$ : 近点引数

#### 楕円軌道上の位置
チャートの左列上3行は楕円軌道上の位置を計算している。この計算は近点を0°として考える。

* 1行目：単純な等速円運動の計算。ある時刻に平均近点角がいくらかを求めている。時刻にこだわってなくアナレンマの軌跡だけ知りたい場合は飛ばしていい。
* 2行目：[**ケプラーの方程式**](https://en.wikipedia.org/wiki/Kepler's_equation)。これによって同時刻における単純な円運動（ $e = 0$ ）と複雑な楕円運動の位置関係を結び付けられる。 $M$ から $E$ を厳密に求めるのは非常に難しいが、ニュートン法などで数値的に求めるか $e \ll 1$ を利用した近似式を用いれば楽。
* 3行目：楕円の幾何を利用。上で求めた $E$ は楕円を真円に引き伸ばした上で中心から見た角度であり、地球が位置する楕円の焦点から見た実際の角度 $\nu$ に直さなければならない。

$E$ は補助線や円を描かないと現れない角度[^trans-stretch]なので意味をつかみにくいが、これを考えるとケプラーの法則に基づいて幾何学的に位置計算するのが楽になる。

<!-- ![keplerian_motions.png](https://qiita-image-store.s3.amazonaws.com/0/127130/a8f0a6d6-c279-3a5c-5264-bcd0a81e85a0.png) -->
![keplerian_motions.png](https://qiita-image-store.s3.amazonaws.com/0/127130/3b6ac089-174e-4acb-7cb3-2a2847ce2654.png)

[^trans-stretch]: 下図では補助線を引く代わりに平行移動と伸長をしている。

#### 赤道座標系への変換
衛星の軌道は赤道から傾いているので、軌道上で測った角度では地球に対する位置が分かりにくい。そこでチャート左列下3行で、地球の赤道を緯度0°・昇交点を経度0°とする赤道座標系に衛星の位置を変換する。

* 4行目：昇交点から測った角度に修正。
* 5,6行目：昇交点まわりに座標回転。それぞれ赤経と赤緯を求めている。

#### 仮想的な衛星との比較
アナレンマを描くには、基準とする仮想的な衛星の位置も計算する必要がある。

* 右列2～5行目：天の赤道上を等速円運動すると仮定した計算。左列と同じ計算を $e=0$, $i=0^\circ$ （逆行の場合は $i=180^\circ$ ）として実行している。 $\omega$ だけは同じ値を使うことに注意。
* 右列6行目：均時差を計算。両列5行目で求めた赤経の差をとる。本来は時角の差として計算するため、一見符号が逆になっている。

### パラメーターの用意
太陽を衛星に見立てたパラメーターは、太陽を回る地球についてのものから求める。

* 軌道離心率 $e$ ：楕円軌道の形は同じなので、地球の軌道離心率を使えばいい。 → 0.0167
* 軌道傾斜角 $i$ ：太陽の通り道（黄道）が傾いて見えるのは、地球が自転軸を傾けたまま公転しているため。地球の赤道傾斜角を使えばいい。 → 23.44°
* 近点引数 $\omega$ ：太陽から見た地球の楕円軌道と地球から見た太陽の楕円軌道は180°逆なので、地球の近日点引数[^lan]に180°足す。 → 283°

なお、どのパラメーターも少しずつ時間変化しているので、気になるようであれば時間変化項を含むデータから値を計算する必要がある。

* [暦Wiki/平均軌道要素 - 国立天文台](http://eco.mtk.nao.ac.jp/koyomi/wiki/CABFB6D1B5B0C6BBCDD7C1C7.html#g79269c4)
	* 瞬時の平均春分点および黄道にもとづく座標系における平均軌道要素#地球－月重心 のeとπ
* [Axial tilt - Wikipedia](https://en.wikipedia.org/wiki/Axial_tilt)

[^lan]: 近日点引数でなく近日点黄経（$\varpi$）が書いてある場合、春分点の移動量が無視できることを確認すれば代用できる。

### コード
ここでは $\bar{u}$ を等間隔にとって計算してみる。時間間隔一定となるが同時刻ではないことに注意。

```ruby:analemma.rb
include Math

# reduced latitude --> geocentric latitude
# * atan(ratio * tan(theta))
def r2gc(ratio, theta)
	if ratio < 0
		ratio = -ratio
		theta = -theta
	end

	c, s = cos(theta), sin(theta)
	theta + atan2((ratio - 1) * s * c, c * c + ratio * s * s)
end

# mean anomaly --> eccentric anomaly
# * inverse Kepler equation
# * fixed-point iteration
# * e << 1
def ma2ea(e, ma)
	3.times.inject(ma) { |ea,_| ma + e * sin(ea) } # +O(e^4)
end

# mean anomaly <-- eccentric anomaly
# * Kepler's equation
def ea2ma(e, ea)
	ea - e * sin(ea)
end

# eccentric anomaly --> true anomaly
def ea2ta(e, ea)
	r2gc(sqrt((1+e)/(1-e)), ea / 2.0) * 2
end

# eccentric anomaly <-- true anomaly
def ta2ea(e, ta)
	r2gc(sqrt((1-e)/(1+e)), ta / 2.0) * 2
end


# orbital parameters
eccentricity = 0.0167
inclination = 23.44 * (PI / 180)
argument_of_periapsis = 283 * (PI / 180)

# prograde or retrograde
direction = cos(inclination) >= 0 ? 1 : -1

puts "# x[deg] EOT[min] dec[deg]"

0.step(to: 360, by: 15) do |x_deg|
	arg_lat_bar = x_deg * (PI / 180)

	mean_anomaly = arg_lat_bar - argument_of_periapsis
	eccentric_anomaly = ma2ea(eccentricity, mean_anomaly)
	true_anomaly = ea2ta(eccentricity, eccentric_anomaly)
	argument_of_latitude = true_anomaly + argument_of_periapsis

	right_ascension = r2gc(cos(inclination), argument_of_latitude)
	declination = asin(sin(inclination) * sin(argument_of_latitude))
	equation_of_time = arg_lat_bar * direction - right_ascension

	printf "%8.3f %8.3f %8.3f\n", x_deg, equation_of_time * (720 / PI), declination * (180 / PI)
end
```

※赤経に関連する値は、角度の単位を時分秒で表す。1周＝2πラジアン＝360°00′00″＝24<sup>h</sup>00<sup>m</sup>00<sup>s</sup>

### グラフ
横軸に均時差 $\varDelta h$、縦軸に赤緯 $\delta$ をとる。この場合は地上から見上げた形のアナレンマになる。

<!-- ![analemma_for_earth.png](https://qiita-image-store.s3.amazonaws.com/0/127130/d4f7c143-eddf-a5c0-9bc0-9c9a35ea737a.png) -->
![analemma_from_earth.png](https://qiita-image-store.s3.amazonaws.com/0/127130/40eeac3e-fb9d-3504-8db4-bee84ad806d6.png)

軌跡の向きはグラフに書いていないが、「8の字ができる理由」を思い出せば判断できる。天球上の太陽の動きは西から東[^e2w]（グラフでは右から左）で、視太陽は北端や南端付近で平均太陽を追い越すので、グラフの軌跡は上下端で左向きということになる。向きが分かれば、例えば春分における均時差はマイナスであることをグラフから読み取れる。

[^e2w]: 地上から見ると太陽が東から西に動くのは、地球の自転のほうが速いから。

## 応用
### 他の惑星におけるアナレンマ
地球を他の惑星に置き換えれば、その惑星から太陽を観察し続けた場合のアナレンマが得られる。例として、[火星に移住](https://ja.wikipedia.org/wiki/火星の植民)したつもりになって火星から見た太陽のアナレンマを求めてみる。

しかし、やってみようとしてパラメーターを探したら、**その惑星にとっての春分点**から近日点までの角度が見つからなかった。そもそも春分点の位置が分からない。地球の場合は春分点を基準に他のパラメーターを測るので意識していなかった。

<details>
<summary>（追記）火星用カレンダーからの取得</summary>
<div>

火星は暦が色々と提案されていて、季節（＝1年間の時期）を表すときは、春分の瞬間を0°とした軌道上の角度solar longitude (Ls)を基にすることがある。これは[二十四節気](https://eco.mtk.nao.ac.jp/koyomi/faq/24sekki.html)と同様の方法となっている。こうしたカレンダーに**近日点の時期**が書かれていることがあり、 Ls = 251° という値を直接取得できる。
http://www-mars.lmd.jussieu.fr/mars/time/solar_longitude.html
（追記おわり）

</div>
</details>

春分点は公転軌道と赤道面の交わる場所である。公転軌道を特定するパラメーターはふつう揃っているので、あとは赤道面がどの方向に傾いているかが分かればいい（赤道傾斜角では方向まで分からない）。探したら[各惑星の**天の北極**の座標](https://en.wikipedia.org/wiki/Poles_of_astronomical_bodies)が手に入ったので、それから春分点を割り出し、アナレンマに必要な近点引数を求めることにする。

<!--
|   | 天の北極<br>(赤経,赤緯) | 昇交点黄経 | 黄道に対する<br>軌道傾斜角 | 近日点黄経 | → | 惑星にとっての<br>近日点引数 | 赤道傾斜角 |
|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|
| 水星   | 281.00,  61.41 |  48.309 | 7.004 |  77.484 |   | 230.46 |   0.04 |
| 金星   |  92.76, -67.16 |  76.631 | 3.395 | 131.565 |   | 253.82 | 177.36 |
| 地球   | 359.89,  89.90 | 174.831 | 0.002 | 102.995 |   | 103.24 |  23.44 |
| 火星   | 317.66,  52.88 |  49.506 | 1.848 | 336.139 |   |  71.12 |  25.19 |
| 木星   | 268.06,  64.50 | 100.497 | 1.303 |  14.369 |   | 237.14 |   3.12 |
| 土星   |  40.58,  83.54 | 113.620 | 2.489 |  93.157 |   |  99.56 |  26.73 |
| 天王星 |  77.31,  15.18 |  74.020 | 0.773 | 173.021 |   | 185.48 |  97.77 |
| 海王星 | 299.45,  42.95 | 131.781 | 1.770 |  48.126 |   | 182.47 |  28.34 |

（単位は全て°。左4項目の座標系はJ2000.0に基づく。理科年表に合わせ、天の北極は2017年9月4.0日TTの値を算出。金星と天王星は自転の向きを考慮して北極を選択。）

海王星を除き、赤道傾斜角が理科年表と一致しているので計算は問題ないと思う。
-->

```ruby
include Math

def rot(x1, x2, theta)
	theta *= PI / 180
	c, s = cos(theta), sin(theta)
	[x1 * c - x2 * s, x1 * s + x2 * c]
end

OBLIQUITY_J2000 = 8.4381406e4 / 3600
# https://nssdc.gsfc.nasa.gov/planetary/planetfact.html
np_right_ascension, np_declination, inclination, lon_ascending_node, lon_perihelion =
#	  0.000,  90.000, 0.00005, -11.26064, 102.94719 # Earth
	317.681,  52.887, 1.85061,  49.57854, 336.04084 # Mars

# unit vector to the north celestial pole
x, y, z = 1, 0, 0
x, z = rot(x, z, np_declination)
x, y = rot(x, y, np_right_ascension)

# coordinate transformations
y, z = rot(y, z, -OBLIQUITY_J2000)
x, y = rot(x, y, -lon_ascending_node)
y, z = rot(y, z, -inclination)

lon_vernal_equinox = atan2(y, x) * (180 / PI) - 90 + lon_ascending_node
obliquity = acos(z) * (180 / PI)

puts "arg_perihelion = #{(lon_perihelion - lon_vernal_equinox) % 360}"
puts "obliquity      = #{obliquity}"   # verification
```

これで、火星における春分点から近日点までの角度は約71°と求まる。軌道離心率0.0934、赤道傾斜角25.19°、近点引数71°+180°より、アナレンマは以下のようになる。

![analemma_for_mars.png](https://qiita-image-store.s3.amazonaws.com/0/127130/3a5e3f7a-6efb-c401-cb5d-9614bd3146d8.png)

※横軸を時間とみる場合、火星の1日を24時間としたときの分であり、地球より約 2.7% 長い。詳細は[Timekeeping on Mars](https://en.wikipedia.org/wiki/Timekeeping_on_Mars)を参照。

火星は軌道離心率が大きいことが影響し、アナレンマが8の字にならず**涙滴型**になる。春分と秋分では太陽の位置が80分ほど異なり、冬至前後には均時差が大きく変化するため、季節に気をつけながら時間を考えないと暮らしに支障が出るかもしれない。

### 月のアナレンマ
人工衛星のアナレンマが求められるのなら、自然衛星である月についても求められそうである。しかし月は公転軌道が10年程度の周期で複雑に回転しているため、いくつか困難がある。

* 軌道要素の値が直接は書いてない。[Wikipedia英語版](https://en.wikipedia.org/wiki/Moon)だと、昇交点黄経と近地点引数の部分に周期だけ書いてある。
* 1ヶ月後にはパラメーターが1°以上も変わっている。本記事のコードではパラメーターが刻々と変化する影響を反映できない。

前者については、[軌道要素は**時間の多項式**として入手できる](http://eco.mtk.nao.ac.jp/koyomi/wiki/CABFB6D1B5B0C6BBCDD7C1C7.html#g72209bd)ので時刻を指定して求めればいい。ただし、月の軌道要素は黄道面を基準としたものなので、地球の赤道面基準に直す必要がある。この計算には[球面三角法](https://ja.wikipedia.org/wiki/球面三角法)が便利。

![analemma_moon.png](https://qiita-image-store.s3.amazonaws.com/0/127130/06961d60-ec0f-5ce6-5adb-8715cc3a5f50.png)

* 月のパラメーター
	* 角B：軌道傾斜角 $i$
	* 弧AB：昇交点黄経 $\Omega$ （-19.34°/年）
	* 弧BD：近地点引数 $\omega$ 
	* 弧AB+BD：近地点黄経 $\varpi$ （+40.69°/年）
* 地球のパラメーター
	* 角A：黄道傾斜角 $\varepsilon$
* アナレンマに必要なパラメーター（未知）
	* 角Cの補角（鋭角側）：軌道傾斜角 $i'$
	* 弧CD：近地点引数 $\omega'$

```ruby
require 'date'
include Math

# spherical trigonometry
def st_ABc2abC(_A, _B, _c)
	_A *= PI / 180
	_B *= PI / 180
	_c *= PI / 180
	_a = atan2(sin(_A) * sin(_c), cos(_A) * sin(_B) + sin(_A) * cos(_B) * cos(_c)) * (180 / PI)
	_b = atan2(sin(_B) * sin(_c), sin(_A) * cos(_B) + cos(_A) * sin(_B) * cos(_c)) * (180 / PI)
	_C = acos(- cos(_A) * cos(_B) + sin(_A) * sin(_B) * cos(_c)) * (180 / PI)
	[_a, _b, _C]
end

# calc c0[deg] + c1[sec] * t + c2[sec] * t^2 + ...
def calc_polynomial(t, *coefs)
	coef0 = coefs.shift
	coefs.each.with_index(1).inject(coef0) { |sum,(c,i)| sum + c/3600.0 * t**i % 360 } % 360
end


datetime = DateTime.parse(ARGV[0] || '2005-06-20T19:00:00-07:00') + (ARGV[1] || 0).to_i * 1.03505
t = (datetime - DateTime.parse('2000-01-01T12:00:00Z')) / 36525.0

# https://en.wikipedia.org/wiki/Axial_tilt
obliquity          = calc_polynomial(t,  23.43927944, -46.836769, -0.0001831, 0.00200340, -5.76e-7, -4.34e-8)
# http://eco.mtk.nao.ac.jp/koyomi/wiki/CABFB6D1B5B0C6BBCDD7C1C7.html#g72209bd
inclination        = calc_polynomial(t,   5.15668983,  -0.00008)
lon_perigee        = calc_polynomial(t,  83.35324312, 14648449.0869, -37.1582, -0.044970,  0.00018948)
lon_ascending_node = calc_polynomial(t, 125.04455501, -6962890.5431,   7.4722,  0.007702, -0.00005939)

_a, _b, _C = st_ABc2abC(obliquity, inclination, lon_ascending_node)
puts "inclination = #{180 - _C}"
puts "arg_perigee = #{(lon_perigee - lon_ascending_node + _a) % 360}"
```

後者については元々本記事の範囲外なので、代わりに1ヶ月後のパラメーターでもアナレンマを描いて変化の度合いを見てみることにする。月が1ヶ月かけてアナレンマをなぞる際に実はアナレンマ自体が変形していっている、と考える。

実際のアナレンマを再現できるか試してみる。[このページ](https://apod.nasa.gov/apod/ap050713.html)の写真は2005年の夏至頃にニューメキシコ州あたりで東の空を撮影したものと思われる。

| 日時 | 黄道傾斜角 $\varepsilon$ | 軌道傾斜角 $i$ | 近地点黄経 $\varpi$ | 昇交点黄経 $\Omega$ | → | 軌道傾斜角 $i'$ | 近地点引数 $\omega'$ |
|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|
| 2005/06/20 19:00:00 (UTC-7) | 23.44 | 5.157 | 305.89 | 19.27 |   | 28.35 | 302.67 |
| 2005/07/20 19:23:41 (UTC-7) | 23.44 | 5.157 | 309.23 | 17.68 |   | 28.39 | 306.27 |
<!--
| 2018/03/21 15:00:00 (UTC+9) | 23.44 | 5.157 | 104.60 | 132.70 |   | 20.28 | 94.40 |
| 2018/04/20 15:23:41 (UTC+9) | 23.44 | 5.157 | 107.95 | 131.11 |   | 20.40 | 97.55 |
-->

これと軌道離心率 0.0555 を使うと、アナレンマは以下のようになる。

![analemma_of_moon.png](https://qiita-image-store.s3.amazonaws.com/0/127130/12f56c61-0ea4-8db2-23f5-d42377b0cba6.png)

※横軸を時間とみる場合、[平均太陰日](http://eco.mtk.nao.ac.jp/koyomi/wiki/B7EEA4CEBDD0C6FEA4EAA4C8C6EEC3E62FCABFB6D1C2C0B1A2C6FC.html)を24時間としたときの分であり、平均太陽日より約 3.5% 長い。

上側のループが小さい様子が写真とよく似ている。1ヶ月後のアナレンマ（点線）は、ずれが月1個分以内に収まっていて、慎重に観察しない限りずれに気付かなそうである。

<!--
軌道傾斜角 $i'$ は $\varepsilon \pm i$ の間で変化する。傾斜角が小さく近地点が南北端に近いときは、軌道離心率の影響により8の字でなくなる。
-->

## 付録
### 英単語表
英語の説明文を読むために。見慣れない単語がよくあるが、語源を知ると覚えやすいものもある。

| 英語 | 日本語 | 備考 |
|:--|:--|:--|
| [analemma](https://www.etymonline.com/word/analemma) | アナレンマ | 元は "日時計の台座" を表す言葉<br>語源はギリシア語で "支え" |
| anomaly<br>\* eccentric —<br>\* mean —<br>\* true — | 近点角、近点離角<br>\* 離心—<br>\* 平均—<br>\* 真— |   |
| celestial<br>\* — equator<br>\* — sphere<br>\* north/south — pole | 天の<br>\* 天の赤道<br>\* 天球<br>\* 天の北極/南極 |   |
| declination | 赤緯 |   |
| eccentricity | 離心率 | 楕円など円錐曲線の形状を表す数値<br>"orbital —" として軌道要素であることを明示 |
| ecliptic<br>\* — latitude<br>\* — longitude | 黄道、黄道面、食の<br>\* 黄緯<br>\* 黄経 | 日食（solar eclipse）などはこの面上で起きる |
| epoch | 元期 |   |
| equation of time | 均時差 | この "equation" は "差" という意味合い |
| equator | 赤道 |   |
| equinox<br>\* autumnal —<br>\* vernal — | 分点<br>\* 秋分点、秋分<br>\* 春分点、春分 | "nox" はラテン語で "夜"<br>季節名もラテン語由来 |
| hour angle | 時角 |   |
| inclination | 軌道傾斜角 |   |
| Kepler's equation | ケプラーの方程式 |   |
| node<br>\* ascending —<br>\* descending — | 交点<br>\* 昇交点<br>\* 降交点 |   |
| nutation | 章動 |   |
| obliquity | 赤道傾斜角、黄道傾斜角 | 自転軸の傾き（axial tilt）と同じ |
| periapsis, pericenter<br>\* perigee<br>\* perihelion | 近点<br>\* 近地点<br>\* 近日点 | "argument of —" で "—引数" という軌道要素<br>遠点の類は "peri-" → "apo-" または "ap-" |
| precession<br>\* apsidal —, perihelion — | 歳差<br>\* 近点移動、近日点移動 | 重力起因の遅く継続的なパラメーター変化<br>自転軸のすりこぎ運動（"赤道の歳差"）が有名 |
| right ascension | 赤経 | 略称は "RA"<br>天の赤道から直角（right angle）に昇る位置 |
| solar time<br>\* apparent —, true —<br>\* mean — | 太陽時<br>\* 視—、真—<br>\* 平均— |   |
| solstice<br>\* summer —<br>\* winter — | 至点<br>\* 夏至点、夏至<br>\* 冬至点、冬至 | ラテン語で太陽が静止するという意味 |
| year<br>\* anomalistic —<br>\* sidereal —<br>\* solar —, tropical — | 年<br>\* 近点年<br>\* 恒星年<br>\* 太陽年、回帰年 |   |
|   |   |   |

### ケプラーの方程式の数値解法
方程式 $M = E - e \sin{E}$ について、与えられた $M$ に対して $E$ を求めたい。

#### 不動点反復法
方程式を変形すると $E = M + e \sin{E}$ であり、この右辺の関数は $E$ が不動点になっている。以下のように反復計算すると不動点に近づいていく。

```math
\begin{align}
E_0 &= M \\
E_{n+1} &= M + e \sin{E_n}
\end{align}
```

実際に繰り返し代入してみると以下の通り。 $e \ll 1$ を利用して高次の項を打ち切ると近似式が得られる。

```math
\begin{align}
E_0 &= M \\
E_1 &= M + e \sin{M} \\
E_2 &= M + e \sin{\left( M + e \sin{M} \right)} \\
    &\simeq M + e \sin{M} + \frac{1}{2} e^2 \sin{2M} \\
\end{align}
```

#### ニュートン法
方程式の解を求めるのによく使われる方法。こちらのほうが収束が速い。

```math
\begin{align}
f(E) &= E - e \sin{E} - M = 0 \\
\frac{df}{dE} &= 1 - e \cos{E} \\
E_{n+1} &= E_n - \frac{E_n - e \sin{E_n} - M}{1 - e \cos{E_n}}
\end{align}
```

### tanの計算
$\tan \theta' = p \tan \theta \;\; (p > 0)$ という形の計算が何度か出てくるが、 `atan` で単純に計算すると $\theta$ が連続的でも $\theta'$ が不連続になる。これだと逆算したら値が変わったり、均時差が12時間ずれたりして多少面倒なことになる。

簡単な解決方法は、周回数に関する情報を別に管理すること。事前に $\theta$ から $n\pi$ を引いておいて、後でその分を $\theta'$ に足す。

```math
\begin{align}
\tan (\theta' - n\pi)
&= p \tan (\theta - n\pi) \quad (|\theta - n\pi| \le \pi/2) \\

\theta'
&= \arctan \left( p \tan (\theta - n\pi) \right) + n\pi
\end{align}
```

別の方法は、 $\theta' - \theta$ を計算するよう式変形すること。上の方法で境界値を正しく扱う自信が無かったので、コードではこちらを使った。

```math
\begin{align}
\theta'
&= \arctan \left( p \tan \theta \right) \\
&= \arctan \left( p \tan \theta \right) - \arctan \left( \tan \theta \right) + \theta \\
&= \theta + \arctan \frac{(p-1) \tan \theta}{1 + p \tan^2 \theta} \\
&= \theta + \arctan \frac{(p-1) \sin \theta \cos \theta}{\cos^2 \theta + p \sin^2 \theta}
\end{align}
```

## 参考

* 説明
	* [太陽の描く曲線、アナレンマ - こよみのページ](http://koyomi8.com/reki_doc/doc_0553.htm)
	* [Equation of time - Wikipedia](https://en.wikipedia.org/wiki/Equation_of_time)
* パラメーター
	* 理科年表 平成29年版
	* [暦Wiki/平均軌道要素 - 国立天文台](http://eco.mtk.nao.ac.jp/koyomi/wiki/CABFB6D1B5B0C6BBCDD7C1C7.html)
	* [Planetary Fact Sheets](https://nssdc.gsfc.nasa.gov/planetary/planetfact.html)
	* [Axial tilt - Wikipedia](https://en.wikipedia.org/wiki/Axial_tilt)
	* [Poles of astronomical bodies - Wikipedia](https://en.wikipedia.org/wiki/Poles_of_astronomical_bodies)