OpsLog/internal/sat/sat.go

// 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
}

// GetNORAD returns one satellite's elements by catalog number.
//
// The exact join, and the only one that stays exact. A name is written
// differently by every party involved — the feed says "RADFXSAT (FOX-1B)", the
// operator says "AO-91", AMSAT's chart says both — and matching on letters and
// digits gets most of them and quietly misses the rest. The catalog number is
// in the TLE itself and is what a frequency plan should carry.
func (s *Store) GetNORAD(n int) (Element, bool) {
	if n <= 0 {
		return Element{}, false
	}
	s.mu.RLock()
	defer s.mu.RUnlock()
	for _, k := range s.order {
		if e := s.byKey[k]; e.NORAD == n {
			return e, true
		}
	}
	return Element{}, false
}

// 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)
	}
	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: e.rangeRate(loc, at.UTC()),
	}, nil
}

// rangeRate is how fast the satellite is closing or opening, in km/s.
//
// MEASURED, not taken from the propagator. The library reports a range rate
// that is wrong by a factor of some 250 AND has the wrong sign — the ISS at
// −5.5 km/s (closing) came back as +2036 km/s — which put the Doppler
// correction hundreds of kilohertz out and moved it the wrong way. The
// difference between two ranges a second apart cannot be wrong in either
// respect: it differentiates the very number the panel displays.
//
// Two extra propagations per call. SGP4 costs microseconds and this runs at
// most a few hundred times a second across every satellite on screen, so the
// price of being right here is not worth optimising away.
func (e Element) rangeRate(loc *sgp4.Location, at time.Time) float64 {
	const dt = time.Second // ±1 s: far below any curvature in the range, far above float noise
	before, ok1 := e.rangeAt(loc, at.Add(-dt))
	after, ok2 := e.rangeAt(loc, at.Add(dt))
	if !ok1 || !ok2 {
		return 0
	}
	return (after - before) / (2 * dt.Seconds())
}

// rangeAt is the distance to the satellite at one instant, in km.
func (e Element) rangeAt(loc *sgp4.Location, at time.Time) (float64, bool) {
	eci, err := e.tle.FindPositionAtTime(at.UTC())
	if err != nil {
		return 0, false
	}
	sv := &sgp4.StateVector{X: eci.Position.X, Y: eci.Position.Y, Z: eci.Position.Z}
	o, err := sv.GetLookAngle(loc, at.UTC())
	if err != nil {
		return 0, false
	}
	return o.LookAngles.Range, true
}

// 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
}