How to calculate Local Sun Time

This is a continuation of Last time. Now also calculate the average sun.

The time difference between the true sun and the average sun, which is the average of the time speeds based on the sun of the year, is called the equation of time. The factors of the equation of time are the elliptical revolution of the earth and the tilt of the earth's axis. Below is the Local Sun Time (LST) program that calculates both the true sun and the average sun.

Preparation

Prepare the orbit information. This time, we will use ALOS-2 (DAICHI-2) TLE (Two Line Element). Save the TLE as a text file. ("./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

Calculation

TLE is calculated by converting it into 6 elements of osculating orbit using skyfield. Simply take the difference between the right ascension of the sun vector and the right ascension ra and convert it to time. In skyfield, it is displayed in degrees, minutes, and seconds like 184deg 09 '18.9" , so I forcibly manipulated the character string.

For the average sun, use equation (2) in A Study on Declination of the Sun and Equation of Time. Did. The Julian century number (JC) from J2000 is calculated from the Julian day, and the right ascension $ α_ {m} $ [h] of the average sun is calculated.

#!/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("At the intersection passing local time(LST True Sun)",true_loc_time)
    print("At the intersection passing local time(LST Mean Sun)",mean_loc_time)

    return

result

The equation of time is about 9 minutes when the equation of time of 2019/9/28 is calculated on this site, so it seems that it is about the same. I understand.

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
At the intersection passing local time(LST True Sun) 00:08
At the intersection passing local time(LST Mean Sun) 23:58

Calculation of satellite orbit LST using python (true sun, average sun)

How to calculate Local Sun Time

Preparation

Calculation

result