OpsLog/internal/geo/geo.go

// Package geo is the one place that turns Maidenhead locators into positions,
// and positions into distances and bearings.
//
// It exists because there were about to be three copies. These functions lived
// in package main, which the internal packages cannot import, so the PSK
// Reporter watcher had its geometry injected from main and the web publisher
// was about to grow its own. A bearing that disagrees with itself between two
// panels is the kind of fault nobody reports, because each screen looks
// plausible on its own.
package geo

import (
	"math"
	"sort"
	"strings"
)

// GridToLatLon parses a Maidenhead locator (4 or 6 characters) and returns the
// centre of that square in degrees. ok=false on malformed input.
func GridToLatLon(grid string) (lat, lon float64, ok bool) {
	g := strings.ToUpper(strings.TrimSpace(grid))
	if len(g) < 4 {
		return 0, 0, false
	}
	A := g[0] - 'A'
	B := g[1] - 'A'
	C := g[2] - '0'
	D := g[3] - '0'
	if A > 17 || B > 17 || C > 9 || D > 9 {
		return 0, 0, false
	}
	lon = -180 + float64(A)*20 + float64(C)*2
	lat = -90 + float64(B)*10 + float64(D)*1
	if len(g) >= 6 {
		E := g[4] - 'A'
		F := g[5] - 'A'
		if E <= 23 && F <= 23 {
			lon += float64(E)*(5.0/60.0) + 2.5/60.0
			lat += float64(F)*(2.5/60.0) + 1.25/60.0
			return lat, lon, true
		}
	}
	// 4-character locator: aim at the centre of the square.
	lon += 1
	lat += 0.5
	return lat, lon, true
}

// HaversineKm returns the great-circle distance between two positions in
// kilometres. Mean Earth radius 6371 km.
func HaversineKm(lat1, lon1, lat2, lon2 float64) float64 {
	const R = 6371.0
	rad := math.Pi / 180.0
	dLat := (lat2 - lat1) * rad
	dLon := (lon2 - lon1) * rad
	a := math.Sin(dLat/2)*math.Sin(dLat/2) +
		math.Cos(lat1*rad)*math.Cos(lat2*rad)*math.Sin(dLon/2)*math.Sin(dLon/2)
	return R * 2 * math.Atan2(math.Sqrt(a), math.Sqrt(1-a))
}

// DistanceBetweenGrids is the distance in kilometres between two locators,
// ok=false when either cannot be parsed.
func DistanceBetweenGrids(a, b string) (km float64, ok bool) {
	lat1, lon1, ok1 := GridToLatLon(a)
	lat2, lon2, ok2 := GridToLatLon(b)
	if !ok1 || !ok2 {
		return 0, false
	}
	return HaversineKm(lat1, lon1, lat2, lon2), true
}

// LatLonToGrid returns the 4-character Maidenhead square for a position.
func LatLonToGrid(lat, lon float64) string {
	lon = math.Mod(lon+180, 360)
	if lon < 0 {
		lon += 360
	}
	lat = lat + 90
	if lat < 0 {
		lat = 0
	} else if lat > 180 {
		lat = 180
	}
	return string([]byte{
		byte('A' + int(lon/20)),
		byte('A' + int(lat/10)),
		byte('0' + int(math.Mod(lon, 20)/2)),
		byte('0' + int(math.Mod(lat, 10)/1)),
	})
}

// NeighbourGrids returns the square holding (lat, lon) and the ring of squares
// around it — 9 squares for ring 1, 25 for ring 2.
//
// Used to filter the PSK Reporter feed at the BROKER rather than in OpsLog. A
// square is about 111 km tall and 150 km wide at mid latitudes, so one ring is
// roughly the 300 km "around here" the feed already meant, and the traffic that
// used to be received and discarded is never sent.
func NeighbourGrids(lat, lon float64, ring int) []string {
	if ring < 0 {
		ring = 0
	}
	seen := map[string]bool{}
	out := []string{}
	for dLat := -ring; dLat <= ring; dLat++ {
		for dLon := -ring; dLon <= ring; dLon++ {
			// One square step: 1° of latitude, 2° of longitude.
			la := lat + float64(dLat)
			lo := lon + float64(dLon)*2
			if la > 90 || la < -90 {
				continue // past a pole there is no square, not a wrapped one
			}
			g := LatLonToGrid(la, lo)
			if !seen[g] {
				seen[g] = true
				out = append(out, g)
			}
		}
	}
	return out
}

// GridsWithin returns every Maidenhead square whose centre lies within km of
// (lat, lon), nearest first, at most max of them.
//
// NeighbourGrids answers "the ring around here", which is the right shape for a
// few hundred kilometres and the wrong one past that: a ring is a square, so
// asking for 2000 km through it means 1369 squares, most of them further away
// than the ones it left out. This measures instead, and the count then follows
// the AREA asked for rather than the corner of a box.
//
// Nearest first because the caller has to be able to trim: these become one
// broker subscription each, and when there are more than can be afforded, the
// squares to keep are the close ones.
func GridsWithin(lat, lon, km float64, max int) []string {
	if km <= 0 || max <= 0 {
		return nil
	}
	// One square is 1° of latitude and 2° of longitude. Sweep a box big enough
	// to hold the circle — a degree of latitude is ~111 km everywhere, and a
	// degree of longitude never MORE than that, so this cannot cut the circle.
	steps := int(km/111.0) + 1
	type cand struct {
		grid string
		d    float64
	}
	seen := map[string]bool{}
	out := []cand{}
	for dLat := -steps; dLat <= steps; dLat++ {
		for dLon := -2 * steps; dLon <= 2*steps; dLon++ {
			la := lat + float64(dLat)
			lo := lon + float64(dLon)*2
			if la > 90 || la < -90 {
				continue
			}
			g := LatLonToGrid(la, lo)
			if seen[g] {
				continue
			}
			// Measured to the square's own centre, not to the sample point, so
			// two samples landing in one square agree about how far it is.
			cLat, cLon, ok := GridToLatLon(g)
			if !ok {
				continue
			}
			d := HaversineKm(lat, lon, cLat, cLon)
			if d > km {
				continue
			}
			seen[g] = true
			out = append(out, cand{g, d})
		}
	}
	sort.Slice(out, func(i, j int) bool { return out[i].d < out[j].d })
	if len(out) > max {
		out = out[:max]
	}
	grids := make([]string, len(out))
	for i, c := range out {
		grids[i] = c.grid
	}
	return grids
}