pythonを用いた衛星軌道LSTの計算（真太陽、平均太陽）

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