feat(sat): the sky engine — elements, look angles, passes, Doppler

The foundation of the satellite branch, and nothing above it yet: where a satellite is (SGP4 from akhenakh/sgp4, Apache-2.0 and pure Go, so the no-cgo rule holds), where it will be (passes with an elevation floor, because a three-degree scrape is a line in a table that will never be a QSO), and what its motion does to a frequency. The Doppler pair is the part worth being careful about: the two corrections go in OPPOSITE directions. The downlink arrives shifted and we tune to meet it; the uplink must LEAVE shifted the other way to land on the transponder's nominal input. A test pins the signs and the size — 7 km/s on 2 m is about 3.4 kHz. Elements keep their raw lines beside the parsed form: that is what the cache stores and what an operator pastes by hand for a bird no feed carries yet, which is exactly when everyone wants to hear it.
2026-09-07 10:36:56 +02:00
parent b0f76a8ba1
commit 2d71351080
4 changed files with 547 additions and 0 deletions
@@ -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
+}
@@ -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")
+	}
+}