1// Package geo resolves ZIP codes to coordinates and does great-circle math.2package geo34import (5 "context"6 "encoding/json"7 "fmt"8 "math"9 "net/http"10 "strconv"1112 "airwaves/internal/web"13)1415const earthRadiusKm = 6371.01617// Point is a WGS84 coordinate in decimal degrees.18type Point struct {19 Lat float64 `json:"lat"`20 Lon float64 `json:"lon"`21}2223// Place is a resolved postal code.24type Place struct {25 ZIP string `json:"zip"`26 Name string `json:"name"`27 State string `json:"state"`28 Point Point `json:"point"`29}3031// LookupZIP resolves a US ZIP code to its centroid using zippopotam.us.32func LookupZIP(ctx context.Context, c *http.Client, zip string) (Place, error) {33 raw, err := web.Get(ctx, c, "https://api.zippopotam.us/us/"+zip)34 if err != nil {35 return Place{}, fmt.Errorf("lookup zip %s: %w", zip, err)36 }37 var body struct {38 Places []struct {39 Name string `json:"place name"`40 State string `json:"state abbreviation"`41 Lat string `json:"latitude"`42 Lon string `json:"longitude"`43 } `json:"places"`44 }45 if err := json.Unmarshal(raw, &body); err != nil {46 return Place{}, fmt.Errorf("decode zip %s: %w", zip, err)47 }48 if len(body.Places) == 0 {49 return Place{}, fmt.Errorf("zip %s not found", zip)50 }51 p := body.Places[0]52 lat, err := strconv.ParseFloat(p.Lat, 64)53 if err != nil {54 return Place{}, fmt.Errorf("zip %s latitude: %w", zip, err)55 }56 lon, err := strconv.ParseFloat(p.Lon, 64)57 if err != nil {58 return Place{}, fmt.Errorf("zip %s longitude: %w", zip, err)59 }60 return Place{ZIP: zip, Name: p.Name, State: p.State, Point: Point{Lat: lat, Lon: lon}}, nil61}6263func rad(d float64) float64 { return d * math.Pi / 180 }64func deg(r float64) float64 { return r * 180 / math.Pi }6566// DistanceKm returns the great-circle distance between a and b.67func DistanceKm(a, b Point) float64 {68 dLat := rad(b.Lat - a.Lat)69 dLon := rad(b.Lon - a.Lon)70 h := math.Sin(dLat/2)*math.Sin(dLat/2) +71 math.Cos(rad(a.Lat))*math.Cos(rad(b.Lat))*math.Sin(dLon/2)*math.Sin(dLon/2)72 return 2 * earthRadiusKm * math.Asin(math.Min(1, math.Sqrt(h)))73}7475// BearingDeg returns the initial bearing from a to b, clockwise from true76// north.77func BearingDeg(a, b Point) float64 {78 φ1, φ2 := rad(a.Lat), rad(b.Lat)79 Δλ := rad(b.Lon - a.Lon)80 y := math.Sin(Δλ) * math.Cos(φ2)81 x := math.Cos(φ1)*math.Sin(φ2) - math.Sin(φ1)*math.Cos(φ2)*math.Cos(Δλ)82 return math.Mod(deg(math.Atan2(y, x))+360, 360)83}8485// Intermediate returns the point a fraction f of the way from a to b along86// the great circle.87func Intermediate(a, b Point, f float64) Point {88 φ1, λ1 := rad(a.Lat), rad(a.Lon)89 φ2, λ2 := rad(b.Lat), rad(b.Lon)90 δ := DistanceKm(a, b) / earthRadiusKm91 if δ == 0 {92 return a93 }94 sa := math.Sin((1-f)*δ) / math.Sin(δ)95 sb := math.Sin(f*δ) / math.Sin(δ)96 x := sa*math.Cos(φ1)*math.Cos(λ1) + sb*math.Cos(φ2)*math.Cos(λ2)97 y := sa*math.Cos(φ1)*math.Sin(λ1) + sb*math.Cos(φ2)*math.Sin(λ2)98 z := sa*math.Sin(φ1) + sb*math.Sin(φ2)99 return Point{Lat: deg(math.Atan2(z, math.Hypot(x, y))), Lon: deg(math.Atan2(y, x))}100}