Local Sun Timeの計算方法
地球観測衛星でよく軌道面の特徴を表す指標となるLocal Sun Time（LST)は昇交点または降交点における地方時のことです。これには真太陽と、平均太陽があります。
同じ太陽条件で繰り返し観測することが多い地球観測衛星では、太陽同期準回帰軌道をよく使います。太陽同期準回帰軌道では、同一地点に何日かの周期で戻ってくる習性（回帰性）がありますが、太陽同期性により、その地点においては常にこのLSTの太陽時になります。この衛星は、〇時の軌道、というように言ったりします。なお、光学衛星ではLSTによって午前衛星、午後衛星に分かれ（正午はSunGrint効果があるので避ける）、SARレーダ衛星ではDawn-Dusk軌道と言って、常に夜明けか日の出、LSTでいうと6時か18時の軌道を取ることが多いです。
真太陽と、一年の太陽を基にした時刻のスピードを平均化した平均太陽の時差を均時差といいます。均時差の要因は地球の楕円公転と地軸傾斜です。以下は、Local Sun Time（LST)の、真太陽と平均太陽を両方を計算するプログラムです。
準備
軌道情報を用意します。今回もALOS-2（だいち2号）のTLE(Two Line Element,2行要素)を用います。TLEはテキストファイルで保存しておきます。（"./TLE/ALOS-2.tle"）
ALOS-2
1 39766U 14029A 19271.21643628 .00000056 00000-0 14378-4 0 9997
2 39766 97.9225 6.5909 0001371 85.1486 274.9890 14.79472450288774
計算
TLEを、skyfieldを駆使して、接触軌道6要素に変換して計算します。太陽方向ベクトルの赤経と昇交点赤経raの差分をとって、時刻に変換するだけです。
skyfieldだと184deg 09' 18.9"のように度分秒で表示されてしまうので、無理やり文字列操作しました。
平均太陽は、太陽視赤緯と均時差計算に関する一考察にある式(2)を利用しました。ユリウス日からJ2000からのユリウス世紀数(JC)を計算し、平均太陽の赤経$α_{m}$[h]を計算しています。
#!/usr/bin/env python3.3
# -*- coding: utf-8 -*-
import numpy as np
import math
import sys
import datetime
from skyfield.api import Loader, Star, Topos, load, JulianDate
from skyfield.api import EarthSatellite
from skyfield.elementslib import osculating_elements_of
TLE_NAME = './TLE/ALOS-2.tle'
def loc_time_calc(elm):
ra = elm[5] # raan[rad]
## True Sun ###
planets = load('de421.bsp')
earth = planets['earth']
sun = planets['sun']
astrometric = earth.at(elm[0]).observe(sun) #earth->sun vector
true_sun, dec, distance = astrometric.radec()
true_sun = math.radians(_to_10deg(true_sun.dstr(warn=False))) # get angle
true_sun_angle = ra - true_sun
if true_sun_angle * 12/math.pi + 12 < 0:
true_loc_time = _to_str(24 + true_sun_angle * 12 / math.pi + 12)
else:
true_loc_time = _to_str(true_sun_angle * 12/math.pi + 12)
## Mean Sun ###
JD = elm[0].ut1 # Julian date
Tu = (JD-2451545)/36525 # Julian Julian Century, JC
# mean sun right ascension[h]
alpha_m = 18 +(41+50.54841/60)/60 + Tu*8640184.812866/60/60 + (0.093104/60/60)*(Tu**2) - (0.0000062/60/60)*(Tu**3)
alpha_m = alpha_m % 24
mean_sun = (alpha_m/24) * 2 * np.pi
mean_sun_angle = ra - mean_sun
if mean_sun_angle * 12/math.pi + 12 < 0:
mean_loc_time = _to_str(24 + mean_sun_angle * 12 / math.pi + 12)
else:
mean_loc_time = _to_str(mean_sun_angle * 12/math.pi + 12)
return true_loc_time, mean_loc_time
def _to_10deg(val):
spl1 = val.split()
spl2 = spl1[0].split("deg",1)
spl3 = spl1[1].split("'",1)
spl4 = spl1[2].split('"',1)
degrees = (float(spl4[0]) / 3600) + (float(spl3[0]) / 60) + float(spl2[0])
return degrees
def _to_str(hour):
h_str = datetime.datetime.strftime(datetime.datetime.strptime((str(int(hour))), "%H"),"%H")
m_str = datetime.datetime.strftime(datetime.datetime.strptime(str(int((hour-int(hour))*60)), "%M"),"%M")
return h_str + ":" + m_str
def main():
with open(TLE_NAME) as f:
lines = f.readlines()
sat = EarthSatellite(lines[1], lines[2], lines[0])
print(lines[0], sat.epoch.utc_jpl())
pos = sat.at(sat.epoch)
print(sat.epoch)
elm = osculating_elements_of(pos)
i = elm.inclination.degrees
e = elm.eccentricity
a = elm.semi_major_axis.km
omega = elm.argument_of_periapsis.degrees
ra = elm.longitude_of_ascending_node.degrees
M = elm.mean_anomaly.degrees
print(i,e,a,omega,ra,M)
# Osculating Orbit
osc_elm = [0 for i in range(7)]
osc_elm[0] = sat.epoch
osc_elm[1] = a
osc_elm[2] = e
osc_elm[3] = np.radians(i)
osc_elm[4] = np.radians(omega)
osc_elm[5] = np.radians(ra)
osc_elm[6] = np.radians(M)
true_loc_time, mean_loc_time = loc_time_calc(osc_elm)
print("交点通過地方時(LST True Sun)",true_loc_time)
print("交点通過地方時(LST Mean Sun)",mean_loc_time)
return
結果
均時差は、このサイトで2019/9/28の均時差を計算すると約9分強なので、だいたい合ってそうなことがわかります。
ALOS-2
A.D. 2019-Sep-28 05:11:40.0946 UT
<Time tt=2458754.717237021>
97.92987988479214 0.0012760645968402655 7015.220028255803 69.31302830191312 6.32305263767209 290.71632630644746
交点通過地方時(LST True Sun) 00:08
交点通過地方時(LST Mean Sun) 23:58