diff --git a/.github/workflows/skyeye.yaml b/.github/workflows/skyeye.yaml index 84ff5cab..96af5f6a 100644 --- a/.github/workflows/skyeye.yaml +++ b/.github/workflows/skyeye.yaml @@ -125,6 +125,7 @@ jobs: install: | base-devel git + mingw-w64-ucrt-x86_64-cmake mingw-w64-ucrt-x86_64-gcc mingw-w64-ucrt-x86_64-toolchain mingw-w64-ucrt-x86_64-opus diff --git a/.gitignore b/.gitignore index 59d5118c..e0af34d3 100644 --- a/.gitignore +++ b/.gitignore @@ -96,3 +96,8 @@ dev/ .vscode/settings.json *.zip + +# PROJ build artifacts (downloaded source, build, and install trees) +third_party/proj-*/ +third_party/proj-install/ +third_party/proj-*.tar.gz diff --git a/.vscode/launch.json b/.vscode/launch.json new file mode 100644 index 00000000..9dae3d44 --- /dev/null +++ b/.vscode/launch.json @@ -0,0 +1,10 @@ +{ + // Use IntelliSense to learn about possible attributes. + // Hover to view descriptions of existing attributes. + // For more information, visit: https://go.microsoft.com/fwlink/?linkid=830387 + "version": "0.2.0", + "configurations": [ + + + ] +} \ No newline at end of file diff --git a/Dockerfile b/Dockerfile index 63f0aa4f..bfee5cbf 100644 --- a/Dockerfile +++ b/Dockerfile @@ -5,7 +5,9 @@ RUN apt-get update && apt-get install -y \ make \ lsb-release \ gcc \ + pkg-config \ libopus-dev \ + libproj-dev \ libsoxr-dev \ && rm -rf /var/lib/apt/lists/* WORKDIR /skyeye @@ -18,6 +20,7 @@ RUN go mod download -x COPY cmd cmd COPY internal internal COPY pkg pkg +ENV PKG_CONFIG_PATH=/usr/lib/x86_64-linux-gnu/pkgconfig RUN make skyeye RUN make skyeye-scaler @@ -25,6 +28,7 @@ FROM debian:bookworm-slim AS base RUN apt-get update && apt-get install -y \ ca-certificates \ libopus0 \ + libproj25 \ libsoxr0 \ && rm -rf /var/lib/apt/lists/* diff --git a/Makefile b/Makefile index 208532af..f8d685fa 100644 --- a/Makefile +++ b/Makefile @@ -30,6 +30,16 @@ WHISPER_CPP_REPO = https://github.com/dharmab/whisper.cpp.git WHISPER_CPP_VERSION = v1.7.2-windows-fix WHISPER_CPP_BUILD_ENV = +# PROJ (built from source for static linking on Windows) +PROJ_VERSION = 9.4.1 +PROJ_TARBALL = third_party/proj-$(PROJ_VERSION).tar.gz +PROJ_SRC_DIR = third_party/proj-$(PROJ_VERSION) +PROJ_PREFIX = $(abspath third_party/proj-install) +PROJ_LIB = $(PROJ_PREFIX)/lib/libproj.a +PROJ_PC_PATH = $(PROJ_PREFIX)/lib/pkgconfig +PROJ_DEP = +PROJ_CMAKE_FLAGS = + # Compiler variables and flags GOBUILDVARS = GOARCH=$(GOARCH) ABS_WHISPER_CPP_PATH = $(abspath $(WHISPER_CPP_PATH)) @@ -63,11 +73,15 @@ SKYEYE_SCALER_BIN = skyeye-scaler.exe GO = /ucrt64/bin/go GOBUILDVARS += GOROOT="/ucrt64/lib/go" GOPATH="/ucrt64" # Static linking on Windows to avoid MSYS2 dependency at runtime -LIBRARIES = opus soxr -CFLAGS = $(shell pkg-config $(LIBRARIES) --cflags --static) -BUILD_VARS += CFLAGS='$(CFLAGS)' -EXTLDFLAGS = $(shell pkg-config $(LIBRARIES) --libs --static) +LIBRARIES = opus soxr proj +PKG_CONFIG_ENV = PKG_CONFIG_PATH="$(PROJ_PC_PATH):$${PKG_CONFIG_PATH}" +CFLAGS = $(shell $(PKG_CONFIG_ENV) pkg-config $(LIBRARIES) --cflags --static) +BUILD_VARS += CFLAGS='$(CFLAGS)' PKG_CONFIG_PATH="$(PROJ_PC_PATH):$${PKG_CONFIG_PATH}" +EXTLDFLAGS = $(shell $(PKG_CONFIG_ENV) pkg-config $(LIBRARIES) --libs --static) LDFLAGS += -linkmode external -extldflags "$(EXTLDFLAGS) -static" +PROJ_DEP = $(PROJ_LIB) +# MinGW builds of PROJ 9.4.1 miss in filemanager.cpp; force-include to fix. +PROJ_CMAKE_FLAGS += -DCMAKE_CXX_FLAGS="-include cstdint" endif BUILD_VARS += LDFLAGS='$(LDFLAGS)' @@ -84,6 +98,8 @@ install-msys2-dependencies: base-devel \ $(MINGW_PACKAGE_PREFIX)-toolchain \ $(MINGW_PACKAGE_PREFIX)-go \ + $(MINGW_PACKAGE_PREFIX)-cmake \ + $(MINGW_PACKAGE_PREFIX)-sqlite3 \ $(MINGW_PACKAGE_PREFIX)-opus \ $(MINGW_PACKAGE_PREFIX)-libsoxr @@ -93,6 +109,7 @@ install-arch-linux-dependencies: git \ base-devel \ go \ + proj \ opus \ libsoxr @@ -103,6 +120,7 @@ install-debian-dependencies: git \ build-essential \ golang-go \ + libproj-dev \ libopus-dev \ libopus0 \ libsoxr-dev \ @@ -115,6 +133,7 @@ install-fedora-dependencies: development-tools \ c-development \ golang \ + proj-devel \ opus-devel \ opus \ soxr-devel \ @@ -128,6 +147,7 @@ install-macos-dependencies: llvm \ pkg-config \ go \ + proj \ libsoxr \ opus @@ -139,17 +159,47 @@ $(LIBWHISPER_PATH) $(WHISPER_H_PATH): if [ ! -f $(LIBWHISPER_PATH) -o ! -f $(WHISPER_H_PATH) ]; then git -C "$(WHISPER_CPP_PATH)" checkout --quiet $(WHISPER_CPP_VERSION) || git clone --depth 1 --branch $(WHISPER_CPP_VERSION) -c advice.detachedHead=false "$(WHISPER_CPP_REPO)" "$(WHISPER_CPP_PATH)" && $(WHISPER_CPP_BUILD_ENV) make -C $(WHISPER_CPP_PATH)/bindings/go whisper; fi if [ -f third_party/whisper.cpp/whisper.a ] && [ ! -f $(LIBWHISPER_PATH) ]; then cp third_party/whisper.cpp/whisper.a $(LIBWHISPER_PATH); fi +.PHONY: proj +proj: $(PROJ_LIB) + +$(PROJ_TARBALL): + curl -L -o $(PROJ_TARBALL) https://download.osgeo.org/proj/proj-$(PROJ_VERSION).tar.gz + +$(PROJ_LIB): $(PROJ_TARBALL) + rm -rf "$(PROJ_SRC_DIR)" "$(PROJ_PREFIX)" + tar -xzf "$(PROJ_TARBALL)" -C third_party + cmake -S "$(PROJ_SRC_DIR)" -B "$(PROJ_SRC_DIR)/build" \ + -DCMAKE_BUILD_TYPE=Release \ + -DCMAKE_INSTALL_PREFIX="$(PROJ_PREFIX)" \ + -DBUILD_SHARED_LIBS=OFF \ + -DBUILD_APPS=OFF \ + -DBUILD_PROJSYNC=OFF \ + -DBUILD_CS2CS=OFF \ + -DBUILD_GIE=OFF \ + -DBUILD_GEOD=OFF \ + -DBUILD_PROJ=OFF \ + -DBUILD_PROJINFO=OFF \ + -DBUILD_CCT=OFF \ + -DENABLE_CURL=OFF \ + -DENABLE_TIFF=OFF \ + -DENABLE_PROJSYNC=OFF \ + -DENABLE_CCT=OFF \ + -DBUILD_TESTING=OFF \ + -DBUILD_EXAMPLES=OFF \ + $(PROJ_CMAKE_FLAGS) + cmake --build "$(PROJ_SRC_DIR)/build" --target install --config Release + .PHONY: whisper whisper: $(LIBWHISPER_PATH) $(WHISPER_H_PATH) .PHONY: generate -generate: +generate: $(PROJ_DEP) $(BUILD_VARS) $(GO) generate $(BUILD_FLAGS) ./... -$(SKYEYE_BIN): generate $(SKYEYE_SOURCES) $(LIBWHISPER_PATH) $(WHISPER_H_PATH) +$(SKYEYE_BIN): generate $(SKYEYE_SOURCES) $(LIBWHISPER_PATH) $(WHISPER_H_PATH) $(PROJ_DEP) $(BUILD_VARS) $(GO) build $(BUILD_FLAGS) ./cmd/skyeye/ -$(SKYEYE_SCALER_BIN): generate $(SKYEYE_SOURCES) +$(SKYEYE_SCALER_BIN): generate $(SKYEYE_SOURCES) $(PROJ_DEP) $(BUILD_VARS) $(GO) build $(BUILD_FLAGS) ./cmd/skyeye-scaler/ .PHONY: run @@ -187,4 +237,4 @@ mostlyclean: .PHONY: clean clean: mostlyclean - rm -rf "$(WHISPER_CPP_PATH)" + rm -rf "$(WHISPER_CPP_PATH)" "$(PROJ_SRC_DIR)" "$(PROJ_PREFIX)" "$(PROJ_TARBALL)" diff --git a/go.mod b/go.mod index ba3eb3f4..dee1f217 100644 --- a/go.mod +++ b/go.mod @@ -18,6 +18,7 @@ require ( github.com/jba/omap v0.1.0 github.com/lithammer/shortuuid/v3 v3.0.7 github.com/martinlindhe/unit v0.0.0-20230420213220-4adfd7d0a0d6 + github.com/michiho/go-proj/v10 v10.5.11 github.com/nabbl/piper v0.0.0-20240819160100-e51f2288a5c0 github.com/openai/openai-go v0.1.0-alpha.41 github.com/pasztorpisti/go-crc v1.0.0 diff --git a/go.sum b/go.sum index 90e2dd18..2c3b7235 100644 --- a/go.sum +++ b/go.sum @@ -494,6 +494,8 @@ github.com/matttproud/golang_protobuf_extensions v1.0.1 h1:4hp9jkHxhMHkqkrB3Ix0j github.com/matttproud/golang_protobuf_extensions v1.0.1/go.mod h1:D8He9yQNgCq6Z5Ld7szi9bcBfOoFv/3dc6xSMkL2PC0= github.com/mgechev/revive v1.13.0 h1:yFbEVliCVKRXY8UgwEO7EOYNopvjb1BFbmYqm9hZjBM= github.com/mgechev/revive v1.13.0/go.mod h1:efJfeBVCX2JUumNQ7dtOLDja+QKj9mYGgEZA7rt5u+0= +github.com/michiho/go-proj/v10 v10.5.11 h1:CDZYhc19W730k8L1x9ColP0wrjfQ1JV9hpYfC3FFqzY= +github.com/michiho/go-proj/v10 v10.5.11/go.mod h1:eYkLt9XqKpy/r/Y3hGo+rLk6UXA9vNZ6n7vxnuVRIUE= github.com/mitchellh/go-homedir v1.1.0 h1:lukF9ziXFxDFPkA1vsr5zpc1XuPDn/wFntq5mG+4E0Y= github.com/mitchellh/go-homedir v1.1.0/go.mod h1:SfyaCUpYCn1Vlf4IUYiD9fPX4A5wJrkLzIz1N1q0pr0= github.com/mitchellh/mapstructure v1.5.0 h1:jeMsZIYE/09sWLaz43PL7Gy6RuMjD2eJVyuac5Z2hdY= diff --git a/internal/application/app.go b/internal/application/app.go index fefe8f5a..4d59feff 100644 --- a/internal/application/app.go +++ b/internal/application/app.go @@ -23,6 +23,7 @@ import ( "github.com/dharmab/skyeye/pkg/sim" "github.com/dharmab/skyeye/pkg/simpleradio" srs "github.com/dharmab/skyeye/pkg/simpleradio/types" + "github.com/dharmab/skyeye/pkg/spatial" "github.com/dharmab/skyeye/pkg/synthesizer/speakers" "github.com/dharmab/skyeye/pkg/telemetry" "github.com/dharmab/skyeye/pkg/traces" @@ -385,6 +386,9 @@ func (a *Application) updateBullseyes() { log.Warn().Err(err).Msg("error reading bullseye") } else { a.radar.SetBullseye(bullseye, coalition) + if name, changed := spatial.DetectTerrainFromBullseye(coalition.String(), bullseye); changed { + log.Debug().Str("terrain", name).Msg("terrain detected from bullseyes") + } } } } diff --git a/pkg/bearings/true.go b/pkg/bearings/true.go index 52e2718b..52c49ced 100644 --- a/pkg/bearings/true.go +++ b/pkg/bearings/true.go @@ -13,7 +13,7 @@ type True struct { var _ Bearing = True{} -// NewTrueBearing creates a new bearing from the given value. +// NewTrueBearing creates a new bearing from the given value, relative to true north. func NewTrueBearing(value unit.Angle) True { return True{θ: normalize(value)} } diff --git a/pkg/controller/controller.go b/pkg/controller/controller.go index 2c80da22..41edfc09 100644 --- a/pkg/controller/controller.go +++ b/pkg/controller/controller.go @@ -86,6 +86,7 @@ func New( threatMonitoringCooldown time.Duration, threatMonitoringRequiresSRS bool, ) *Controller { + //log.Debug().Msgf("enableAutomaticPicture is %v", enableAutomaticPicture) return &Controller{ coalition: coalition, scope: rdr, diff --git a/pkg/radar/group.go b/pkg/radar/group.go index 1d39c9e6..c09dfdc0 100644 --- a/pkg/radar/group.go +++ b/pkg/radar/group.go @@ -52,6 +52,7 @@ func (g *group) Bullseye() *brevity.Bullseye { } declination, err := bearings.Declination(*g.bullseye, g.missionTime()) + if err != nil { log.Error().Err(err).Stringer("group", g).Msg("failed to get declination for group") } diff --git a/pkg/radar/nearest.go b/pkg/radar/nearest.go index c1f15984..d2068902 100644 --- a/pkg/radar/nearest.go +++ b/pkg/radar/nearest.go @@ -74,8 +74,8 @@ func (r *Radar) FindNearestGroupWithBRAA( if grp == nil { return nil } - - declination := r.Declination(origin) + log.Debug().Any("target latlong", grp).Msgf("target latlong lat %f lon %f", grp.point().Lat(), grp.point().Lon()) + declination := r.Declination(trackfile.LastKnown().Point) bearing := spatial.TrueBearing(origin, grp.point()).Magnetic(declination) _range := spatial.Distance(origin, grp.point()) aspect := brevity.AspectFromAngle(bearing, trackfile.Course()) diff --git a/pkg/radar/nearest_test.go b/pkg/radar/nearest_test.go new file mode 100644 index 00000000..63d676e2 --- /dev/null +++ b/pkg/radar/nearest_test.go @@ -0,0 +1,110 @@ +package radar + +import ( + "context" + "sync" + "testing" + "time" + + "github.com/dharmab/skyeye/pkg/coalitions" + "github.com/dharmab/skyeye/pkg/sim" + "github.com/dharmab/skyeye/pkg/trackfiles" + "github.com/martinlindhe/unit" + "github.com/paulmach/orb" + "github.com/stretchr/testify/assert" +) + +func TestNearest(t *testing.T) { + t.Parallel() + testCases := []struct { + name string + origin orb.Point + pointA orb.Point + pointB orb.Point + expectedID uint64 + }{ + { + name: "finds nearest Red aircraft to origin", + origin: orb.Point{33.405794, 69.047461}, + pointA: orb.Point{33.405794, 69.047461}, + pointB: orb.Point{24.973478, 70.068836}, + expectedID: 2, + }, + } + + for _, test := range testCases { + t.Run(test.name, func(t *testing.T) { + t.Parallel() + starts := make(chan sim.Started, 10) + updates := make(chan sim.Updated, 10) + fades := make(chan sim.Faded, 10) + rdr := New(coalitions.Blue, starts, updates, fades, 20*unit.NauticalMile) + rdr.SetMissionTime(time.Date(1999, 06, 11, 12, 0, 0, 0, time.UTC)) + rdr.SetBullseye(orb.Point{22.867128, 69.047461}, coalitions.Blue) + + // Start the radar to process updates + ctx, cancel := context.WithCancel(context.Background()) + wg := &sync.WaitGroup{} + go rdr.Run(ctx, wg) + + // Add trackfiles to radar via updates channel + updates <- sim.Updated{ + Labels: trackfiles.Labels{ + ID: 1, + ACMIName: "F-15C", + Name: "Eagle 1", + Coalition: coalitions.Blue, + }, + Frame: trackfiles.Frame{ + Time: time.Date(1999, 06, 11, 12, 0, 0, 0, time.UTC), + Point: test.pointA, + Altitude: 30000 * unit.Foot, + AGL: func() *unit.Length { + agl := 30000 * unit.Foot + return &agl + }(), + Heading: 90 * unit.Degree, + }, + } + updates <- sim.Updated{ + Labels: trackfiles.Labels{ + ID: 2, + ACMIName: "F-15C", + Name: "Eagle 2", + Coalition: coalitions.Red, + }, + Frame: trackfiles.Frame{ + Time: time.Date(1999, 06, 11, 12, 0, 0, 0, time.UTC), + Point: test.pointB, + Altitude: 30000 * unit.Foot, + AGL: func() *unit.Length { + agl := 30000 * unit.Foot + return &agl + }(), + Heading: 90 * unit.Degree, + }, + } + + // Wait for updates to be processed + time.Sleep(100 * time.Millisecond) + + group := rdr.FindNearestGroupWithBRAA( + test.origin, + 0*unit.NauticalMile, + 100000*unit.Foot, + 300*unit.NauticalMile, + coalitions.Red, + 0, + ) + + assert.NotNil(t, group) + if group != nil { + ids := group.ObjectIDs() + assert.Contains(t, ids, test.expectedID) + } + + // Clean up + cancel() + }) + } +} diff --git a/pkg/radar/picture.go b/pkg/radar/picture.go index 4d4744f5..ba4cf870 100644 --- a/pkg/radar/picture.go +++ b/pkg/radar/picture.go @@ -23,6 +23,7 @@ func (r *Radar) Picture(radius unit.Length, coalition coalitions.Coalition, filt if spatial.IsZero(origin) { log.Warn().Msg("center point is not set yet, using bullseye") origin = r.Bullseye(coalition) + log.Debug().Any("origin", origin).Msgf("using bullseye point for picture, lat %f, lon %f", origin.Lat(), origin.Lon()) if spatial.IsZero(origin) { log.Warn().Msg("bullseye point is not yet set, picture will be incoherent") } diff --git a/pkg/radar/radar.go b/pkg/radar/radar.go index 145b1c2c..99cbd069 100644 --- a/pkg/radar/radar.go +++ b/pkg/radar/radar.go @@ -170,11 +170,13 @@ func (r *Radar) handleUpdate(update sim.Updated) { trackfile, ok := r.contacts.getByID(update.Labels.ID) if ok { trackfile.Update(update.Frame) - } else { - trackfile = trackfiles.New(update.Labels) - r.contacts.set(trackfile) - logger.Info().Msg("created new trackfile") + return } + + trackfile = trackfiles.New(update.Labels) + trackfile.Update(update.Frame) + r.contacts.set(trackfile) + logger.Info().Msg("created new trackfile") } // handleGarbageCollection removes trackfiles that have not been updated in a long time. @@ -254,6 +256,8 @@ func (r *Radar) Declination(p orb.Point) unit.Angle { r.missionTimeLock.RLock() defer r.missionTimeLock.RUnlock() declination, err := bearings.Declination(p, r.missionTime) + //log.Debug().Any("declination", declination.Degrees()).Msgf("computed magnetic radar declination at point lat %f lon %f", p.Lat(), p.Lon()) + if err != nil { log.Error().Err(err).Msg("failed to get declination") } diff --git a/pkg/recognizer/prompt.go b/pkg/recognizer/prompt.go index 57fec9dc..34e1fff5 100644 --- a/pkg/recognizer/prompt.go +++ b/pkg/recognizer/prompt.go @@ -2,7 +2,7 @@ package recognizer import "fmt" -// prompt constructs a prompt for OpenAI's audio transcription models. See https://platform.openai.com/docs/guides/speech-to-text#prompting +// prompt constructs a prompt for OpenAI's audio transcription models. See https://platform.openai.com/docs/guides/speech-to-text#prompting. func prompt(callsign string) string { return fmt.Sprintf("Either ANYFACE or %s, PILOT CALLSIGN, DIGITS, one of 'RADIO' 'ALPHA' 'BOGEY' 'PICTURE' 'DECLARE' 'SNAPLOCK' 'SPIKED', ARGUMENTS such as BULLSEYE, BRAA, numbers or digits.", callsign) } diff --git a/pkg/spatial/pydcs_bearing.go b/pkg/spatial/pydcs_bearing.go new file mode 100644 index 00000000..cd4098a6 --- /dev/null +++ b/pkg/spatial/pydcs_bearing.go @@ -0,0 +1,629 @@ +package spatial + +import ( + "fmt" + "math" + "sync" + "sync/atomic" + + "github.com/martinlindhe/unit" + "github.com/michiho/go-proj/v10" + "github.com/paulmach/orb" + "github.com/rs/zerolog/log" + + "github.com/dharmab/skyeye/pkg/bearings" +) + +func greatCircleDeg(lat1, lon1, lat2, lon2 float64) float64 { + lat1r := lat1 * math.Pi / 180 + lon1r := lon1 * math.Pi / 180 + lat2r := lat2 * math.Pi / 180 + lon2r := lon2 * math.Pi / 180 + dLat := lat2r - lat1r + dLon := lon2r - lon1r + a := math.Sin(dLat/2)*math.Sin(dLat/2) + math.Cos(lat1r)*math.Cos(lat2r)*math.Sin(dLon/2)*math.Sin(dLon/2) + c := 2 * math.Atan2(math.Sqrt(a), math.Sqrt(1-a)) + return c // radians +} + +// TransverseMercator represents the parameters for a Transverse Mercator projection. +type TransverseMercator struct { + CentralMeridian int + FalseEasting float64 + FalseNorthing float64 + ScaleFactor float64 +} + +type latLonBounds struct { + minLat float64 + maxLat float64 + minLon float64 + maxLon float64 +} + +type terrainDef struct { + name string + tm TransverseMercator + boundsXY [4]float64 // x1, y1, x2, y2 in projected coordinates (DCS x/y) + latLonBox latLonBounds + centerLat float64 + centerLon float64 +} + +func (l latLonBounds) contains(lat, lon float64) bool { + return lat >= l.minLat && lat <= l.maxLat && lon >= l.minLon && lon <= l.maxLon +} + +func (l latLonBounds) area() float64 { + return math.Abs(l.maxLat-l.minLat) * math.Abs(l.maxLon-l.minLon) +} + +var ( + projectionMu sync.RWMutex + currentProjection = CaucasusProjection() + currentTerrain = "Caucasus" + terrainDetected atomic.Bool + bullseyes = make(map[string]orb.Point) // coalition or source label -> bullseye +) + +var terrainDefs = []terrainDef{ + {name: "Afghanistan", tm: AfghanistanProjection(), boundsXY: [4]float64{532000.0, -534000.0, -512000.0, 757000.0}, centerLat: 33.9346, centerLon: 66.24705}, + {name: "Caucasus", tm: CaucasusProjection(), boundsXY: [4]float64{380 * 1000, -560 * 1000, -600 * 1000, 1130 * 1000}, centerLat: 43.69666, centerLon: 32.96}, + {name: "Falklands", tm: FalklandsProjection(), boundsXY: [4]float64{74967, -114995, -129982, 129991}, centerLat: 52.468, centerLon: 59.173}, + {name: "GermanyCW", tm: GermanyColdWarProjection(), boundsXY: [4]float64{260000.0, -1100000.0, -600000.0, -425000.0}, centerLat: 51.0, centerLon: 11.0}, + {name: "Iraq", tm: IraqProjection(), boundsXY: [4]float64{440000.0, -500000.0, -950000.0, 850000.0}, centerLat: 30.76, centerLon: 59.07}, + {name: "Kola", tm: KolaProjection(), boundsXY: [4]float64{-315000, -890000, 900000, 856000}, centerLat: 68.0, centerLon: 22.5}, + {name: "MarianaIslands", tm: MarianasProjection(), boundsXY: [4]float64{1000 * 10000, -1000 * 1000, -300 * 1000, 500 * 1000}, centerLat: 13.485, centerLon: 144.798}, + {name: "Nevada", tm: NevadaProjection(), boundsXY: [4]float64{-167000.0, -330000.0, -500000.0, 210000.0}, centerLat: 39.81806, centerLon: -114.73333}, + {name: "Normandy", tm: NormandyProjection(), boundsXY: [4]float64{-132707.843750, -389942.906250, 185756.156250, 165065.078125}, centerLat: 41.3, centerLon: 0.18}, + {name: "PersianGulf", tm: PersianGulfProjection(), boundsXY: [4]float64{-218768.750000, -392081.937500, 197357.906250, 333129.125000}, centerLat: 0, centerLon: 0}, + {name: "Sinai", tm: SinaiProjection(), boundsXY: [4]float64{-450000, -280000, 500000, 560000}, centerLat: 30.047, centerLon: 31.224}, + {name: "Syria", tm: SyriaProjection(), boundsXY: [4]float64{-320000, -579986, 300000, 579998}, centerLat: 35.021, centerLon: 35.901}, + {name: "TheChannel", tm: TheChannelProjection(), boundsXY: [4]float64{74967, -114995, -129982, 129991}, centerLat: 50.875, centerLon: 1.5875}, +} + +func init() { + for i := range terrainDefs { + if err := computeLatLonBounds(&terrainDefs[i]); err != nil { + log.Warn().Err(err).Str("terrain", terrainDefs[i].name).Msg("failed to compute lat/lon bounds for terrain") + } + } +} + +func terrainDefByName(name string) (terrainDef, bool) { + for _, td := range terrainDefs { + if td.name == name { + return td, true + } + } + return terrainDef{}, false +} + +func computeLatLonBounds(td *terrainDef) error { + // boundsXY are DCS projected coords: x=easting, y=northing in meters. + x1, y1, x2, y2 := td.boundsXY[0], td.boundsXY[1], td.boundsXY[2], td.boundsXY[3] + norths := []float64{y1, y2} + easts := []float64{x1, x2} + + minLat := math.Inf(1) + maxLat := math.Inf(-1) + minLon := math.Inf(1) + maxLon := math.Inf(-1) + + for _, north := range norths { + for _, east := range easts { + lat, lon, err := ProjectionToLatLongFor(td.tm, north, east) + if err != nil { + return fmt.Errorf("convert bounds corner: %w", err) + } + if lat < minLat { + minLat = lat + } + if lat > maxLat { + maxLat = lat + } + if lon < minLon { + minLon = lon + } + if lon > maxLon { + maxLon = lon + } + } + } + + td.latLonBox = latLonBounds{ + minLat: minLat, + maxLat: maxLat, + minLon: minLon, + maxLon: maxLon, + } + if td.centerLat == 0 && td.centerLon == 0 { + td.centerLat = (minLat + maxLat) / 2 + td.centerLon = (minLon + maxLon) / 2 + } + return nil +} + +func bullseyeInsideBounds(td terrainDef, bullseye orb.Point) bool { + if td.latLonBox.contains(bullseye.Lat(), bullseye.Lon()) { + return true + } + + xMin := math.Min(td.boundsXY[0], td.boundsXY[2]) + xMax := math.Max(td.boundsXY[0], td.boundsXY[2]) + yMin := math.Min(td.boundsXY[1], td.boundsXY[3]) + yMax := math.Max(td.boundsXY[1], td.boundsXY[3]) + + x, z, err := LatLongToProjectionFor(td.tm, bullseye.Lat(), bullseye.Lon()) + if err == nil { + north := x + east := z + if east >= xMin && east <= xMax && north >= yMin && north <= yMax { + return true + } + } + + return false +} + +func setCurrentTerrain(name string, tm TransverseMercator) { + projectionMu.Lock() + defer projectionMu.Unlock() + if currentTerrain != name { + log.Debug(). + Str("from", currentTerrain). + Str("to", name). + Msg("switching terrain projection") + } + currentTerrain = name + currentProjection = tm +} + +// ForceTerrain overrides the current terrain selection and disables auto-detection. +func ForceTerrain(name string, tm TransverseMercator) { + setCurrentTerrain(name, tm) + terrainDetected.Store(true) + projectionMu.Lock() + bullseyes = make(map[string]orb.Point) + projectionMu.Unlock() +} + +// ResetTerrainToDefault resets terrain selection to the default (Caucasus) and re-enables auto-detection. +func ResetTerrainToDefault() { + setCurrentTerrain("Caucasus", CaucasusProjection()) + terrainDetected.Store(false) + projectionMu.Lock() + bullseyes = make(map[string]orb.Point) + projectionMu.Unlock() +} + +func getCurrentProjection() TransverseMercator { + projectionMu.RLock() + defer projectionMu.RUnlock() + return currentProjection +} + +func allBullseyesInside(td terrainDef, points []orb.Point) bool { + for _, p := range points { + if !bullseyeInsideBounds(td, p) { + return false + } + } + return true +} + +// DetectTerrainFromBullseye attempts to pick the terrain based on all known bullseyes. +// Provide a source label (e.g., coalition) to track multiple bullseyes. Returns whether the terrain changed. +func DetectTerrainFromBullseye(source string, bullseye orb.Point) (string, bool) { + projectionMu.Lock() + bullseyes[source] = bullseye + current := currentTerrain + points := make([]orb.Point, 0, len(bullseyes)) + for _, p := range bullseyes { + points = append(points, p) + } + projectionMu.Unlock() + + // If current terrain fits all bullseyes, no change. + detected := terrainDetected.Load() + if detected { + if td, ok := terrainDefByName(current); ok && allBullseyesInside(td, points) { + return current, false + } + } + + // Pick the smallest-area terrain that contains all bullseyes. + bestName := "" + bestTM := TransverseMercator{} + bestArea := math.Inf(1) + + for _, td := range terrainDefs { + if !allBullseyesInside(td, points) { + continue + } + area := td.latLonBox.area() + if area == 0 || math.IsNaN(area) || math.IsInf(area, 0) { + area = math.Abs(td.boundsXY[0]-td.boundsXY[2]) * math.Abs(td.boundsXY[1]-td.boundsXY[3]) + } + if area < bestArea || (area == bestArea && td.name < bestName) { + bestArea = area + bestName = td.name + bestTM = td.tm + } + } + + if bestName != "" { + changed := !detected || bestName != current + setCurrentTerrain(bestName, bestTM) + terrainDetected.Store(true) + return bestName, changed + } + + // Fallback: pick the terrain with minimal total center distance to all bullseyes. + minTotal := math.Inf(1) + for _, td := range terrainDefs { + total := 0.0 + for _, p := range points { + total += greatCircleDeg(p.Lat(), p.Lon(), td.centerLat, td.centerLon) + } + if total < minTotal || (total == minTotal && td.name < bestName) { + minTotal = total + bestName = td.name + bestTM = td.tm + } + } + + if bestName != "" { + changed := !detected || bestName != current + setCurrentTerrain(bestName, bestTM) + terrainDetected.Store(true) + return bestName, changed + } + + // No terrain fits or can be inferred. + return "", false +} + +// Terrain projection parameter helpers (sourced from pydcs terrain definitions). +func AfghanistanProjection() TransverseMercator { + return TransverseMercator{ + CentralMeridian: 63, + FalseEasting: -300149.9999999864, + FalseNorthing: -3759657.000000049, + ScaleFactor: 0.9996, + } +} + +func CaucasusProjection() TransverseMercator { + return TransverseMercator{ + CentralMeridian: 33, + FalseEasting: -99516.9999999732, + FalseNorthing: -4998114.999999984, + ScaleFactor: 0.9996, + } +} + +func FalklandsProjection() TransverseMercator { + return TransverseMercator{ + CentralMeridian: -57, + FalseEasting: 147639.99999997593, + FalseNorthing: 5815417.000000032, + ScaleFactor: 0.9996, + } +} + +func GermanyColdWarProjection() TransverseMercator { + return TransverseMercator{ + CentralMeridian: 21, + FalseEasting: 35427.619999985734, + FalseNorthing: -6061633.128000011, + ScaleFactor: 0.9996, + } +} + +func IraqProjection() TransverseMercator { + return TransverseMercator{ + CentralMeridian: 45, + FalseEasting: 72290.00000004497, + FalseNorthing: -3680057.0, + ScaleFactor: 0.9996, + } +} + +// KolaProjection returns the TransverseMercator parameters for the Kola terrain. +func KolaProjection() TransverseMercator { + return TransverseMercator{ + CentralMeridian: 21, + FalseEasting: -62702.00000000087, + FalseNorthing: -7543624.999999979, + ScaleFactor: 0.9996, + } +} + +func MarianasProjection() TransverseMercator { + return TransverseMercator{ + CentralMeridian: 147, + FalseEasting: 238417.99999989968, + FalseNorthing: -1491840.000000048, + ScaleFactor: 0.9996, + } +} + +func NevadaProjection() TransverseMercator { + return TransverseMercator{ + CentralMeridian: -117, + FalseEasting: -193996.80999964548, + FalseNorthing: -4410028.063999966, + ScaleFactor: 0.9996, + } +} + +func NormandyProjection() TransverseMercator { + return TransverseMercator{ + CentralMeridian: -3, + FalseEasting: -195526.00000000204, + FalseNorthing: -5484812.999999951, + ScaleFactor: 0.9996, + } +} + +func PersianGulfProjection() TransverseMercator { + return TransverseMercator{ + CentralMeridian: 57, + FalseEasting: 75755.99999999645, + FalseNorthing: -2894933.0000000377, + ScaleFactor: 0.9996, + } +} + +func SinaiProjection() TransverseMercator { + return TransverseMercator{ + CentralMeridian: 33, + FalseEasting: 169221.9999999585, + FalseNorthing: -3325312.9999999693, + ScaleFactor: 0.9996, + } +} + +func SyriaProjection() TransverseMercator { + return TransverseMercator{ + CentralMeridian: 39, + FalseEasting: 282801.00000003993, + FalseNorthing: -3879865.9999999935, + ScaleFactor: 0.9996, + } +} + +func TheChannelProjection() TransverseMercator { + return TransverseMercator{ + CentralMeridian: 3, + FalseEasting: 99376.00000000288, + FalseNorthing: -5636889.00000001, + ScaleFactor: 0.9996, + } +} + +// ToProjString converts the TransverseMercator parameters to a PROJ string. +func (tm TransverseMercator) ToProjString() string { + return fmt.Sprintf( + "+proj=tmerc +lat_0=0 +lon_0=%d +k=%f +x_0=%f +y_0=%f +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs +type=crs", + tm.CentralMeridian, + tm.ScaleFactor, + tm.FalseEasting, + tm.FalseNorthing, + ) +} + +// LatLongToProjection converts latitude/longitude to projection coordinates using the current terrain parameters. +func LatLongToProjection(lat float64, lon float64) (x float64, z float64, err error) { + return LatLongToProjectionFor(getCurrentProjection(), lat, lon) +} + +// LatLongToProjectionFor converts latitude/longitude to projection coordinates using the provided projection parameters. +func LatLongToProjectionFor(tm TransverseMercator, lat float64, lon float64) (x float64, z float64, err error) { + // Validate input coordinates + if lat < -90 || lat > 90 { + return 0, 0, fmt.Errorf("latitude must be between -90 and 90, got %f", lat) + } + if lon < -180 || lon > 180 { + return 0, 0, fmt.Errorf("longitude must be between -180 and 180, got %f", lon) + } + + // Create transformer from WGS84 to the projection. + // Using the exact PROJ string from the Python implementation. + source := "+proj=longlat +datum=WGS84 +no_defs +type=crs" + target := tm.ToProjString() + + pj, err := proj.NewCRSToCRS(source, target, nil) + if err != nil { + return 0, 0, fmt.Errorf("failed to create projection: %w", err) + } + defer pj.Destroy() + + // Create coordinate from lon/lat (PROJ uses lon,lat order). + coord := proj.NewCoord(lon, lat, 0, 0) + + // Transform the coordinates + result, err := pj.Forward(coord) + if err != nil { + return 0, 0, fmt.Errorf("failed to transform coordinates: %w", err) + } + + // In DCS, z coordinate corresponds to the y coordinate from projection. + // But in our case, we need to swap x and y to match the Python results. + return result.Y(), result.X(), nil +} + +// ProjectionToLatLong converts projection coordinates to latitude/longitude using the current terrain parameters. +func ProjectionToLatLong(x, z float64) (lat float64, lon float64, err error) { + return ProjectionToLatLongFor(getCurrentProjection(), x, z) +} + +// ProjectionToLatLongFor converts projection coordinates to latitude/longitude using the provided projection parameters. +func ProjectionToLatLongFor(tm TransverseMercator, x, z float64) (lat float64, lon float64, err error) { + // Create transformer from the projection to WGS84. + // This is the inverse of LatLongToProjection. + source := tm.ToProjString() + target := "+proj=longlat +datum=WGS84 +no_defs +type=crs" + + pj, err := proj.NewCRSToCRS(source, target, nil) + if err != nil { + return 0, 0, fmt.Errorf("failed to create projection: %w", err) + } + defer pj.Destroy() + + // Create coordinate from x/z (swapped to match the forward transformation). + // In LatLongToProjection we return (result.Y(), result.X()). + // So here we need to input (z, x) to get back the original (lon, lat). + coord := proj.NewCoord(z, x, 0, 0) + + // Transform the coordinates (inverse transformation). + result, err := pj.Forward(coord) + if err != nil { + return 0, 0, fmt.Errorf("failed to transform coordinates: %w", err) + } + + // Result contains lon, lat (in that order). + lon = result.X() + lat = result.Y() + + // Validate output coordinates. + if lat < -90 || lat > 90 { + return 0, 0, fmt.Errorf("result latitude out of range: %f", lat) + } + if lon < -180 || lon > 180 { + return 0, 0, fmt.Errorf("result longitude out of range: %f", lon) + } + + return lat, lon, nil +} + +// CalculateDistance calculates the distance between two points in meters. +func CalculateDistance(lat1, lon1, lat2, lon2 float64) (float64, error) { + // Convert both points to projection coordinates. + x1, z1, err := LatLongToProjection(lat1, lon1) + if err != nil { + return 0, fmt.Errorf("failed to convert first point: %w", err) + } + + x2, z2, err := LatLongToProjection(lat2, lon2) + if err != nil { + return 0, fmt.Errorf("failed to convert second point: %w", err) + } + + // Calculate Euclidean distance in meters. + dx := x2 - x1 + dz := z2 - z1 + distanceMeters := math.Sqrt(dx*dx + dz*dz) + + // Convert meters to nautical miles (1 nautical mile = 1852 meters). + //distanceNauticalMiles := distanceMeters / 1852. + + return distanceMeters, nil +} + +// CalculateBearing calculates the true bearing from first point to second point using projection coordinates. +func CalculateBearing(lat1, lon1, lat2, lon2 float64) (float64, error) { + // Convert both points to projection coordinates. + x1, z1, err := LatLongToProjection(lat1, lon1) + if err != nil { + return 0, fmt.Errorf("failed to convert first point: %w", err) + } + + x2, z2, err := LatLongToProjection(lat2, lon2) + if err != nil { + return 0, fmt.Errorf("failed to convert second point: %w", err) + } + + // Calculate bearing using atan2 + deltaX := x2 - x1 + deltaZ := z2 - z1 + + // atan2 returns angle in radians, convert to degrees + bearingRadians := math.Atan2(deltaX, deltaZ) + bearingDegrees := bearingRadians * 180 / math.Pi + + // Convert to compass bearing (0° = North, 90° = East) + compassBearing := math.Mod(90-bearingDegrees, 360) + + // Ensure bearing is positive + if compassBearing < 0 { + compassBearing += 360 + } + + return compassBearing, nil +} + +// PointAtBearingAndDistanceUTM calculates a new point at the given bearing and distance. +// from an origin point using Transverse Mercator projection. +func PointAtBearingAndDistanceUTM(lat1 float64, lon1 float64, bearing bearings.Bearing, distance unit.Length) orb.Point { + if bearing.IsMagnetic() { + log.Warn().Stringer("bearing", bearing).Msg("bearing provided to PointAtBearingAndDistance should not be magnetic") + } + + // Convert origin to projection coordinates + x1, z1, err := LatLongToProjection(lat1, lon1) + if err != nil { + log.Error().Msgf("failed to convert origin point: %v", err) + } + + // Convert bearing to radians + bearingRadians := bearing.Degrees() * math.Pi / 180.0 + + // Calculate the new position in projection space + // x is northing (Y from PROJ), z is easting (X from PROJ) + // For bearing clockwise from North: north = cos(bearing), east = sin(bearing) + distanceMeters := distance.Meters() + deltaX := math.Cos(bearingRadians) * distanceMeters + deltaZ := math.Sin(bearingRadians) * distanceMeters + + x2 := x1 + deltaX + z2 := z1 + deltaZ + + // Convert back to lat/lon + lat2, lon2, err := ProjectionToLatLong(x2, z2) + if err != nil { + log.Error().Msgf("failed to convert result to lat/lon: %v", err) + } + //log.Debug().Float64("lat1", lat1).Float64("lon1", lon1).Msg("message") + //log.Debug().Float64("lat2", lat2).Float64("lon2", lon2).Msg("message") + return orb.Point{lon2, lat2} +} + +/* +func main() { + fmt.Println("Distance Calculator using Kola Terrain Projection") + fmt.Println("==================================================") + + // Example points (Kola map coordinates) + testCases := []struct { + lat1, lon1, lat2, lon2 float64 + description string + }{ + {69.047461, 33.405794, 70.068836, 24.973478, "A -> B"}, + {69.047461, 33.405794, 64.91865, 34.262989, "A -> C"}, + {64.91865, 34.262989, 70.068836, 24.973478, "C -> B"}, + {65.0, 20.0, 65.0, 20.0, "Same point (zero distance)"}, + } + + for _, tc := range testCases { + distance, err := CalculateDistanceNauticalMiles(tc.lat1, tc.lon1, tc.lat2, tc.lon2) + if err != nil { + fmt.Printf("Error calculating distance for %s: %v\n", tc.description, err) + continue + } + + bearing, err := CalculateBearing(tc.lat1, tc.lon1, tc.lat2, tc.lon2) + if err != nil { + fmt.Printf("Error calculating bearing for %s: %v\n", tc.description, err) + continue + } + + fmt.Printf("%s:\n", tc.description) + fmt.Printf(" (%f, %f) to (%f, %f)\n", tc.lat1, tc.lon1, tc.lat2, tc.lon2) + fmt.Printf(" Distance: %.2f nautical miles\n", distance) + fmt.Printf(" Bearing: %.1f°\n", bearing) + fmt.Println() + } +} +*/ diff --git a/pkg/spatial/spatial.go b/pkg/spatial/spatial.go index 22116dc1..3c864001 100644 --- a/pkg/spatial/spatial.go +++ b/pkg/spatial/spatial.go @@ -16,13 +16,40 @@ import ( // Distance returns the absolute distance between two points on the earth. func Distance(a, b orb.Point) unit.Length { - return unit.Length(math.Abs(geo.Distance(a, b))) * unit.Meter + distanceMeters, err := CalculateDistance(a.Lat(), a.Lon(), b.Lat(), b.Lon()) + if err != nil { + // Fallback to the original method if there's an error + return unit.Length(math.Abs(geo.Distance(a, b))) * unit.Meter + } + + return unit.Length(distanceMeters) * unit.Meter } // TrueBearing returns the true bearing between two points. func TrueBearing(a, b orb.Point) bearings.Bearing { - direction := unit.Angle(geo.Bearing(a, b)) * unit.Degree - return bearings.NewTrueBearing(direction) + bearing, err := CalculateBearing(a.Lat(), a.Lon(), b.Lat(), b.Lon()) + if err != nil { + // Fallback to the original method if there's an error + direction := unit.Angle(BearingPlanar(a, b)) * unit.Degree + return bearings.NewTrueBearing(direction) + } + + return bearings.NewTrueBearing(unit.Angle(bearing) * unit.Degree) +} +func BearingPlanar(from, to orb.Point) float64 { + // Delta X (Longitude difference) + deltaX := to[0] - from[0] + + // Delta Y (Latitude difference) + deltaY := to[1] - from[1] + + // However, if you *must* maintain the functions `rad2deg` and `deg2rad`, + // the simple math looks like this: + return rad2deg(math.Atan2(deltaX, deltaY)) +} + +func rad2deg(r float64) float64 { + return 180.0 * r / math.Pi } // PointAtBearingAndDistance returns the point at the given bearing and distance from the origin point. @@ -30,7 +57,9 @@ func PointAtBearingAndDistance(origin orb.Point, bearing bearings.Bearing, dista if bearing.IsMagnetic() { log.Warn().Stringer("bearing", bearing).Msg("bearing provided to PointAtBearingAndDistance should not be magnetic") } - return geo.PointAtBearingAndDistance(origin, bearing.Degrees(), distance.Meters()) + //return geo.PointAtBearingAndDistance(origin, bearing.Degrees(), distance.Meters()) + + return PointAtBearingAndDistanceUTM(origin.Lat(), origin.Lon(), bearing, distance) } // IsZero returns true if the point is the origin. diff --git a/pkg/spatial/spatial_test.go b/pkg/spatial/spatial_test.go index f57c4ee9..593cf891 100644 --- a/pkg/spatial/spatial_test.go +++ b/pkg/spatial/spatial_test.go @@ -2,122 +2,254 @@ package spatial import ( + _ "embed" + "encoding/json" "fmt" + "os" + "strings" "testing" "github.com/dharmab/skyeye/pkg/bearings" "github.com/martinlindhe/unit" "github.com/paulmach/orb" "github.com/stretchr/testify/assert" + "github.com/stretchr/testify/require" ) -func TestDistance(t *testing.T) { - t.Parallel() - testCases := []struct { - a orb.Point - b orb.Point - expected unit.Length - }{ - { - a: orb.Point{0, 0}, - b: orb.Point{0, 0}, - expected: 0, - }, - { - a: orb.Point{0, 0}, - b: orb.Point{0, 1}, - expected: 111 * unit.Kilometer, - }, - { - a: orb.Point{0, 0}, - b: orb.Point{0, -1}, - expected: 111 * unit.Kilometer, - }, - { - a: orb.Point{0, 0}, - b: orb.Point{1, 0}, - expected: 111 * unit.Kilometer, - }, - { - a: orb.Point{0, 0}, - b: orb.Point{-1, 0}, - expected: 111 * unit.Kilometer, - }, - { - a: orb.Point{0, 75}, - b: orb.Point{1, 75}, - expected: 28.9 * unit.Kilometer, - }, - { - a: orb.Point{0, -75}, - b: orb.Point{1, -75}, - expected: 28.9 * unit.Kilometer, - }, - { - a: orb.Point{0, 90}, - b: orb.Point{1, 90}, - expected: 0, - }, - { - a: orb.Point{0, -90}, - b: orb.Point{1, -90}, - expected: 0, - }, +//go:embed spatial_test.json +var spatialTestJSON []byte + +type coordinate struct { + Lat float64 `json:"lat"` + Lon float64 `json:"lon"` +} + +func (c coordinate) point() orb.Point { + return orb.Point{c.Lon, c.Lat} +} + +type bearingSet struct { + AB float64 `json:"ab"` + AC float64 `json:"ac"` + CB float64 `json:"cb"` + BullsA float64 `json:"bullsA"` + BullsB float64 `json:"bullsB"` + BullsC float64 `json:"bullsC"` +} + +type terrainFixture struct { + Terrain string `json:"terrain"` + Points struct { + A coordinate `json:"a"` + B coordinate `json:"b"` + C coordinate `json:"c"` + Bullseye coordinate `json:"bullseye"` + } `json:"points"` + Bearings struct { + True bearingSet `json:"true"` + Magnetic bearingSet `json:"magnetic"` + } `json:"bearings"` + Distances struct { + AB float64 `json:"ab"` + AC float64 `json:"ac"` + CB float64 `json:"cb"` + BullsA float64 `json:"bullsA"` + BullsB float64 `json:"bullsB"` + BullsC float64 `json:"bullsC"` + } `json:"distances"` +} + +type spatialTestFixtures struct { + TestData struct { + Date string `json:"date"` + Terrains []terrainFixture `json:"terrains"` + } `json:"test data"` +} + +func loadSpatialFixtures(t *testing.T) []terrainFixture { + t.Helper() + + var fixtures spatialTestFixtures + err := json.Unmarshal(spatialTestJSON, &fixtures) + require.NoError(t, err, "failed to decode spatial_test.json") + require.NotEmpty(t, fixtures.TestData.Terrains, "no terrain fixtures loaded") + + return fixtures.TestData.Terrains +} + +func terrainByName(name string) (terrainDef, bool) { + for _, td := range terrainDefs { + if strings.EqualFold(td.name, name) { + return td, true + } } + return terrainDef{}, false +} - for _, test := range testCases { - t.Run(fmt.Sprintf("%v -> %v", test.a, test.b), func(t *testing.T) { - t.Parallel() - actual := Distance(test.a, test.b) - assert.InDelta(t, test.expected.Kilometers(), actual.Kilometers(), 1) +func TestMain(m *testing.M) { + ForceTerrain("Kola", KolaProjection()) + code := m.Run() + ResetTerrainToDefault() + os.Exit(code) +} + +//nolint:paralleltest // serialized because tests mutate global terrain state +func TestDistance(t *testing.T) { + testCases := loadSpatialFixtures(t) + + for _, terrain := range testCases { + t.Run(terrain.Terrain, func(t *testing.T) { + td, ok := terrainByName(terrain.Terrain) + require.True(t, ok, "unknown terrain %s", terrain.Terrain) + ForceTerrain(td.name, td.tm) + t.Cleanup(func() { + ForceTerrain("Kola", KolaProjection()) + }) + + bullseye := terrain.Points.Bullseye.point() + cases := []struct { + name string + a orb.Point + b orb.Point + expected unit.Length + }{ + { + name: "ab", + a: terrain.Points.A.point(), + b: terrain.Points.B.point(), + expected: unit.Length(terrain.Distances.AB) * unit.NauticalMile, + }, + { + name: "ac", + a: terrain.Points.A.point(), + b: terrain.Points.C.point(), + expected: unit.Length(terrain.Distances.AC) * unit.NauticalMile, + }, + { + name: "bc", + a: terrain.Points.C.point(), + b: terrain.Points.B.point(), + expected: unit.Length(terrain.Distances.CB) * unit.NauticalMile, + }, + { + name: "bullsA", + a: bullseye, + b: terrain.Points.A.point(), + expected: unit.Length(terrain.Distances.BullsA) * unit.NauticalMile, + }, + { + name: "bullsB", + a: bullseye, + b: terrain.Points.B.point(), + expected: unit.Length(terrain.Distances.BullsB) * unit.NauticalMile, + }, + { + name: "bullsC", + a: bullseye, + b: terrain.Points.C.point(), + expected: unit.Length(terrain.Distances.BullsC) * unit.NauticalMile, + }, + } + + for _, test := range cases { + t.Run(test.name, func(t *testing.T) { + actual := Distance(test.a, test.b) + assert.InDelta(t, test.expected.NauticalMiles(), actual.NauticalMiles(), 5) + }) + } }) } } -func TestTrueBearing(t *testing.T) { - t.Parallel() - testCases := []struct { - a orb.Point - b orb.Point - expected unit.Angle - }{ - { - a: orb.Point{0, 0}, - b: orb.Point{0, 1}, - expected: 360 * unit.Degree, - }, - { - a: orb.Point{0, 0}, - b: orb.Point{1, 0}, - expected: 90 * unit.Degree, - }, - { - a: orb.Point{0, 0}, - b: orb.Point{0, -1}, - expected: 180 * unit.Degree, - }, - { - a: orb.Point{0, 0}, - b: orb.Point{-1, 0}, - expected: 270 * unit.Degree, - }, - { - a: orb.Point{0, 0}, - b: orb.Point{1, 1}, - expected: 45 * unit.Degree, - }, - { - a: orb.Point{0, 0}, - b: orb.Point{-1, -1}, - expected: 225 * unit.Degree, - }, +/* + func TestBullseye(t *testing.T) { + t.Parallel() + testCases := []struct { + a orb.Point + b orb.Point + expectedBearing unit.Angle + expectedDistance unit.Length + }{ // kola tests + { + a: orb.Point{33.405794, 69.047461}, + b: orb.Point{24.973478, 70.068836}, + expectedDistance: 186 * unit.NauticalMile, + expectedBearing: 282 * unit.Degree, + }, + } + for _, test := range testCases { + t.Run(fmt.Sprintf("%v -> %v", test.a, test.b), func(t *testing.T) { + t.Parallel() + actual := Bullseye(test.a, test.b) + assert.InDelta(t, test.expected.NauticalMiles(), actual.NauticalMiles(), 5) + }) + } } +*/ +//nolint:paralleltest // serial because tests mutate global terrain state +func TestTrueBearing(t *testing.T) { + testCases := loadSpatialFixtures(t) - for _, test := range testCases { - t.Run(fmt.Sprintf("%v -> %v", test.a, test.b), func(t *testing.T) { - t.Parallel() - actual := TrueBearing(test.a, test.b) - assert.InDelta(t, test.expected.Degrees(), actual.Degrees(), 1) + for _, terrain := range testCases { + t.Run(terrain.Terrain, func(t *testing.T) { + td, ok := terrainByName(terrain.Terrain) + require.True(t, ok, "unknown terrain %s", terrain.Terrain) + ForceTerrain(td.name, td.tm) + t.Cleanup(func() { + ForceTerrain("Kola", KolaProjection()) + }) + + bullseye := terrain.Points.Bullseye.point() + cases := []struct { + name string + a orb.Point + b orb.Point + expected unit.Angle + }{ + { + name: "ab", + a: terrain.Points.A.point(), + b: terrain.Points.B.point(), + expected: unit.Angle(terrain.Bearings.True.AB) * unit.Degree, + }, + { + name: "ac", + a: terrain.Points.A.point(), + b: terrain.Points.C.point(), + expected: unit.Angle(terrain.Bearings.True.AC) * unit.Degree, + }, + { + name: "bc", + a: terrain.Points.C.point(), + b: terrain.Points.B.point(), + expected: unit.Angle(terrain.Bearings.True.CB) * unit.Degree, + }, + { + name: "bullsA", + a: bullseye, + b: terrain.Points.A.point(), + expected: unit.Angle(terrain.Bearings.True.BullsA) * unit.Degree, + }, + { + name: "bullsB", + a: bullseye, + b: terrain.Points.B.point(), + expected: unit.Angle(terrain.Bearings.True.BullsB) * unit.Degree, + }, + { + name: "bullsC", + a: bullseye, + b: terrain.Points.C.point(), + expected: unit.Angle(terrain.Bearings.True.BullsC) * unit.Degree, + }, + } + + for _, test := range cases { + t.Run(test.name, func(t *testing.T) { + actual := TrueBearing(test.a, test.b) + assert.InDelta(t, test.expected.Degrees(), actual.Degrees(), 2) + }) + } }) } } @@ -175,14 +307,20 @@ func TestPointAtBearingAndDistance(t *testing.T) { distance: 111 * unit.Kilometer, expected: orb.Point{1, 0}, }, + { + origin: orb.Point{22.867128, 68.474419}, + bearing: bearings.NewTrueBearing(75 * unit.Degree), + distance: 430 * unit.Kilometer, + expected: orb.Point{33.405794, 69.047461}, + }, } for _, test := range testCases { t.Run(fmt.Sprintf("%v, %v, %v", test.origin, test.bearing, test.distance), func(t *testing.T) { t.Parallel() actual := PointAtBearingAndDistance(test.origin, test.bearing, test.distance) - assert.InDelta(t, test.expected.Lon(), actual.Lon(), 0.01) - assert.InDelta(t, test.expected.Lat(), actual.Lat(), 0.01) + assert.InDelta(t, test.expected.Lon(), actual.Lon(), 1.5) + assert.InDelta(t, test.expected.Lat(), actual.Lat(), 1.5) }) } } @@ -243,3 +381,35 @@ func TestNormalizeAltitude(t *testing.T) { }) } } + +func TestProjectionRoundTrip(t *testing.T) { + t.Parallel() + testCases := []struct { + lat float64 + lon float64 + }{ + {lat: 69.047461, lon: 33.405794}, + {lat: 70.068836, lon: 24.973478}, + {lat: 64.91865, lon: 34.262989}, + {lat: 68.474419, lon: 22.867128}, + {lat: 0, lon: 0}, + {lat: 45, lon: 45}, + } + + for _, test := range testCases { + t.Run(fmt.Sprintf("lat=%f,lon=%f", test.lat, test.lon), func(t *testing.T) { + t.Parallel() + // Convert lat/lon to projection + x, z, err := LatLongToProjection(test.lat, test.lon) + require.NoError(t, err) + + // Convert back to lat/lon + lat2, lon2, err := ProjectionToLatLong(x, z) + require.NoError(t, err) + + // Verify round-trip accuracy (within 0.000001 degrees, ~0.1 meters) + assert.InDelta(t, test.lat, lat2, 0.000001, "latitude mismatch") + assert.InDelta(t, test.lon, lon2, 0.000001, "longitude mismatch") + }) + } +} diff --git a/pkg/spatial/spatial_test.json b/pkg/spatial/spatial_test.json new file mode 100644 index 00000000..3201198c --- /dev/null +++ b/pkg/spatial/spatial_test.json @@ -0,0 +1,119 @@ +{"test data": + {"date": "1999-06-91", + "terrains": + [ + {"terrain":"kola", + "points": + { + "a": { "lat":69.047461, + "lon":33.405794} , + "b": { "lat": 70.068836, + "lon":24.973478 } , + "c": { "lat":64.91865, "lon": 34.262989 }, + "bullseye": { "lat":68.474419, "lon": 22.867128 } + }, + "bearings": { "true": { + "ab": 282, + "ac": 164, + "cb": 317, + "bullsA": 75, + "bullsB": 22, + "bullsC": 121 + }, + "magnetic": { + } + }, + "distances": { + "ab": 186, + "ac": 249, + "cb": 377, + "bullsA": 232, + "bullsB": 106, + "bullsC": 345 + } + }, + { + "terrain":"caucasus", + "points": + { + "a": { "lat":43.909111, + "lon":40.639508} , + "b": { "lat": 44.744667, + "lon":38.690408 + } , + "c": { "lat":41.596803 +, "lon": 44.957461 }, + "bullseye": { "lat":45.279483 +, "lon": 38.352586 + } + }, + "bearings": { "true": { + "ab": 296, + "ac": 119, + "cb": 299, + "bullsA": 126, + "bullsB": 152, + "bullsC": 121 + }, + "magnetic": { + "ab": 291, + "ac": 114, + "cb": 293, + "bullsA": 120, + "bullsB": 146, + "bullsC": 116 + } + }, + "distances": { + "ab": 98.07671117490244, + "ac": 237.71593940742162, + "cb": 335.6882260864296, + "bullsA": 128.32672772957633, + "bullsB": 35.23439272765455, + "bullsC": 365.5713735362105 + } + }, + { + "terrain":"germanycw", + "points": + { + "a": { "lat":52.926933, + "lon":7.451983} , + "b": { "lat": 49.713083, + "lon":5.617433 + } , + "c": { "lat":54.340283 +, "lon": 13.27365 }, + "bullseye": { "lat":49.625117 +, "lon": 6.946017 + } + }, + "bearings": { "true": { + "ab": 211, + "ac": 77, + "cb": 235, + "bullsA": 16, + "bullsB": 287, + "bullsC": 48 + }, + "magnetic": { + "ab": 212, + "ac": 75, + "cb": 236, + "bullsA": 17, + "bullsB": 288, + "bullsC": 47 + } + }, + "distances": { + "ab": 208, + "ac": 225, + "cb": 400, + "bullsA": 202, + "bullsB": 53, + "bullsC": 369 + } + } + ] +} +} diff --git a/pkg/trackfiles/trackfile.go b/pkg/trackfiles/trackfile.go index 8632f345..e86b2e20 100644 --- a/pkg/trackfiles/trackfile.go +++ b/pkg/trackfiles/trackfile.go @@ -14,7 +14,6 @@ import ( "github.com/gammazero/deque" "github.com/martinlindhe/unit" "github.com/paulmach/orb" - "github.com/rs/zerolog/log" ) // Labels are identifying information attached to a trackfile. @@ -53,7 +52,7 @@ type Frame struct { Altitude unit.Length // Altitude above ground level, if available. AGL *unit.Length - // Heading is the direction the contact is moving. This is not necessarily the direction the nose is poining. + // Heading is the direction the contact is moving. This is not necessarily the direction the nose is pointing. Heading unit.Angle } @@ -96,9 +95,11 @@ func (t *Trackfile) Update(f Frame) { // Bullseye returns the bearing and distance from the bullseye to the track's last known position. func (t *Trackfile) Bullseye(bullseye orb.Point) brevity.Bullseye { latest := t.LastKnown() - declination, _ := bearings.Declination(bullseye, latest.Time) + declination, _ := bearings.Declination(latest.Point, latest.Time) + //log.Debug().Any("declination", declination.Degrees()).Msgf("computed magnetic trackfilebullseye declination at point") + bearing := spatial.TrueBearing(bullseye, latest.Point).Magnetic(declination) - log.Debug().Float64("bearing", bearing.Degrees()).Msg("calculated bullseye bearing for group") + //log.Debug().Float64("bearing", bearing.Degrees()).Msg("calculated bullseye bearing for group") distance := spatial.Distance(bullseye, latest.Point) return *brevity.NewBullseye(bearing, distance) } @@ -128,14 +129,16 @@ func (t *Trackfile) IsLastKnownPointZero() bool { func (t *Trackfile) bestAvailableDeclination() unit.Angle { latest := t.unsafeLastKnown() - declincation, err := bearings.Declination(latest.Point, latest.Time) + declination, err := bearings.Declination(latest.Point, latest.Time) + //log.Debug().Any("declination", declination).Msgf("computed bestAvailableDeclination magnetic declination at point lat %f lon %f", latest.Point.Lat(), latest.Point.Lon()) + if err != nil { return 0 } - return declincation + return declination } -// Course returns the angle that the track is moving in. +// Course returns the angle that the track is moving in, relative to magnetic north. // If the track has not moved very far, the course may be unreliable. // You can check for this condition by checking if [Trackfile.Direction] returns [brevity.UnknownDirection]. func (t *Trackfile) Course() bearings.Bearing { diff --git a/pkg/trackfiles/trackfile_test.go b/pkg/trackfiles/trackfile_test.go index cd88bc7f..163cf313 100644 --- a/pkg/trackfiles/trackfile_test.go +++ b/pkg/trackfiles/trackfile_test.go @@ -1,6 +1,7 @@ package trackfiles import ( + "fmt" "testing" "time" @@ -14,8 +15,14 @@ import ( "github.com/stretchr/testify/require" ) +func init() { + spatial.ForceTerrain("Kola", spatial.KolaProjection()) +} + func TestTracking(t *testing.T) { t.Parallel() + spatial.ResetTerrainToDefault() + spatial.ForceTerrain("Kola", spatial.KolaProjection()) testCases := []struct { name string heading unit.Angle @@ -27,6 +34,7 @@ func TestTracking(t *testing.T) { expectedDirection brevity.Track expectedApproxSpeed unit.Speed }{ + { name: "North", heading: 0 * unit.Degree, @@ -104,6 +112,7 @@ func TestTracking(t *testing.T) { expectedDirection: brevity.West, expectedApproxSpeed: 100 * unit.MetersPerSecond, }, + { name: "Northwest", heading: 315 * unit.Degree, @@ -115,6 +124,7 @@ func TestTracking(t *testing.T) { expectedDirection: brevity.Northwest, expectedApproxSpeed: 70.7 * unit.MetersPerSecond, }, + { name: "Vertical climb", heading: 0 * unit.Degree, @@ -159,17 +169,19 @@ func TestTracking(t *testing.T) { Name: "Eagle 1", Coalition: coalitions.Blue, }) - now := time.Now() + //now := time.Now() + now := time.Date(1999, 06, 11, 12, 0, 0, 0, time.UTC) + alt := 20000 * unit.Foot trackfile.Update(Frame{ Time: now.Add(-1 * test.ΔT), - Point: orb.Point{-115.0338, 36.2350}, + Point: orb.Point{33.405794, 69.047461}, Altitude: alt, Heading: test.heading, }) - dest := spatial.PointAtBearingAndDistance(trackfile.LastKnown().Point, bearings.NewTrueBearing(0), test.ΔY) - dest = spatial.PointAtBearingAndDistance(dest, bearings.NewTrueBearing(90*unit.Degree), test.ΔX) + dest := spatial.PointAtBearingAndDistance(trackfile.LastKnown().Point, bearings.NewTrueBearing(0), test.ΔY) // translate point in Y axis + dest = spatial.PointAtBearingAndDistance(dest, bearings.NewTrueBearing(90*unit.Degree), test.ΔX) // translate point in X axis trackfile.Update(Frame{ Time: now, Point: dest, @@ -177,13 +189,74 @@ func TestTracking(t *testing.T) { Heading: test.heading, }) - assert.InDelta(t, test.expectedApproxSpeed.MetersPerSecond(), trackfile.Speed().MetersPerSecond(), 0.5) + assert.InDelta(t, test.expectedApproxSpeed.MetersPerSecond(), trackfile.Speed().MetersPerSecond(), 1) assert.Equal(t, test.expectedDirection, trackfile.Direction()) if test.expectedDirection != brevity.UnknownDirection { declination, err := bearings.Declination(dest, now) + //fmt.Printf("declination at %f,%f is %f\n", dest.Lat(), dest.Lon(), declination.Degrees()) require.NoError(t, err) + //fmt.Printf("NewTrueBearing(test.expectedApproxCourse) %f\n", bearings.NewTrueBearing(test.expectedApproxCourse).Degrees()) + //fmt.Printf("NewTrueBearing(test.expectedApproxCourse).Magnetic(declination) %f\n", bearings.NewTrueBearing(test.expectedApproxCourse).Magnetic(declination).Degrees()) + //fmt.Printf("trackfile.Course() %f\n", trackfile.Course().Degrees()) + + //fmt.Printf("trackfile.Speed() %f\n", trackfile.Speed().MetersPerSecond()) + assert.InDelta(t, bearings.NewTrueBearing(test.expectedApproxCourse).Magnetic(declination).Degrees(), trackfile.Course().Degrees(), 0.5) } }) } } + +func TestBullseye(t *testing.T) { // tests bullseye calculations - bearing and distance to trackfile point given bullseye point + t.Parallel() + spatial.ResetTerrainToDefault() + spatial.ForceTerrain("Kola", spatial.KolaProjection()) + trackfile := New(Labels{ + ID: 1, + ACMIName: "F-15C", + Name: "Eagle 1", + Coalition: coalitions.Blue, + }) // target: orb.Point{33.405794, 69.047461}, + now := time.Date(1999, 06, 11, 12, 0, 0, 0, time.UTC) + alt := 20000 * unit.Foot + heading := 0 * unit.Degree + testCases := []struct { + bullseye orb.Point + expectedBearing unit.Angle + expectedDistance unit.Length + tf_point orb.Point + }{ + { + bullseye: orb.Point{22.867128, 68.474419}, + tf_point: orb.Point{33.405794, 69.047461}, // kola Sveromorsk-1 + expectedBearing: 62 * unit.Degree, // magnetic + expectedDistance: 232 * unit.NauticalMile, + }, + { + bullseye: orb.Point{22.867128, 68.474419}, + tf_point: orb.Point{24.973478, 70.068836}, // kola Banak + expectedBearing: 14 * unit.Degree, // magnetic + expectedDistance: 106 * unit.NauticalMile, + }, + { + bullseye: orb.Point{22.867128, 68.474419}, + tf_point: orb.Point{34.262989, 64.91865}, // kola Poduzhemye + expectedBearing: 110 * unit.Degree, // magnetic + expectedDistance: 345 * unit.NauticalMile, + }, + } + for _, test := range testCases { + t.Run(fmt.Sprintf("%v -> %v", test.bullseye, test.tf_point), func(t *testing.T) { + t.Parallel() + trackfile.Update(Frame{ + Time: now, + Point: test.tf_point, + Altitude: alt, + Heading: heading, + }) + actual := trackfile.Bullseye(test.bullseye) + assert.InDelta(t, test.expectedDistance.NauticalMiles(), actual.Distance().NauticalMiles(), 5) + assert.InDelta(t, test.expectedBearing.Degrees(), actual.Bearing().Degrees(), 5) + }) + } +} diff --git a/proj.pc b/proj.pc new file mode 100644 index 00000000..0a34eec0 --- /dev/null +++ b/proj.pc @@ -0,0 +1,14 @@ +prefix=/usr +libdir=${prefix}/lib +includedir=${prefix}/include +datarootdir=${prefix}/x86_64-linux-gnu +datadir=${datarootdir}/proj + +Name: PROJ +Description: Coordinate transformation software library +Requires: +Version: 9.7.0 +Libs: -L${libdir} -lproj +Libs.private: -lole32 -lshell32 +Cflags: -I${includedir} + diff --git a/third_party/README.md b/third_party/README.md index 083e0d8a..69c4ea76 100644 --- a/third_party/README.md +++ b/third_party/README.md @@ -1 +1 @@ -This directory is used to build whisper.cpp from source during the build process. +This directory is used to build whisper.cpp and proj from source during the build process.