diff --git a/go.mod b/go.mod index 1a37adf..e860a78 100644 --- a/go.mod +++ b/go.mod @@ -21,6 +21,7 @@ require ( require ( filippo.io/edwards25519 v1.2.0 // indirect + github.com/akhenakh/sgp4 v0.0.0-20260314155803-8ee03fc877eb // indirect github.com/bep/debounce v1.2.1 // indirect github.com/dustin/go-humanize v1.0.1 // indirect github.com/godbus/dbus/v5 v5.1.0 // indirect diff --git a/go.sum b/go.sum index 6705811..1a1d46b 100644 --- a/go.sum +++ b/go.sum @@ -1,5 +1,7 @@ filippo.io/edwards25519 v1.2.0 h1:crnVqOiS4jqYleHd9vaKZ+HKtHfllngJIiOpNpoJsjo= filippo.io/edwards25519 v1.2.0/go.mod h1:xzAOLCNug/yB62zG1bQ8uziwrIqIuxhctzJT18Q77mc= +github.com/akhenakh/sgp4 v0.0.0-20260314155803-8ee03fc877eb h1:d9tZ7tJrssgs7Va9j8iu9vl7BlK2rmIs5RiOU7WQJrs= +github.com/akhenakh/sgp4 v0.0.0-20260314155803-8ee03fc877eb/go.mod h1:JfAepWD223Cel6uRpzYdip/xijWZ2FT457YFLWy8Md4= github.com/bep/debounce v1.2.1 h1:v67fRdBA9UQu2NhLFXrSg0Brw7CexQekrBwDMM8bzeY= github.com/bep/debounce v1.2.1/go.mod h1:H8yggRPQKLUhUoqrJC1bO2xNya7vanpDl7xR3ISbCJ0= github.com/braheezy/shine-mp3 v0.1.0 h1:N2wZhv6ipCFduTSftaPNdDgZ5xFmQAPvB7JcqA4sSi8= diff --git a/internal/sat/sat.go b/internal/sat/sat.go new file mode 100644 index 0000000..0206e14 --- /dev/null +++ b/internal/sat/sat.go @@ -0,0 +1,372 @@ +// Package sat is where a satellite is, where it will be, and what that does to +// a frequency. +// +// Three things live here and nothing else: the orbital elements a station keeps +// (Store), the sky as seen from that station (Track, Passes), and the Doppler +// shift the motion imposes (Shift). The radio, the rotator and the screen are +// all somebody else's business — they are handed numbers by the app layer. +// +// The propagation itself is SGP4 from github.com/akhenakh/sgp4 (Apache-2.0, +// pure Go): the model everyone in this hobby uses, and the one the TLEs are +// built for. Writing it again would be writing it worse. +package sat + +import ( + "fmt" + "math" + "sort" + "strings" + "sync" + "time" + + "github.com/akhenakh/sgp4" +) + +// speedOfLightKmS is the constant every Doppler correction here is built on. +const speedOfLightKmS = 299792.458 + +// Observer is the ground station: where the antenna is, in degrees and metres. +type Observer struct { + Lat, Lon float64 + AltM float64 +} + +// Position is a satellite seen from the ground at one instant. +// +// The two halves answer different questions and both are wanted: where the +// thing IS (for the map) and where to POINT (for the rotator and the Doppler). +type Position struct { + Name string `json:"name"` + At time.Time `json:"at"` + + // Sub-satellite point and height — the map's half. + Lat float64 `json:"lat"` + Lon float64 `json:"lon"` + AltKm float64 `json:"alt_km"` + Footprint float64 `json:"footprint_km"` // radius of the visibility circle + + // Look angles — the station's half. + Az float64 `json:"az"` + El float64 `json:"el"` + RangeKm float64 `json:"range_km"` + RangeRate float64 `json:"range_rate"` // km/s, positive = receding +} + +// Visible reports whether the satellite is above the horizon. +// +// Zero degrees, not a courtesy margin: an operator with a clear take-off works +// a pass from the moment it rises, and a station in a valley knows its own +// horizon better than this package ever will. +func (p Position) Visible() bool { return p.El > 0 } + +// Pass is one crossing of the sky, from rise to set. +type Pass struct { + Name string `json:"name"` + AOS time.Time `json:"aos"` + LOS time.Time `json:"los"` + AOSAz float64 `json:"aos_az"` + LOSAz float64 `json:"los_az"` + MaxEl float64 `json:"max_el"` + MaxElAz float64 `json:"max_el_az"` + MaxElAt time.Time `json:"max_el_at"` + Duration float64 `json:"duration_s"` +} + +// Element is one satellite's orbital elements, as they were published. +// +// The raw lines are kept beside the parsed form because they are what gets +// written to the cache and what an operator pastes in by hand for a bird that +// is not in any feed yet — a freshly launched one, above all, which is exactly +// when everybody wants to hear it. +type Element struct { + Name string `json:"name"` + NORAD int `json:"norad"` + Line1 string `json:"line1"` + Line2 string `json:"line2"` + // Epoch is when these elements were computed. Their accuracy falls away + // from it, which is why the store knows how old they are. + Epoch time.Time `json:"epoch"` + + tle *sgp4.TLE +} + +// Age is how long ago these elements were computed. +func (e Element) Age() time.Duration { + if e.Epoch.IsZero() { + return 0 + } + return time.Since(e.Epoch) +} + +// ParseElement reads one satellite from its two or three TLE lines. +func ParseElement(name, line1, line2 string) (Element, error) { + name = strings.TrimSpace(name) + line1 = strings.TrimSpace(line1) + line2 = strings.TrimSpace(line2) + if line1 == "" || line2 == "" { + return Element{}, fmt.Errorf("sat: %q has no orbital elements", name) + } + raw := line1 + "\n" + line2 + if name != "" { + raw = name + "\n" + raw + } + t, err := sgp4.ParseTLE(raw) + if err != nil { + return Element{}, fmt.Errorf("sat: %q: %w", name, err) + } + if name == "" { + name = strings.TrimSpace(t.Name) + } + return Element{ + Name: name, + NORAD: t.SatelliteNumber, + Line1: line1, + Line2: line2, + Epoch: tleEpoch(t), + tle: t, + }, nil +} + +// tleEpoch turns the two-digit year and fractional day of a TLE into a time. +// +// The pivot is the one the format itself defines: 57 and above is the twentieth +// century, below it the twenty-first. It matters for the AGE of the elements, +// which is how an operator knows whether to trust a prediction. +func tleEpoch(t *sgp4.TLE) time.Time { + if t == nil || t.EpochDay <= 0 { + return time.Time{} + } + year := t.EpochYear + switch { + case year >= 57 && year <= 99: + year += 1900 + case year < 57: + year += 2000 + } + start := time.Date(year, 1, 1, 0, 0, 0, 0, time.UTC) + return start.Add(time.Duration((t.EpochDay - 1) * float64(24*time.Hour))) +} + +// Store holds the elements a station tracks. Safe for concurrent use: the app +// refreshes it from a feed while the tracking loop reads it several times a +// second. +type Store struct { + mu sync.RWMutex + byKey map[string]Element + order []string // insertion order, so a listing reads like the feed + fetch time.Time +} + +func NewStore() *Store { return &Store{byKey: map[string]Element{}} } + +// key is how a satellite is addressed. Case and spacing vary between feeds and +// between the operator's typing; the NORAD number would be exact but is not +// what anybody says out loud. +func key(name string) string { return strings.ToUpper(strings.TrimSpace(name)) } + +// Put adds or replaces one satellite's elements. +func (s *Store) Put(e Element) { + if e.tle == nil || e.Name == "" { + return + } + k := key(e.Name) + s.mu.Lock() + defer s.mu.Unlock() + if _, had := s.byKey[k]; !had { + s.order = append(s.order, k) + } + s.byKey[k] = e +} + +// Get returns one satellite's elements. +func (s *Store) Get(name string) (Element, bool) { + s.mu.RLock() + defer s.mu.RUnlock() + e, ok := s.byKey[key(name)] + return e, ok +} + +// Names lists what the store holds, in the order it arrived. +func (s *Store) Names() []string { + s.mu.RLock() + defer s.mu.RUnlock() + out := make([]string, 0, len(s.order)) + for _, k := range s.order { + out = append(out, s.byKey[k].Name) + } + return out +} + +// Len is how many satellites are known. +func (s *Store) Len() int { + s.mu.RLock() + defer s.mu.RUnlock() + return len(s.byKey) +} + +// FetchedAt is when the elements were last loaded from a feed, zero if never. +func (s *Store) FetchedAt() time.Time { + s.mu.RLock() + defer s.mu.RUnlock() + return s.fetch +} + +// Replace swaps the whole set — what a feed refresh does. The order of the new +// set is kept, and the fetch time is stamped. +func (s *Store) Replace(els []Element, at time.Time) { + byKey := make(map[string]Element, len(els)) + order := make([]string, 0, len(els)) + for _, e := range els { + if e.tle == nil || e.Name == "" { + continue + } + k := key(e.Name) + if _, had := byKey[k]; !had { + order = append(order, k) + } + byKey[k] = e + } + s.mu.Lock() + defer s.mu.Unlock() + s.byKey, s.order, s.fetch = byKey, order, at +} + +// Track is where one satellite is, seen from one station, at one instant. +func (s *Store) Track(name string, obs Observer, at time.Time) (Position, error) { + e, ok := s.Get(name) + if !ok { + return Position{}, fmt.Errorf("sat: %q is not in the element set", name) + } + return e.Track(obs, at) +} + +// Track is the same for elements already in hand. +func (e Element) Track(obs Observer, at time.Time) (Position, error) { + if e.tle == nil { + return Position{}, fmt.Errorf("sat: %q has no usable elements", e.Name) + } + loc := &sgp4.Location{Latitude: obs.Lat, Longitude: obs.Lon, Altitude: obs.AltM} + eci, err := e.tle.FindPositionAtTime(at.UTC()) + if err != nil { + return Position{}, fmt.Errorf("sat: %q: %w", e.Name, err) + } + // The state vector carries the position AND the velocity, which is what the + // look angle needs for the range rate — and the range rate is the whole of + // the Doppler shift. + sv := &sgp4.StateVector{ + X: eci.Position.X, Y: eci.Position.Y, Z: eci.Position.Z, + VX: eci.Velocity.X, VY: eci.Velocity.Y, VZ: eci.Velocity.Z, + } + o, err := sv.GetLookAngle(loc, at.UTC()) + if err != nil { + return Position{}, fmt.Errorf("sat: %q look angle: %w", e.Name, err) + } + return Position{ + Name: e.Name, + At: at.UTC(), + Lat: o.SatellitePos.Latitude, + Lon: o.SatellitePos.Longitude, + AltKm: o.SatellitePos.Altitude, + Footprint: footprintKm(o.SatellitePos.Altitude), + Az: o.LookAngles.Azimuth, + El: o.LookAngles.Elevation, + RangeKm: o.LookAngles.Range, + RangeRate: o.LookAngles.RangeRate, + }, nil +} + +// earthRadiusKm is the mean radius — the footprint is a circle drawn on a +// sphere, and a metre of flattening does not show at that scale. +const earthRadiusKm = 6371.0 + +// footprintKm is the radius of the circle from which the satellite is above the +// horizon: the ground distance to where it sits exactly on it. +func footprintKm(altKm float64) float64 { + if altKm <= 0 { + return 0 + } + return earthRadiusKm * math.Acos(earthRadiusKm/(earthRadiusKm+altKm)) +} + +// Passes lists the crossings of the sky between two instants. +// +// minEl drops the passes not worth waiting for: a bird that scrapes three +// degrees over the horizon is a line in a table that will never be a QSO, and +// on a busy evening those are most of the list. +func (s *Store) Passes(name string, obs Observer, from, to time.Time, minEl float64) ([]Pass, error) { + e, ok := s.Get(name) + if !ok { + return nil, fmt.Errorf("sat: %q is not in the element set", name) + } + if !to.After(from) { + return nil, fmt.Errorf("sat: the window ends before it starts") + } + // Thirty seconds: fine enough that the rise and set times are right to a few + // seconds, coarse enough that a day of predictions for a dozen satellites + // stays instant. + details, err := e.tle.GeneratePasses(obs.Lat, obs.Lon, obs.AltM, from.UTC(), to.UTC(), 30) + if err != nil { + return nil, fmt.Errorf("sat: %q passes: %w", e.Name, err) + } + out := make([]Pass, 0, len(details)) + for _, d := range details { + if d.MaxElevation < minEl { + continue + } + out = append(out, Pass{ + Name: e.Name, + AOS: d.AOS.UTC(), + LOS: d.LOS.UTC(), + AOSAz: d.AOSAzimuth, + LOSAz: d.LOSAzimuth, + MaxEl: d.MaxElevation, + MaxElAz: d.MaxElevationAz, + MaxElAt: d.MaxElevationTime.UTC(), + Duration: d.Duration.Seconds(), + }) + } + return out, nil +} + +// NextPasses is Passes over several satellites at once, in time order — the +// question an operator actually asks: what is coming, and when. +func (s *Store) NextPasses(names []string, obs Observer, from time.Time, window time.Duration, minEl int) []Pass { + var all []Pass + for _, n := range names { + ps, err := s.Passes(n, obs, from, from.Add(window), float64(minEl)) + if err != nil { + continue // a satellite whose elements are missing is simply not listed + } + all = append(all, ps...) + } + sort.Slice(all, func(i, j int) bool { return all[i].AOS.Before(all[j].AOS) }) + return all +} + +// Shift is the Doppler-corrected pair for one moment. +type Shift struct { + DownHz int64 `json:"down_hz"` // where to LISTEN for a nominal downlink + UpHz int64 `json:"up_hz"` // where to TRANSMIT for a nominal uplink +} + +// Doppler corrects a nominal uplink/downlink pair for the satellite's motion. +// +// Two corrections, opposite in sign, and that is the part worth being careful +// about: the DOWNLINK is what we receive, so it arrives shifted and we tune to +// meet it — approaching (negative range rate) means a higher frequency. The +// UPLINK is what the satellite receives, so we must transmit shifted the other +// way for it to land on the transponder's nominal input. +// +// Zero in, zero out: a satellite with no uplink (a beacon) is not given an +// invented one. +func Doppler(p Position, downHz, upHz int64) Shift { + f := -p.RangeRate / speedOfLightKmS // fraction, positive when approaching + var s Shift + if downHz > 0 { + s.DownHz = downHz + int64(math.Round(float64(downHz)*f)) + } + if upHz > 0 { + s.UpHz = upHz - int64(math.Round(float64(upHz)*f)) + } + return s +} diff --git a/internal/sat/sat_test.go b/internal/sat/sat_test.go new file mode 100644 index 0000000..300b73d --- /dev/null +++ b/internal/sat/sat_test.go @@ -0,0 +1,172 @@ +package sat + +import ( + "math" + "testing" + "time" +) + +// A real ISS element set, and the answers a second tracker agrees with. The +// point is not the third decimal — it is that the observer, the epoch and the +// look angle are wired the right way round, which is exactly what silently +// comes out mirrored or an hour late. +const ( + issName = "ISS (ZARYA)" + issLine1 = "1 25544U 98067A 24298.54791435 .00016717 00000+0 30074-3 0 9991" + issLine2 = "2 25544 51.6392 121.4587 0007976 86.1587 27.9639 15.50126585478227" +) + +func issElement(t *testing.T) Element { + t.Helper() + e, err := ParseElement(issName, issLine1, issLine2) + if err != nil { + t.Fatalf("parse: %v", err) + } + return e +} + +func TestElementCarriesItsIdentityAndEpoch(t *testing.T) { + e := issElement(t) + if e.NORAD != 25544 { + t.Errorf("NORAD = %d, want 25544", e.NORAD) + } + // Day 298.548 of 2024 — the day the elements were computed. + want := time.Date(2024, 10, 24, 13, 9, 0, 0, time.UTC) + if d := e.Epoch.Sub(want); d > time.Minute || d < -time.Minute { + t.Errorf("epoch = %s, want about %s", e.Epoch.Format(time.RFC3339), want.Format(time.RFC3339)) + } +} + +// The satellite is somewhere, that somewhere is on Earth's scale, and the look +// angles are self-consistent: a bird below the horizon is further away than one +// overhead, and the footprint is a plausible circle. +func TestTrackIsSaneFromAKnownStation(t *testing.T) { + e := issElement(t) + obs := Observer{Lat: 48.85, Lon: 2.35, AltM: 35} // JN18, Paris + at := time.Date(2024, 10, 24, 14, 0, 0, 0, time.UTC) + + p, err := e.Track(obs, at) + if err != nil { + t.Fatalf("track: %v", err) + } + if p.Lat < -90 || p.Lat > 90 || p.Lon < -180 || p.Lon > 180 { + t.Errorf("sub-satellite point off the planet: %.3f %.3f", p.Lat, p.Lon) + } + if p.AltKm < 300 || p.AltKm > 500 { + t.Errorf("altitude %.1f km — the ISS is not there", p.AltKm) + } + if p.Az < 0 || p.Az >= 360 || p.El < -90 || p.El > 90 { + t.Errorf("look angles out of range: az %.1f el %.1f", p.Az, p.El) + } + // A satellite on the FAR side of the planet is still at a distance — up to + // two Earth radii plus its height — so the useful invariant is the one that + // holds when it is actually up: above the horizon it cannot be further away + // than the slant range to its own footprint edge. + if p.RangeKm < 300 || p.RangeKm > 13200 { + t.Errorf("range %.0f km is not this orbit seen from the ground", p.RangeKm) + } + if p.El > 0 && p.RangeKm > 2600 { + t.Errorf("visible at %.1f° yet %.0f km away", p.El, p.RangeKm) + } + // ~2000 km of visibility circle at 420 km up. + if p.Footprint < 1500 || p.Footprint > 2600 { + t.Errorf("footprint %.0f km", p.Footprint) + } +} + +// Twelve hours of ISS passes over a European station: there are always several, +// they rise before they set, and the filter keeps its promise. +func TestPassesRiseBeforeTheySetAndRespectTheFloor(t *testing.T) { + s := NewStore() + s.Put(issElement(t)) + obs := Observer{Lat: 48.85, Lon: 2.35, AltM: 35} + from := time.Date(2024, 10, 24, 12, 0, 0, 0, time.UTC) + + all, err := s.Passes(issName, obs, from, from.Add(12*time.Hour), 0) + if err != nil { + t.Fatalf("passes: %v", err) + } + if len(all) == 0 { + t.Fatal("no ISS pass in twelve hours over Paris") + } + for _, p := range all { + if !p.LOS.After(p.AOS) { + t.Errorf("%s: sets (%s) before it rises (%s)", p.Name, p.LOS, p.AOS) + } + if p.MaxEl <= 0 || p.MaxEl > 90 { + t.Errorf("max elevation %.1f", p.MaxEl) + } + if p.MaxElAt.Before(p.AOS) || p.MaxElAt.After(p.LOS) { + t.Errorf("the highest point falls outside the pass") + } + } + high, err := s.Passes(issName, obs, from, from.Add(12*time.Hour), 30) + if err != nil { + t.Fatalf("passes: %v", err) + } + if len(high) > len(all) { + t.Error("the elevation floor let MORE passes through") + } + for _, p := range high { + if p.MaxEl < 30 { + t.Errorf("a %.1f° pass survived a 30° floor", p.MaxEl) + } + } +} + +// The two corrections go in OPPOSITE directions, and that is the whole of it: +// the downlink arrives shifted so we tune to meet it, while the uplink has to +// leave shifted the other way to land on the transponder's nominal input. +func TestDopplerCorrectsBothWaysRoundTheRightWay(t *testing.T) { + const down, up = 145_950_000, 435_250_000 + + approaching := Position{RangeRate: -7.0} // km/s, coming towards us + receding := Position{RangeRate: +7.0} + + a := Doppler(approaching, down, up) + if a.DownHz <= down { + t.Errorf("approaching: listen at %d, expected above %d", a.DownHz, down) + } + if a.UpHz >= up { + t.Errorf("approaching: transmit at %d, expected below %d", a.UpHz, up) + } + + r := Doppler(receding, down, up) + if r.DownHz >= down { + t.Errorf("receding: listen at %d, expected below %d", r.DownHz, down) + } + if r.UpHz <= up { + t.Errorf("receding: transmit at %d, expected above %d", r.UpHz, up) + } + + // Size, not just sign: 7 km/s on 145.950 MHz is about 3.4 kHz. + if d := math.Abs(float64(a.DownHz - down)); d < 3000 || d > 3800 { + t.Errorf("shift of %.0f Hz on 2 m at 7 km/s", d) + } + // Stationary is untouched, and an absent uplink is not invented. + if s := Doppler(Position{}, down, 0); s.DownHz != down || s.UpHz != 0 { + t.Errorf("a still satellite was corrected: %+v", s) + } +} + +func TestStoreReplaceKeepsOrderAndStampsTheFetch(t *testing.T) { + s := NewStore() + e := issElement(t) + at := time.Date(2026, 9, 7, 10, 0, 0, 0, time.UTC) + s.Replace([]Element{e}, at) + + if s.Len() != 1 || s.Names()[0] != issName { + t.Errorf("store holds %v", s.Names()) + } + if !s.FetchedAt().Equal(at) { + t.Errorf("fetched at %s", s.FetchedAt()) + } + // Case and spacing vary between feeds and typists; the name is not a + // password. + if _, ok := s.Get("iss (zarya)"); !ok { + t.Error("a satellite could not be found under its own name in another case") + } + if _, err := s.Track("NOTHING", Observer{}, at); err == nil { + t.Error("an unknown satellite was tracked anyway") + } +}