diff --git a/changelog.json b/changelog.json index 420bf13..47ef53e 100644 --- a/changelog.json +++ b/changelog.json @@ -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": "", diff --git a/cmd/satdiag/main.go b/cmd/satdiag/main.go new file mode 100644 index 0000000..02a79c6 --- /dev/null +++ b/cmd/satdiag/main.go @@ -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 +// +// 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 +} diff --git a/internal/sat/sat.go b/internal/sat/sat.go index fcf8eb3..008ee7c 100644 --- a/internal/sat/sat.go +++ b/internal/sat/sat.go @@ -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 diff --git a/internal/sat/sat_test.go b/internal/sat/sat_test.go index 300b73d..183e5a4 100644 --- a/internal/sat/sat_test.go +++ b/internal/sat/sat_test.go @@ -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 +}