275 lines
8.2 KiB
Go
275 lines
8.2 KiB
Go
package orbitviewer
|
|
|
|
import (
|
|
"fmt"
|
|
"math"
|
|
"time"
|
|
)
|
|
|
|
// PassPoint represents a single point in time during a pass.
|
|
type PassPoint struct {
|
|
Time time.Time `json:"time"`
|
|
TimeFormatted string `json:"time_formatted"`
|
|
Azimuth float64 `json:"azimuth"`
|
|
AzimuthCardinal string `json:"azimuth_cardinal"`
|
|
Elevation float64 `json:"elevation"`
|
|
RangeKm float64 `json:"range_km"`
|
|
SunElevation float64 `json:"sun_elevation"`
|
|
IsSunlit bool `json:"is_sunlit"`
|
|
}
|
|
|
|
// SatellitePass describes a predicted pass over an observer's location.
|
|
type SatellitePass struct {
|
|
StartTimeFormatted string `json:"start_time_formatted"`
|
|
Start PassPoint `json:"start"`
|
|
Max PassPoint `json:"max"`
|
|
End PassPoint `json:"end"`
|
|
DurationSec int `json:"duration_sec"`
|
|
DurationFormatted string `json:"duration_formatted"`
|
|
Visibility string `json:"visibility"` // "visible", "daylight", "eclipsed"
|
|
VisibilityLabel string `json:"visibility_label"`
|
|
EstimatedMag float64 `json:"estimated_mag"`
|
|
Trajectory []PassPoint `json:"trajectory"`
|
|
}
|
|
|
|
// PredictPassesOptions configures the pass prediction query.
|
|
type PredictPassesOptions struct {
|
|
Observer Geodetic `json:"observer"`
|
|
StartTime time.Time `json:"start_time"`
|
|
Days int `json:"days"` // default 10
|
|
MinElevationDeg float64 `json:"min_elevation_deg"` // default 10.0
|
|
VisibleOnly bool `json:"visible_only"` // default false
|
|
StandardMag float64 `json:"standard_mag"`
|
|
}
|
|
|
|
// PredictPasses scans the given time window and returns predicted satellite passes.
|
|
func PredictPasses(model *SGP4Model, opts PredictPassesOptions) ([]SatellitePass, error) {
|
|
if opts.Days <= 0 || opts.Days > 14 {
|
|
opts.Days = 10
|
|
}
|
|
if opts.MinElevationDeg < 0 {
|
|
opts.MinElevationDeg = 10.0
|
|
}
|
|
if opts.StartTime.IsZero() {
|
|
opts.StartTime = time.Now()
|
|
}
|
|
|
|
endTime := opts.StartTime.Add(time.Duration(opts.Days) * 24 * time.Hour)
|
|
step := 45 * time.Second // Coarse search step
|
|
|
|
var passes []SatellitePass
|
|
t := opts.StartTime
|
|
|
|
inPass := false
|
|
var passPoints []PassPoint
|
|
|
|
for t.Before(endTime) {
|
|
posECEF, _, err := model.PositionECEF(t)
|
|
if err != nil {
|
|
// In case of propagation issue, advance
|
|
t = t.Add(step)
|
|
continue
|
|
}
|
|
|
|
look := CalculateLookAngles(opts.Observer, posECEF)
|
|
isAbove := look.Elevation >= opts.MinElevationDeg
|
|
|
|
if isAbove {
|
|
if !inPass {
|
|
// Just entered pass: backtrack and find exact AOS (elevation crosses minElevation)
|
|
inPass = true
|
|
passPoints = nil
|
|
aosPoint := refineAOS(model, opts.Observer, t.Add(-step), t, opts.MinElevationDeg)
|
|
passPoints = append(passPoints, aosPoint)
|
|
}
|
|
jd := JulianDate(t)
|
|
sunEl := ObserverSunElevation(opts.Observer, t)
|
|
sunlit := IsSatelliteIlluminated(posECEF, jd)
|
|
passPoints = append(passPoints, PassPoint{
|
|
Time: t,
|
|
TimeFormatted: t.Format("15:04:05"),
|
|
Azimuth: math.Round(look.Azimuth*10) / 10,
|
|
AzimuthCardinal: CardinalDirection(look.Azimuth),
|
|
Elevation: math.Round(look.Elevation*10) / 10,
|
|
RangeKm: math.Round(look.Range),
|
|
SunElevation: math.Round(sunEl*10) / 10,
|
|
IsSunlit: sunlit,
|
|
})
|
|
// Inside pass, take finer 15s steps to capture peak elevation accurately
|
|
t = t.Add(15 * time.Second)
|
|
} else {
|
|
if inPass {
|
|
// Just exited pass: refine exact LOS (elevation drops below minElevation)
|
|
losPoint := refineLOS(model, opts.Observer, t.Add(-15*time.Second), t, opts.MinElevationDeg)
|
|
passPoints = append(passPoints, losPoint)
|
|
inPass = false
|
|
|
|
// Analyze the completed pass
|
|
if pass, ok := evaluatePass(passPoints, opts); ok {
|
|
if !opts.VisibleOnly || pass.Visibility == "visible" {
|
|
passes = append(passes, pass)
|
|
}
|
|
}
|
|
passPoints = nil
|
|
}
|
|
// Normal coarse stepping
|
|
t = t.Add(step)
|
|
}
|
|
}
|
|
|
|
return passes, nil
|
|
}
|
|
|
|
func evaluatePass(pts []PassPoint, opts PredictPassesOptions) (SatellitePass, bool) {
|
|
if len(pts) < 2 {
|
|
return SatellitePass{}, false
|
|
}
|
|
|
|
start := pts[0]
|
|
end := pts[len(pts)-1]
|
|
|
|
// Find max elevation point
|
|
maxPt := pts[0]
|
|
for _, p := range pts {
|
|
if p.Elevation > maxPt.Elevation {
|
|
maxPt = p
|
|
}
|
|
}
|
|
|
|
durationSec := int(end.Time.Sub(start.Time).Seconds())
|
|
if durationSec <= 30 {
|
|
return SatellitePass{}, false
|
|
}
|
|
|
|
// Classify visibility based on observer sun elevation and satellite illumination
|
|
visibility := "daylight"
|
|
label := "Daylight"
|
|
estimatedMag := 99.0
|
|
|
|
// Sun elevation < -6° means civil twilight / night for observer
|
|
if maxPt.SunElevation < -6.0 {
|
|
if maxPt.IsSunlit {
|
|
visibility = "visible"
|
|
label = "Visible"
|
|
|
|
// Calculate apparent magnitude at max elevation
|
|
jd := JulianDate(maxPt.Time)
|
|
satECEF, _, _ := CalculateSatECEF(opts.Observer, maxPt)
|
|
obsECEF := GeodeticToECEF(opts.Observer)
|
|
estimatedMag = EstimateVisualMagnitude(opts.StandardMag, satECEF, obsECEF, jd)
|
|
} else {
|
|
visibility = "eclipsed"
|
|
label = "In Shadow"
|
|
}
|
|
}
|
|
|
|
durMin := durationSec / 60
|
|
durSec := durationSec % 60
|
|
|
|
return SatellitePass{
|
|
StartTimeFormatted: start.Time.Local().Format("Jan 02, 15:04:05"),
|
|
Start: start,
|
|
Max: maxPt,
|
|
End: end,
|
|
DurationSec: durationSec,
|
|
DurationFormatted: fmt.Sprintf("%dm %ds", durMin, durSec),
|
|
Visibility: visibility,
|
|
VisibilityLabel: label,
|
|
EstimatedMag: estimatedMag,
|
|
Trajectory: pts,
|
|
}, true
|
|
}
|
|
|
|
// CalculateSatECEF reconstitutes sat position from look angles or propagation
|
|
func CalculateSatECEF(obs Geodetic, p PassPoint) (Vector3, Vector3, error) {
|
|
obsECEF := GeodeticToECEF(obs)
|
|
latRad := obs.Latitude * Deg2Rad
|
|
lonRad := obs.Longitude * Deg2Rad
|
|
azRad := p.Azimuth * Deg2Rad
|
|
elRad := p.Elevation * Deg2Rad
|
|
|
|
sinLat := math.Sin(latRad)
|
|
cosLat := math.Cos(latRad)
|
|
sinLon := math.Sin(lonRad)
|
|
cosLon := math.Cos(lonRad)
|
|
|
|
// Topocentric ENU
|
|
east := p.RangeKm * math.Cos(elRad) * math.Sin(azRad)
|
|
north := p.RangeKm * math.Cos(elRad) * math.Cos(azRad)
|
|
up := p.RangeKm * math.Sin(elRad)
|
|
|
|
rx := -sinLon*east - sinLat*cosLon*north + cosLat*cosLon*up
|
|
ry := cosLon*east - sinLat*sinLon*north + cosLat*sinLon*up
|
|
rz := cosLat*north + sinLat*up
|
|
|
|
return Vector3{
|
|
X: obsECEF.X + rx,
|
|
Y: obsECEF.Y + ry,
|
|
Z: obsECEF.Z + rz,
|
|
}, Vector3{}, nil
|
|
}
|
|
|
|
// refineAOS uses binary search to pinpoint when satellite rises above minEl
|
|
func refineAOS(model *SGP4Model, obs Geodetic, tLow, tHigh time.Time, minEl float64) PassPoint {
|
|
for i := 0; i < 6; i++ {
|
|
tMid := tLow.Add(tHigh.Sub(tLow) / 2)
|
|
posECEF, _, err := model.PositionECEF(tMid)
|
|
if err != nil {
|
|
break
|
|
}
|
|
look := CalculateLookAngles(obs, posECEF)
|
|
if look.Elevation >= minEl {
|
|
tHigh = tMid
|
|
} else {
|
|
tLow = tMid
|
|
}
|
|
}
|
|
posECEF, _, _ := model.PositionECEF(tHigh)
|
|
look := CalculateLookAngles(obs, posECEF)
|
|
jd := JulianDate(tHigh)
|
|
sunEl := ObserverSunElevation(obs, tHigh)
|
|
sunlit := IsSatelliteIlluminated(posECEF, jd)
|
|
return PassPoint{
|
|
Time: tHigh,
|
|
TimeFormatted: tHigh.Format("15:04:05"),
|
|
Azimuth: math.Round(look.Azimuth*10) / 10,
|
|
AzimuthCardinal: CardinalDirection(look.Azimuth),
|
|
Elevation: math.Round(look.Elevation*10) / 10,
|
|
RangeKm: math.Round(look.Range),
|
|
SunElevation: math.Round(sunEl*10) / 10,
|
|
IsSunlit: sunlit,
|
|
}
|
|
}
|
|
|
|
// refineLOS uses binary search to pinpoint when satellite drops below minEl
|
|
func refineLOS(model *SGP4Model, obs Geodetic, tLow, tHigh time.Time, minEl float64) PassPoint {
|
|
for i := 0; i < 6; i++ {
|
|
tMid := tLow.Add(tHigh.Sub(tLow) / 2)
|
|
posECEF, _, err := model.PositionECEF(tMid)
|
|
if err != nil {
|
|
break
|
|
}
|
|
look := CalculateLookAngles(obs, posECEF)
|
|
if look.Elevation >= minEl {
|
|
tLow = tMid
|
|
} else {
|
|
tHigh = tMid
|
|
}
|
|
}
|
|
posECEF, _, _ := model.PositionECEF(tLow)
|
|
look := CalculateLookAngles(obs, posECEF)
|
|
jd := JulianDate(tLow)
|
|
sunEl := ObserverSunElevation(obs, tLow)
|
|
sunlit := IsSatelliteIlluminated(posECEF, jd)
|
|
return PassPoint{
|
|
Time: tLow,
|
|
TimeFormatted: tLow.Format("15:04:05"),
|
|
Azimuth: math.Round(look.Azimuth*10) / 10,
|
|
AzimuthCardinal: CardinalDirection(look.Azimuth),
|
|
Elevation: math.Round(look.Elevation*10) / 10,
|
|
RangeKm: math.Round(look.Range),
|
|
SunElevation: math.Round(sunEl*10) / 10,
|
|
IsSunlit: sunlit,
|
|
}
|
|
}
|