KeiruaProd

import datetime as dt
from pytz import timezone
from skyfield import almanac
from skyfield.api import N, W, wgs84, load
from skyfield.api import load, Topos, EarthSatellite
from tabulate import tabulate
import math

# Coordinates of the city of Lyon
LONGITUDE = 4.834605661405483
LATITUDE = 45.763459143320446
place = Topos(LATITUDE, LONGITUDE)

# TLE data for ISS on may 24
ISS_TLE = """ISS (ZARYA)             
1 25544U 98067A   21144.23167620  .00102477  00000-0  18465-2 0  9996
2 25544  51.6424 101.2130 0002848  30.0669  67.4570 15.48960510284886"""

# We want the predictions between may 24th and june 2nd
ts = load.timescale(builtin=True)
t1 = ts.utc(2021, 5, 24)
t2 = ts.utc(2021, 6, 2+1)
zone = timezone('Europe/Paris')

def find_all_twilights(place, t1, t2, eph):
    """
    Returns all the Civil twilights during interval t1, t2
    https://www.odysseymagazine.com/astronomical-twilight/

    note that almanac.sunrise_sunset
    (https://rhodesmill.org/skyfield/almanac.html#sunrise-and-sunset)
    will provide too many false positives so we discard it
    """
    f = almanac.dark_twilight_day(eph, place)
    times, events = almanac.find_discrete(t1, t2, f)

    twilights = []
    for i in range(0, len(events)-8, 8):
        twilights.append((times[i + 4], times[i + 3+8]))
    return twilights


def satellite_from_tle(tle):
    """Creates an EarthSatellite out of a TLE string"""
    ts = load.timescale()
    lines = tle.split("\n")

    return EarthSatellite(lines[1], lines[2], lines[0], ts)

def compute_magnitude():
    # https://astronomy.stackexchange.com/questions/28744/calculating-the-apparent-magnitude-of-a-satellite/28765#28765
    # int_mag =  intrinsic magnitude. visual Magnitude at 1000km away
    # Mag = int_mag - 15 + 5*LOG(Range) - 2.5*LOG(SIN(phase_angle) + (pi-phase_angle)*COS(phase_angle)
    int_mag = -1.8
    phase_angle = 113 * math.pi/180
    distance = 483
    
    term1 = int_mag
    term2 = 5*math.log(distance/1000.0, 10)

    arg3 = math.sin(phase_angle) + (math.pi - phase_angle) * math.cos(phase_angle)
    term3 = - 2.5*math.log(arg3, 10)

    mag = term1 + term2 + term3
    
    return mag


def extract_events(satellite, position, t_sunrise, t_sunset):
    t, events = satellite.find_events(place, t_sunset, t_sunrise, altitude_degrees=10)
    table = []

    difference = satellite - position

    for event_index in range(0, len(events), 3):
        rise = events[event_index]
        t_rise = t[event_index]
        t_culminate = t[event_index+1]
        t_end = t[event_index+2]

        observer_to_sat_topocentric_rise = difference.at(t_rise)
        observer_to_sat_topocentric_culminate = difference.at(t_culminate)
        observer_to_sat_topocentric_end = difference.at(t_end)
        alt_rise, az_rise, dist_rise = observer_to_sat_topocentric_rise.altaz()
        alt_culminate, az_culminate, dist_culminate = observer_to_sat_topocentric_culminate.altaz()
        alt_end, az_end, dist_end = observer_to_sat_topocentric_end.altaz()
        mag = compute_magnitude()

        table.append([
            t_rise.utc_strftime('%Y %b %d'),
            mag,
            t_rise.astimezone(zone).strftime('%H:%M:%S'),
            alt_rise.degrees,
            t_culminate.astimezone(zone).strftime('%H:%M:%S'),
            alt_culminate.degrees,
            t_end.astimezone(zone).strftime('%H:%M:%S'),
            alt_end.degrees
        ])
    return table

satellite = satellite_from_tle(ISS_TLE)

eph = load('de421.bsp')
twilights = find_all_twilights(place, t1, t2, eph)

table = []
for twilight in twilights:
    t_sunset, t_sunrise = twilight
    for e in extract_events(satellite, place, t_sunrise, t_sunset):
        table.append(e)

headers = [
    "observation day",
    "magnitude",
    "rise",
    "rise °",
    "culminate",
    "culminate °",
    "set",
    "set °",
]

# https://www.heavens-above.com/PassSummary.aspx?satid=25544
# https://pypi.org/project/tabulate/
print(tabulate(table, headers=headers))