first commit
This commit is contained in:
337
backend/propagator/engine.go
Normal file
337
backend/propagator/engine.go
Normal file
@@ -0,0 +1,337 @@
|
||||
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
|
||||
}
|
||||
Reference in New Issue
Block a user