## 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:

- find the night and day period − we need to find the twilight times
- find the passes above our head at a given latitude/longitude
- compute the magnitude

## 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:

- https://stackoverflow.com/questions/53806966/calculating-satellite-apparent-magnitude
- https://astronomy.stackexchange.com/questions/28744/calculating-the-apparent-magnitude-of-a-satellite/28765#28765
- https://github.com/skyfielders/python-skyfield/issues/327

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.