338 lines
8.7 KiB
Go
338 lines
8.7 KiB
Go
package propagator
|
|
|
|
import (
|
|
"fmt"
|
|
"math"
|
|
"sort"
|
|
"sync"
|
|
"time"
|
|
|
|
"SatMaster/backend/tle"
|
|
|
|
"github.com/akhenakh/sgp4"
|
|
)
|
|
|
|
const (
|
|
EarthRadius = 6371.0 // km
|
|
RadToDeg = 180.0 / math.Pi
|
|
)
|
|
|
|
// SatPosition holds the computed position of a satellite at a given moment.
|
|
type SatPosition struct {
|
|
Name string `json:"name"`
|
|
Latitude float64 `json:"lat"`
|
|
Longitude float64 `json:"lon"`
|
|
Altitude float64 `json:"alt"` // km
|
|
Azimuth float64 `json:"az"` // degrees 0=N
|
|
Elevation float64 `json:"el"` // degrees above horizon
|
|
Range float64 `json:"range"` // km from observer
|
|
RangeRate float64 `json:"rangeRate"` // km/s positive=receding
|
|
Footprint float64 `json:"footprint"` // km radius on ground
|
|
}
|
|
|
|
// Pass represents one overhead pass of a satellite.
|
|
type Pass struct {
|
|
SatName string `json:"satName"`
|
|
AOS time.Time `json:"aos"`
|
|
LOS time.Time `json:"los"`
|
|
MaxEl float64 `json:"maxEl"`
|
|
MaxElTime time.Time `json:"maxElTime"`
|
|
AosAz float64 `json:"aosAz"`
|
|
LosAz float64 `json:"losAz"`
|
|
MaxElAz float64 `json:"maxElAz"`
|
|
Duration float64 `json:"duration"` // seconds
|
|
Points []PassPoint `json:"points"`
|
|
}
|
|
|
|
// PassPoint is one az/el sample in a pass track.
|
|
type PassPoint struct {
|
|
Time time.Time `json:"t"`
|
|
Az float64 `json:"az"`
|
|
El float64 `json:"el"`
|
|
}
|
|
|
|
// Observer holds the QTH location.
|
|
type Observer struct {
|
|
Lat float64 // degrees
|
|
Lon float64 // degrees
|
|
Alt float64 // meters above sea level
|
|
}
|
|
|
|
// Engine propagates satellite positions using SGP4.
|
|
type Engine struct {
|
|
mu sync.RWMutex
|
|
observer Observer
|
|
}
|
|
|
|
func NewEngine() *Engine {
|
|
return &Engine{
|
|
observer: Observer{Lat: 45.75, Lon: 4.85, Alt: 200}, // Default: Lyon area
|
|
}
|
|
}
|
|
|
|
func (e *Engine) SetObserver(lat, lon, altM float64) {
|
|
e.mu.Lock()
|
|
defer e.mu.Unlock()
|
|
e.observer = Observer{Lat: lat, Lon: lon, Alt: altM}
|
|
}
|
|
|
|
func (e *Engine) ObserverPosition() Observer {
|
|
e.mu.RLock()
|
|
defer e.mu.RUnlock()
|
|
return e.observer
|
|
}
|
|
|
|
// Position computes the current position of one satellite.
|
|
func (e *Engine) Position(sat *tle.Satellite, t time.Time) *SatPosition {
|
|
e.mu.RLock()
|
|
obs := e.observer
|
|
e.mu.RUnlock()
|
|
return computePosition(sat, obs, t)
|
|
}
|
|
|
|
// AllPositions computes positions for all loaded satellites.
|
|
func (e *Engine) AllPositions(sats []*tle.Satellite, t time.Time) []SatPosition {
|
|
e.mu.RLock()
|
|
obs := e.observer
|
|
e.mu.RUnlock()
|
|
|
|
result := make([]SatPosition, 0, len(sats))
|
|
for _, sat := range sats {
|
|
if pos := computePosition(sat, obs, t); pos != nil {
|
|
result = append(result, *pos)
|
|
}
|
|
}
|
|
return result
|
|
}
|
|
|
|
// ComputePasses predicts upcoming passes over `hours` hours.
|
|
func (e *Engine) ComputePasses(sat *tle.Satellite, start time.Time, hours float64) []Pass {
|
|
e.mu.RLock()
|
|
obs := e.observer
|
|
e.mu.RUnlock()
|
|
|
|
var passes []Pass
|
|
end := start.Add(time.Duration(float64(time.Hour) * hours))
|
|
step := 10 * time.Second
|
|
t := start
|
|
|
|
var inPass bool
|
|
var cur *Pass
|
|
|
|
for t.Before(end) {
|
|
pos := computePosition(sat, obs, t)
|
|
if pos == nil {
|
|
t = t.Add(step)
|
|
continue
|
|
}
|
|
|
|
if pos.Elevation > 0 && !inPass {
|
|
inPass = true
|
|
aosT := bisectAOS(sat, obs, t.Add(-step), t)
|
|
aosPos := computePosition(sat, obs, aosT)
|
|
if aosPos == nil {
|
|
aosPos = pos
|
|
}
|
|
cur = &Pass{
|
|
SatName: sat.Name,
|
|
AOS: aosT,
|
|
AosAz: aosPos.Azimuth,
|
|
Points: []PassPoint{},
|
|
}
|
|
}
|
|
|
|
if inPass && cur != nil {
|
|
cur.Points = append(cur.Points, PassPoint{Time: t, Az: pos.Azimuth, El: pos.Elevation})
|
|
if pos.Elevation > cur.MaxEl {
|
|
cur.MaxEl = pos.Elevation
|
|
cur.MaxElTime = t
|
|
cur.MaxElAz = pos.Azimuth
|
|
}
|
|
}
|
|
|
|
if pos.Elevation <= 0 && inPass {
|
|
inPass = false
|
|
if cur != nil {
|
|
losT := bisectLOS(sat, obs, t.Add(-step), t)
|
|
losPos := computePosition(sat, obs, losT)
|
|
if losPos == nil {
|
|
losPos = pos
|
|
}
|
|
cur.LOS = losT
|
|
cur.LosAz = losPos.Azimuth
|
|
cur.Duration = losT.Sub(cur.AOS).Seconds()
|
|
if cur.MaxEl >= 1.0 {
|
|
passes = append(passes, *cur)
|
|
}
|
|
cur = nil
|
|
}
|
|
}
|
|
|
|
if pos.Elevation < -15 {
|
|
step = 30 * time.Second
|
|
} else {
|
|
step = 10 * time.Second
|
|
}
|
|
t = t.Add(step)
|
|
}
|
|
|
|
sort.Slice(passes, func(i, j int) bool {
|
|
return passes[i].AOS.Before(passes[j].AOS)
|
|
})
|
|
return passes
|
|
}
|
|
|
|
// ─── Core SGP4 computation ────────────────────────────────────────────────────
|
|
|
|
// parsedTLE caches the orbital elements for a satellite TLE.
|
|
// akhenakh/sgp4 ParseTLE accepts the 3-line TLE as a single string.
|
|
func makeTLEString(sat *tle.Satellite) string {
|
|
return fmt.Sprintf("%s\n%s\n%s", sat.Name, sat.TLE1, sat.TLE2)
|
|
}
|
|
|
|
func computePosition(sat *tle.Satellite, obs Observer, t time.Time) *SatPosition {
|
|
// Parse TLE — akhenakh/sgp4 accepts the 3-line block as one string
|
|
tleObj, err := sgp4.ParseTLE(makeTLEString(sat))
|
|
if err != nil {
|
|
return nil
|
|
}
|
|
|
|
// Propagate to time t → ECI state with geodetic fields
|
|
eciState, err := tleObj.FindPositionAtTime(t)
|
|
if err != nil {
|
|
return nil
|
|
}
|
|
|
|
// Convert ECI to geodetic (lat deg, lon deg, alt km)
|
|
lat, lon, altKm := eciState.ToGeodetic()
|
|
|
|
// Build observer location
|
|
location := &sgp4.Location{
|
|
Latitude: obs.Lat, // degrees
|
|
Longitude: obs.Lon, // degrees
|
|
Altitude: obs.Alt, // meters
|
|
}
|
|
|
|
// Build state vector for look angle calculation
|
|
sv := &sgp4.StateVector{
|
|
X: eciState.Position.X,
|
|
Y: eciState.Position.Y,
|
|
Z: eciState.Position.Z,
|
|
VX: eciState.Velocity.X,
|
|
VY: eciState.Velocity.Y,
|
|
VZ: eciState.Velocity.Z,
|
|
}
|
|
|
|
// GetLookAngle uses eciState.DateTime (the propagated time) not t
|
|
observation, err := sv.GetLookAngle(location, eciState.DateTime)
|
|
if err != nil {
|
|
return nil
|
|
}
|
|
|
|
la := observation.LookAngles
|
|
footprint := EarthRadius * math.Acos(EarthRadius/(EarthRadius+altKm))
|
|
|
|
// Compute range rate by finite difference (2s) - library value is unreliable
|
|
rangeRate := finiteRangeRate(tleObj, location, t)
|
|
|
|
return &SatPosition{
|
|
Name: sat.Name,
|
|
Latitude: lat,
|
|
Longitude: lon,
|
|
Altitude: altKm,
|
|
Azimuth: la.Azimuth,
|
|
Elevation: la.Elevation,
|
|
Range: la.Range,
|
|
RangeRate: rangeRate,
|
|
Footprint: footprint,
|
|
}
|
|
}
|
|
|
|
// finiteRangeRate computes range rate (km/s) by finite difference over 2 seconds.
|
|
// Positive = satellite receding, negative = approaching.
|
|
func finiteRangeRate(tleObj *sgp4.TLE, loc *sgp4.Location, t time.Time) float64 {
|
|
dt := 2.0 // seconds
|
|
t2 := t.Add(time.Duration(dt * float64(time.Second)))
|
|
|
|
e1, err1 := tleObj.FindPositionAtTime(t)
|
|
e2, err2 := tleObj.FindPositionAtTime(t2)
|
|
if err1 != nil || err2 != nil {
|
|
return 0
|
|
}
|
|
|
|
sv1 := &sgp4.StateVector{X: e1.Position.X, Y: e1.Position.Y, Z: e1.Position.Z,
|
|
VX: e1.Velocity.X, VY: e1.Velocity.Y, VZ: e1.Velocity.Z}
|
|
sv2 := &sgp4.StateVector{X: e2.Position.X, Y: e2.Position.Y, Z: e2.Position.Z,
|
|
VX: e2.Velocity.X, VY: e2.Velocity.Y, VZ: e2.Velocity.Z}
|
|
|
|
obs1, err1 := sv1.GetLookAngle(loc, e1.DateTime)
|
|
obs2, err2 := sv2.GetLookAngle(loc, e2.DateTime)
|
|
if err1 != nil || err2 != nil {
|
|
return 0
|
|
}
|
|
|
|
return (obs2.LookAngles.Range - obs1.LookAngles.Range) / dt
|
|
}
|
|
|
|
// bisectAOS finds exact AOS time (precision: 1 second).
|
|
func bisectAOS(sat *tle.Satellite, obs Observer, lo, hi time.Time) time.Time {
|
|
for hi.Sub(lo) > time.Second {
|
|
mid := lo.Add(hi.Sub(lo) / 2)
|
|
pos := computePosition(sat, obs, mid)
|
|
if pos != nil && pos.Elevation > 0 {
|
|
hi = mid
|
|
} else {
|
|
lo = mid
|
|
}
|
|
}
|
|
return hi
|
|
}
|
|
|
|
// bisectLOS finds exact LOS time (precision: 1 second).
|
|
func bisectLOS(sat *tle.Satellite, obs Observer, lo, hi time.Time) time.Time {
|
|
for hi.Sub(lo) > time.Second {
|
|
mid := lo.Add(hi.Sub(lo) / 2)
|
|
pos := computePosition(sat, obs, mid)
|
|
if pos != nil && pos.Elevation > 0 {
|
|
lo = mid
|
|
} else {
|
|
hi = mid
|
|
}
|
|
}
|
|
return lo
|
|
}
|
|
|
|
// GroundtrackPoint is a lat/lon position at a given time.
|
|
type GroundtrackPoint struct {
|
|
Lat float64 `json:"lat"`
|
|
Lon float64 `json:"lon"`
|
|
Time time.Time `json:"t"`
|
|
}
|
|
|
|
// ComputeGroundtrack returns the future orbit track for the next `minutes` minutes.
|
|
// Points every 30 seconds. Handles the antimeridian by splitting segments.
|
|
func (e *Engine) ComputeGroundtrack(sat *tle.Satellite, start time.Time, minutes float64) []GroundtrackPoint {
|
|
e.mu.RLock()
|
|
defer e.mu.RUnlock()
|
|
|
|
var points []GroundtrackPoint
|
|
end := start.Add(time.Duration(float64(time.Minute) * minutes))
|
|
step := 30 * time.Second
|
|
|
|
// Parse TLE once outside the loop for efficiency
|
|
tleObj, err := sgp4.ParseTLE(fmt.Sprintf("%s\n%s\n%s", sat.Name, sat.TLE1, sat.TLE2))
|
|
if err != nil {
|
|
return points
|
|
}
|
|
|
|
for t := start; t.Before(end); t = t.Add(step) {
|
|
eciState, err := tleObj.FindPositionAtTime(t)
|
|
if err != nil {
|
|
continue
|
|
}
|
|
lat, lon, _ := eciState.ToGeodetic()
|
|
points = append(points, GroundtrackPoint{Lat: lat, Lon: lon, Time: t})
|
|
}
|
|
return points
|
|
}
|