from skyfield.api import load, EarthSatellite, wgs84
# Load timescale and planetary ephemeris
ts = load.timescale()
eph = load('de421.bsp')
earth = eph['earth']
# Ground station at San Diego
observer = earth + wgs84.latlon(32.8, -117.1)
# Starlink satellite
sat = EarthSatellite(
"1 44713U 19074P 24129.25000000 .00002185 00000+0 14192-3 0 9992",
"2 44713 53.0000 110.0000 0001250 100.0000 260.0000 15.05500000 15",
"STARLINK-3001", ts
)
# Observation times
times = ts.utc(2025, 5, 10, range(60))
# ✅ This is the only valid way to get az/el
topocentric = observer.at(times).observe(sat).apparent()
az, el, distance = topocentric.altaz()
# Output result
for t, a, e, d in zip(times, az.degrees, el.degrees, distance.km):
print(f"{t.utc_strftime('%H:%M')} - Az: {a:.1f}°, El: {e:.1f}°, Dist: {d:.1f} km")
