fix(sat): the Doppler correction was 250 times too big, and backwards
Reported as "the Doppler moves the frequency enormously", and it did: a 2 m
downlink was being shifted two megahertz across a pass instead of three
kilohertz, in the wrong direction.
The propagator library reports a range rate that is not one. Measured against
the range it is meant to be the derivative of, on this station's own cached
elements:
PO-101 library -1457.3 km/s measured +6.227 km/s
ISS library +2036.8 km/s measured -5.522 km/s
Wrong by a factor of some 250 and of the wrong sign, so the correction both
overshot and pushed the operator away from the station they could hear. Nothing
else was affected — the elevation and the passes come from the look angle, which
is right — which is why this survived: the satellite was in the correct place on
the map while the radio was told to go megahertz away from it.
So OpsLog computes it itself, as the difference between two ranges a second
apart. That cannot be wrong in either magnitude or sign: it differentiates the
very number the panel displays. Two extra propagations per call, which is
microseconds.
Two tests pin it, and both fail against the old behaviour: a range rate faster
than orbital velocity is a units mistake, and the Doppler on the two bands
satellites are worked on has a textbook size — about ±3.5 kHz on 2 m, ±10 kHz on
70 cm.
cmd/satdiag is the throwaway that found it, kept because the next report of this
shape ("the frequency moves oddly", "that pass is not real") is answered by the
same three numbers: which elements the satellite resolved to, the range rate
reported against the range rate measured, and the Doppler each transponder gets.
Co-Authored-By: Claude Opus 5 (1M context) <[email protected]>
This commit is contained in:
Claude Opus 5
@@ -1,4 +1,14 @@
|
||||
[
|
||||
{
|
||||
"version": "0.27.21",
|
||||
"date": "",
|
||||
"en": [
|
||||
"The Doppler correction was wrong — by a factor of about 250, and in the wrong direction. The SGP4 library reports a range rate that is not one: the ISS closing at 5.5 km/s came back as +2036 km/s, which moved a 2 m downlink two megahertz instead of three kilohertz, and moved it the wrong way. OpsLog now measures the range rate from the range itself, which cannot disagree with physics. A 2 m downlink shifts about ±3.5 kHz across a pass and a 70 cm one about ±10 kHz, as they should."
|
||||
],
|
||||
"fr": [
|
||||
"La correction Doppler était fausse — d’un facteur d’environ 250, et dans le mauvais sens. La bibliothèque SGP4 renvoie une vitesse radiale qui n’en est pas une : l’ISS se rapprochant à 5,5 km/s était rapportée à +2036 km/s, ce qui déplaçait une descente 2 m de deux mégahertz au lieu de trois kilohertz, et dans la mauvaise direction. OpsLog mesure désormais cette vitesse à partir de la distance elle-même, ce qui ne peut pas contredire la physique. Une descente 2 m se décale d’environ ±3,5 kHz sur un passage et une 70 cm d’environ ±10 kHz, comme il se doit."
|
||||
]
|
||||
},
|
||||
{
|
||||
"version": "0.27.20",
|
||||
"date": "",
|
||||
|
||||
@@ -0,0 +1,167 @@
|
||||
// Command satdiag answers "is this pass real, and is that Doppler right?" from
|
||||
// a station's own cached elements, without launching OpsLog.
|
||||
//
|
||||
// go run ./cmd/satdiag <data dir> <locator> <satellite>
|
||||
//
|
||||
// It prints which element set the satellite resolved to and how old it is, the
|
||||
// look angle now, the range rate BOTH as the propagator reports it and as the
|
||||
// range actually changes, the Doppler each transponder would be given, and the
|
||||
// next passes. It exists because a wrong Doppler and a wrong satellite look the
|
||||
// same from the front — an operator saying "the frequency moves enormously" —
|
||||
// and the two are told apart by these numbers in a second.
|
||||
//
|
||||
// It found the range rate the SGP4 library reports being wrong by a factor of
|
||||
// 250 and of the wrong sign. Not part of the build.
|
||||
package main
|
||||
|
||||
import (
|
||||
"fmt"
|
||||
"os"
|
||||
"strings"
|
||||
"time"
|
||||
|
||||
"hamlog/internal/sat"
|
||||
)
|
||||
|
||||
func main() {
|
||||
dir := os.Args[1]
|
||||
grid := os.Args[2]
|
||||
name := os.Args[3]
|
||||
|
||||
f := sat.NewFetcher(dir)
|
||||
els, at, err := f.LoadCache()
|
||||
if err != nil {
|
||||
fmt.Println("cache:", err)
|
||||
os.Exit(1)
|
||||
}
|
||||
store := sat.NewStore()
|
||||
store.Replace(els, at)
|
||||
fmt.Printf("elements: %d, fetched %s (%s ago)\n\n", len(els), at.Format(time.RFC3339), time.Since(at).Round(time.Minute))
|
||||
|
||||
birds, err := sat.LoadBirds(dir)
|
||||
if err != nil {
|
||||
fmt.Println("birds:", err)
|
||||
}
|
||||
b, ok := birds.Find(name)
|
||||
if !ok {
|
||||
fmt.Println("no frequency plan for", name)
|
||||
os.Exit(1)
|
||||
}
|
||||
// The same resolution the app does.
|
||||
var el sat.Element
|
||||
found := false
|
||||
if e, ok := store.GetNORAD(b.NORAD); ok {
|
||||
el, found = e, true
|
||||
fmt.Printf("elements found BY NORAD %d → %q\n", b.NORAD, e.Name)
|
||||
} else if e, ok := store.Get(b.Name); ok {
|
||||
el, found = e, true
|
||||
fmt.Printf("elements found by name → %q (NORAD %d)\n", e.Name, e.NORAD)
|
||||
} else {
|
||||
for _, a := range b.Aliases {
|
||||
if e, ok := store.Get(a); ok {
|
||||
el, found = e, true
|
||||
fmt.Printf("elements found by alias %q → %q (NORAD %d)\n", a, e.Name, e.NORAD)
|
||||
break
|
||||
}
|
||||
}
|
||||
}
|
||||
if !found {
|
||||
fmt.Println("NO ELEMENTS")
|
||||
os.Exit(1)
|
||||
}
|
||||
fmt.Printf("epoch: %s (%s old)\n", el.Epoch.Format(time.RFC3339), time.Since(el.Epoch).Round(time.Hour))
|
||||
fmt.Println("line1:", el.Line1)
|
||||
|
||||
lat, lon, okGrid := gridToLatLon(grid)
|
||||
if !okGrid {
|
||||
fmt.Println("bad locator:", grid)
|
||||
os.Exit(1)
|
||||
}
|
||||
obs := sat.Observer{Lat: lat, Lon: lon}
|
||||
fmt.Printf("observer: %s → %.4f, %.4f\n\n", grid, obs.Lat, obs.Lon)
|
||||
|
||||
now := time.Now().UTC()
|
||||
p, err := el.Track(obs, now)
|
||||
if err != nil {
|
||||
fmt.Println("track:", err)
|
||||
os.Exit(1)
|
||||
}
|
||||
fmt.Printf("NOW %s : az %.1f el %.1f range %.0f km\n", now.Format("15:04:05"), p.Az, p.El, p.RangeKm)
|
||||
fmt.Printf(" range rate REPORTED by the library : %+10.3f km/s\n", p.RangeRate)
|
||||
fmt.Printf(" range rate MEASURED (d range / dt) : %+10.3f km/s\n", numericRate(el, obs, now))
|
||||
|
||||
for _, tp := range b.Transponders {
|
||||
sh := sat.Doppler(p, tp.DownLo, tp.UpLo)
|
||||
fmt.Printf(" %-28s down %d → %d (%+d Hz) up %d → %d (%+d Hz)\n",
|
||||
tp.Label, tp.DownLo, sh.DownHz, sh.DownHz-tp.DownLo, tp.UpLo, sh.UpHz, sh.UpHz-tp.UpLo)
|
||||
}
|
||||
|
||||
fmt.Println("\nnext passes (min el 0):")
|
||||
passes, err := store.Passes(el.Name, obs, now, now.Add(12*time.Hour), 0)
|
||||
if err != nil {
|
||||
fmt.Println("passes:", err)
|
||||
}
|
||||
for i, ps := range passes {
|
||||
if i >= 8 {
|
||||
break
|
||||
}
|
||||
fmt.Printf(" %s → %s max %.1f° az %.0f→%.0f\n",
|
||||
ps.AOS.Format("15:04:05"), ps.LOS.Format("15:04:05"), ps.MaxEl, ps.AOSAz, ps.LOSAz)
|
||||
}
|
||||
|
||||
// The extremes of the Doppler across the next pass, which is the honest
|
||||
// answer to "does it move that much".
|
||||
if len(passes) > 0 {
|
||||
ps := passes[0]
|
||||
var lo, hi int64
|
||||
for tt := ps.AOS; tt.Before(ps.LOS); tt = tt.Add(10 * time.Second) {
|
||||
q, err := el.Track(obs, tt)
|
||||
if err != nil {
|
||||
continue
|
||||
}
|
||||
d := sat.Doppler(q, b.Transponders[0].DownLo, 0).DownHz - b.Transponders[0].DownLo
|
||||
if d < lo {
|
||||
lo = d
|
||||
}
|
||||
if d > hi {
|
||||
hi = d
|
||||
}
|
||||
}
|
||||
fmt.Printf("\ndownlink Doppler across that pass: %+d Hz … %+d Hz (span %d Hz)\n", lo, hi, hi-lo)
|
||||
}
|
||||
}
|
||||
|
||||
// gridToLatLon is the six-character Maidenhead centre.
|
||||
func gridToLatLon(g string) (float64, float64, bool) {
|
||||
g = strings.ToUpper(strings.TrimSpace(g))
|
||||
if len(g) < 4 {
|
||||
return 0, 0, false
|
||||
}
|
||||
lon := float64(g[0]-'A')*20 - 180
|
||||
lat := float64(g[1]-'A')*10 - 90
|
||||
lon += float64(g[2]-'0') * 2
|
||||
lat += float64(g[3]-'0') * 1
|
||||
if len(g) >= 6 {
|
||||
lon += float64(g[4]-'A') * (2.0 / 24)
|
||||
lat += float64(g[5]-'A') * (1.0 / 24)
|
||||
lon += (2.0 / 24) / 2
|
||||
lat += (1.0 / 24) / 2
|
||||
} else {
|
||||
lon += 1
|
||||
lat += 0.5
|
||||
}
|
||||
return lat, lon, true
|
||||
}
|
||||
|
||||
// numericRate is the range rate measured rather than reported: the distance a
|
||||
// second later minus the distance a second earlier, over two seconds. It cannot
|
||||
// disagree with physics, so it is the reference the library's own figure is
|
||||
// checked against.
|
||||
func numericRate(el sat.Element, obs sat.Observer, at time.Time) float64 {
|
||||
a, e1 := el.Track(obs, at.Add(-time.Second))
|
||||
b, e2 := el.Track(obs, at.Add(time.Second))
|
||||
if e1 != nil || e2 != nil {
|
||||
return 0
|
||||
}
|
||||
return (b.RangeKm - a.RangeKm) / 2
|
||||
}
|
||||
+37
-4
@@ -271,9 +271,6 @@ func (e Element) Track(obs Observer, at time.Time) (Position, error) {
|
||||
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,
|
||||
@@ -292,10 +289,46 @@ func (e Element) Track(obs Observer, at time.Time) (Position, error) {
|
||||
Az: o.LookAngles.Azimuth,
|
||||
El: o.LookAngles.Elevation,
|
||||
RangeKm: o.LookAngles.Range,
|
||||
RangeRate: o.LookAngles.RangeRate,
|
||||
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
|
||||
|
||||
@@ -4,6 +4,8 @@ import (
|
||||
"math"
|
||||
"testing"
|
||||
"time"
|
||||
|
||||
"github.com/akhenakh/sgp4"
|
||||
)
|
||||
|
||||
// A real ISS element set, and the answers a second tracker agrees with. The
|
||||
@@ -16,6 +18,9 @@ const (
|
||||
issLine2 = "2 25544 51.6392 121.4587 0007976 86.1587 27.9639 15.50126585478227"
|
||||
)
|
||||
|
||||
// testLoc is the same observer, in the form the internal range helper takes.
|
||||
var testLoc = sgp4.Location{Latitude: 48.5, Longitude: 3.0}
|
||||
|
||||
func issElement(t *testing.T) Element {
|
||||
t.Helper()
|
||||
e, err := ParseElement(issName, issLine1, issLine2)
|
||||
@@ -170,3 +175,71 @@ func TestStoreReplaceKeepsOrderAndStampsTheFetch(t *testing.T) {
|
||||
t.Error("an unknown satellite was tracked anyway")
|
||||
}
|
||||
}
|
||||
|
||||
// The range rate is the whole of the Doppler shift, and it was wrong in both
|
||||
// magnitude and sign — the propagator library reported +2036 km/s for an ISS
|
||||
// that was closing at 5.5, which moved the correction hundreds of kilohertz the
|
||||
// wrong way. These are the two things about it that cannot be argued with.
|
||||
func TestRangeRateIsPhysical(t *testing.T) {
|
||||
e := issElement(t)
|
||||
obs := Observer{Lat: 48.5, Lon: 3.0}
|
||||
// A day's worth, sampled across every geometry a pass goes through.
|
||||
base := e.Epoch.Add(2 * time.Hour)
|
||||
for i := 0; i < 240; i++ {
|
||||
at := base.Add(time.Duration(i) * 6 * time.Minute)
|
||||
p, err := e.Track(obs, at)
|
||||
if err != nil {
|
||||
t.Fatalf("track: %v", err)
|
||||
}
|
||||
// Nothing in low earth orbit closes faster than it flies, and it flies
|
||||
// at about 7.7 km/s. A figure outside this is a units mistake.
|
||||
if math.Abs(p.RangeRate) > 8 {
|
||||
t.Fatalf("%s: range rate %.1f km/s — faster than orbital velocity", at.Format(time.RFC3339), p.RangeRate)
|
||||
}
|
||||
// And it must be the derivative of the range we display, sign included.
|
||||
before, _ := e.rangeAt(&testLoc, at.Add(-2*time.Second))
|
||||
after, _ := e.rangeAt(&testLoc, at.Add(2*time.Second))
|
||||
want := (after - before) / 4
|
||||
if math.Abs(p.RangeRate-want) > 0.05 {
|
||||
t.Errorf("%s: range rate %.3f but the range moves at %.3f km/s",
|
||||
at.Format(time.RFC3339), p.RangeRate, want)
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// The Doppler that comes out of it, on the two bands satellites are worked on.
|
||||
// A LEO gives about ±3.5 kHz on 2 m and ±10 kHz on 70 cm; ten times either is
|
||||
// the bug this pins.
|
||||
func TestDopplerStaysWithinTheTextbookRange(t *testing.T) {
|
||||
e := issElement(t)
|
||||
obs := Observer{Lat: 48.5, Lon: 3.0}
|
||||
base := e.Epoch.Add(2 * time.Hour)
|
||||
var maxVHF, maxUHF int64
|
||||
for i := 0; i < 480; i++ {
|
||||
p, err := e.Track(obs, base.Add(time.Duration(i)*3*time.Minute))
|
||||
if err != nil {
|
||||
continue
|
||||
}
|
||||
vhf := Doppler(p, 145_800_000, 0).DownHz - 145_800_000
|
||||
uhf := Doppler(p, 437_800_000, 0).DownHz - 437_800_000
|
||||
if a := abs64(vhf); a > maxVHF {
|
||||
maxVHF = a
|
||||
}
|
||||
if a := abs64(uhf); a > maxUHF {
|
||||
maxUHF = a
|
||||
}
|
||||
}
|
||||
if maxVHF < 1_500 || maxVHF > 5_000 {
|
||||
t.Errorf("2 m Doppler peaks at %d Hz, expected roughly 3.5 kHz", maxVHF)
|
||||
}
|
||||
if maxUHF < 5_000 || maxUHF > 14_000 {
|
||||
t.Errorf("70 cm Doppler peaks at %d Hz, expected roughly 10 kHz", maxUHF)
|
||||
}
|
||||
}
|
||||
|
||||
func abs64(v int64) int64 {
|
||||
if v < 0 {
|
||||
return -v
|
||||
}
|
||||
return v
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user