KeiruaProd

Replicating "Heavens above" satellite charts

Heavens above is a popular astronomy website that has predictions about when satellits are visible at a given location, as well as some very cool sky charts in order to help locate the less visible satellites. Here is the page about the international space station (ISS):

Can we replicate this table ? We already have Celestrak satellite data, now we need to exploit them. We need to:

Satellite positions

We already have TLE data. My initial goal was to implement SGP4 (Simplified General Perturbations), the standard to compute a satellite position given its position and speed. How hard can it be ?

Sadly, it turns out it’s pretty hard, too much for the small evening project I was looking for, so this project became an API problem about how to use skyfield rather than being solved with first principles.

Twilight

There are multiple twilights in astronomy, depending on how dark you need the night to be. It is not necessarily possible to start observing when . Civil twilight is enough for our needs. Again, skyfield has something here.

Calculating apparent magnitude of ISS

Skyfield will provide us the rise and set times, then we can compute the magnitude (how bright) of the satellite. I’ve found a couple things that I implemented below:

Unfortunately I have not found a general expression for the brightness, since it relies on a parameter for which I’ve not found many data points.

The code

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))

Sample output

The magnitude is off for now, but the rest of the data roughly matches heavens-above prediction.

observation day    rise        rise °  culminate      culminate °  set         set °
-----------------  --------  --------  -----------  -------------  --------  -------
2021 May 24        21:23:01   10.004   21:26:02           30.0895  21:29:05  9.99208
2021 May 24        23:00:19   10.1176  23:03:25           33.8905  23:06:33  9.93682
2021 May 24        00:37:05   10.059   00:40:26           71.0289  00:43:48  9.90397
2021 May 25        02:15:36   10.0107  02:16:42           11.2632  02:17:49  9.98614
2021 May 25        22:12:38   10.0092  22:15:40           29.8908  22:18:42  9.98924
2021 May 25        23:49:26   10.0087  23:52:49           77.231   23:56:11  9.98246
2021 May 25        01:26:53   10.0086  01:29:16           18.0853  01:31:38  9.99862
2021 May 26        21:24:44   10.0002  21:27:43           28.1116  21:30:41  9.99074
2021 May 26        23:01:39   10.0756  23:04:57           53.9947  23:08:16  9.91694
2021 May 26        00:38:39   10.0326  00:41:35           28.2578  00:44:33  9.81892
2021 May 27        22:13:39   10.0065  22:16:52           40.7469  22:20:04  9.98306
2021 May 27        23:50:27   10.0113  23:53:41           45.081   23:56:54  9.94896
2021 May 28        21:25:27   10.0058  21:28:33           33.4156  21:31:40  9.91402
2021 May 28        23:02:11   10.1223  23:05:31           72.8247  23:08:52  9.94095
2021 May 28        00:40:33   10.0036  00:41:45           11.5205  00:43:00  9.89364
2021 May 29        22:13:46   10.0012  22:17:08           74.5721  22:20:28  9.99829
2021 May 29        23:51:07   10.0078  23:53:33           18.7166  23:55:58  9.96924
2021 May 30        21:25:12   10.0085  21:28:29           51.7951  21:31:47  9.87673
2021 May 30        23:02:09   10.0982  23:05:06           29.6972  23:08:06  9.82239
2021 May 31        22:13:10   10.0006  22:16:24           48.3399  22:19:38  9.97611
2021 Jun 01        21:24:07   10.0184  21:27:28           78.7881  21:30:50  9.75563
2021 Jun 01        23:02:16   10.3176  23:03:40           12.5475  23:05:28  9.10785

See a typo ? You can suggest a modification on Github.