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

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

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

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

目的

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

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

8の字ができる理由

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

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

analemma_sun.png

この理屈は、準天頂衛星(みちびきなど)の8の字軌道と同じことである3。人工衛星ならパラメーターを変える想像がしやすいので、3つの軌道要素が8の字の形を決める様子を見ていく。

軌道傾斜角

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

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

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

analemma_satellite_1.png

軌道離心率

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

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

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

analemma_satellite_2.png

近点引数4

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

analemma_satellite_3.png

周期

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

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

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

\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

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

楕円軌道上の位置

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

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

赤道座標系への変換

衛星の軌道は赤道から傾いているので、軌道上で測った角度では地球に対する位置が分かりにくい。そこでチャート左列下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°逆なので、地球の近日点引数5に180°足す。 → 283°

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

コード

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

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″=24h00m00s

グラフ

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

analemma_from_earth.png

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

応用

他の惑星におけるアナレンマ

地球を他の惑星に置き換えれば、その惑星から太陽を観察し続けた場合のアナレンマが得られる。例として、火星に移住したつもりになって火星から見た太陽のアナレンマを求めてみる。

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

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

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

※横軸を時間とみる場合、火星の1日を24時間としたときの分であり、地球より約 2.7% 長い。詳細はTimekeeping on Marsを参照。

月のアナレンマ

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

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

前者については、軌道要素は時間の多項式として入手できるので時刻を指定して求めればいい。ただし、月の軌道要素は黄道面を基準としたものなので、地球の赤道面基準に直す必要がある。この計算には球面三角法が便利。

analemma_moon.png

  • 月のパラメーター
    • 角B:軌道傾斜角 $i$
    • 弧AB:昇交点黄経 $\Omega$ (-19.34°/年)
    • 弧BD:近地点引数 $\omega$
    • 弧AB+BD:近地点黄経 $\varpi$ (+40.69°/年)
  • 地球のパラメーター
    • 角A:黄道傾斜角 $\varepsilon$
  • アナレンマに必要なパラメーター(未知)
    • 角Cの補角(鋭角側):軌道傾斜角 $i'$
    • 弧CD:近地点引数 $\omega'$
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ヶ月かけてアナレンマをなぞる際に実はアナレンマ自体が変形していっている、と考える。

実際のアナレンマを再現できるか試してみる。このページの写真は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

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

analemma_of_moon.png

※横軸を時間とみる場合、平均太陰日を24時間としたときの分であり、平均太陽日より約 3.5% 長い。

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

付録

英単語表

英語の説明文を読むために。

英語 日本語 備考
analemma アナレンマ
anomaly
* eccentric —
* mean —
* true —
近点角、近点離角
* 離心—
* 平均—
* 真—
celestial
* — equator
* — sphere
* north/south — pole
天の
* 天の赤道
* 天球
* 天の北極/南極
declination 赤緯
eccentricity 離心率 楕円など円錐曲線の形状を表す数値
"orbital —" として軌道要素であることを明示
ecliptic
* — latitude
* — longitude
黄道、黄道面
* 黄緯
* 黄経
epoch 元期
equation of time 均時差 この "equation" は "差" という意味合い
equator 赤道
equinox
* autumnal —
* vernal —
分点
* 秋分点、秋分
* 春分点、春分
"nox" はラテン語で "夜"
季節名もラテン語由来
hour angle 時角
inclination 軌道傾斜角
Kepler's equation ケプラーの方程式
node
* ascending —
* descending —
交点
* 昇交点
* 降交点
nutation 章動
obliquity 赤道傾斜角、黄道傾斜角 自転軸の傾き(axial tilt)と同じ
periapsis, pericenter
* perigee
* perihelion
近点
* 近地点
* 近日点
"argument of —" で "—引数" という軌道要素
precession
* apsidal —, perihelion —
歳差
* 近点移動、近日点移動
重力起因の遅く継続的なパラメーター変化
自転軸のすりこぎ運動("赤道の歳差")が有名
right ascension 赤経 略称は "RA"
solar time
* apparent —, true —
* mean —
太陽時
* 視—、真—
* 平均—
solstice
* summer —
* winter —
至点
* 夏至点、夏至
* 冬至点、冬至
ラテン語で太陽が静止するという意味
year
* anomalistic —
* sidereal —
* solar —, tropical —

* 近点年
* 恒星年
* 太陽年、回帰年

ケプラーの方程式の数値解法

方程式 $M = E - e \sin{E}$ について、与えられた $M$ に対して $E$ を求めたい。

不動点反復法

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

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

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

\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}

ニュートン法

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

\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'$ に足す。

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

\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}

参考


  1. 逆に言えば、日時計で正確に時刻を知るには均時差の補正が必要である。 

  2. 1年間の地球の公転運動を地球中心に見た状態で、地球からは太陽が天球上の星々の合間をゆっくり進むように見える(→黄道)。混同してはいけないこととして、毎日の太陽の出没は、太陽の周回でなく地球自身の自転によるものと考え続ける。 

  3. 地表に対して8の字と説明されるが、静止衛星を基準に8の字とも考えることができる。 

  4. 近点は特に、主星が地球なら近地点、主星が太陽なら近日点という。 

  5. 近日点引数でなく近日点黄経($\varpi$)が書いてある場合、春分点の移動量が無視できることを確認すれば代用できる。 

  6. 地上から見ると太陽が東から西に動くのは、地球の自転のほうが速いから。 

Sign up for free and join this conversation.
Sign Up
If you already have a Qiita account log in.