Merge branch 'gnomonic-projection' into mainline

.github/workflows/build.yml

@@ -0,0 +1,24 @@
+# This workflow will build a golang project
+# For more information see: https://docs.github.com/en/actions/automating-builds-and-tests/building-and-testing-go
+
+name: Build
+
+on:
+  push:
+  pull_request:
+
+jobs:
+  build:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+      - name: Set up Go
+        uses: actions/setup-go@v4
+        with:
+          go-version: '1.17'
+      - name: Display Go version
+        run: go version
+      - name: Install dependencies
+        run: go get ./...
+      - name: Build
+        run: go build -v ./...

.github/workflows/tests.yml

@@ -0,0 +1,26 @@
+# This workflow will build a golang project
+# For more information see: https://docs.github.com/en/actions/automating-builds-and-tests/building-and-testing-go
+
+name: Tests
+
+on:
+  push:
+  pull_request:
+
+jobs:
+  build:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+      - name: Set up Go
+        uses: actions/setup-go@v4
+        with:
+          go-version: '1.17'
+      - name: Display Go version
+        run: go version
+      - name: Install dependencies
+        run: go get ./...
+      - name: Build
+        run: go build -v ./...
+      - name: Tests
+        run: go test -v ./...

.gitignore

@@ -13,3 +13,6 @@
 
 # Dependency directories (remove the comment below to include it)
 # vendor/
+
+# Custom ignores
+data

geodesic/direct_test.go

@@ -5,7 +5,6 @@ import (
 	"testing"
 
 	"github.com/pymaxion/geographiclib-go/geodesic/capabilities"
-
 	"github.com/stretchr/testify/assert"
 )
 

geodesic/gnomonic.go

@@ -0,0 +1,202 @@
+package geodesic
+
+import (
+	"math"
+
+	"github.com/pymaxion/geographiclib-go/geodesic/capabilities"
+)
+
+/*
+ Gnomonic is a type representing calculations for gnomonic projections centered
+ at an arbitrary position C on the ellipsoid. This projection is derived in
+ Section 8 of C. F. F. Karney, Algorithms for geodesics, J. Geodesy 87, 43–55
+ (2013); DOI: 10.1007/s00190-012-0578-z; addenda:
+ https://geographiclib.sourceforge.io/geod-addenda.html.
+
+ The gnomonic projection of a point P on the ellipsoid is defined as follows:
+ compute the geodesic line from C to P; compute the reduced length m12, geodesic
+ scale M12, and ρ = m12/M12; finally, this gives the coordinates x and y of P in
+ gnomonic projection with x = ρ sin azi1; y = ρ cos azi1, where azi1 is the
+ azimuth of the geodesic at C. The method Forward(double, double, double, double)
+ performs the forward projection and Reverse(double, double, double, double) is
+ the inverse of the projection. The methods also return the azimuth azi of the
+ geodesic at P and reciprocal scale rk in the azimuthal direction. The scale in
+ the radial direction is 1/(rk^2).
+
+ For a sphere, ρ reduces to a tan(s12/a), where s12 is the length of the geodesic
+ from C to P, and the gnomonic projection has the property that all geodesics
+ appear as straight lines. For an ellipsoid, this property holds only for
+ geodesics intersecting the centers. However geodesic segments close to the center
+ are approximately straight.
+
+ Consider a geodesic segment of length l. Let T be the point on the geodesic
+ (extended if necessary) closest to C, the center of the projection, and t, be
+ the distance CT. To lowest order, the maximum deviation (as a true distance) of
+ the corresponding gnomonic line segment (i.e., with the same end points) from
+ the geodesic is:
+
+  K(T) − K(C)) l^2 t / 32.
+
+ where K is the Gaussian curvature.
+
+ This result applies for any surface. For an ellipsoid of revolution, consider
+ all geodesics whose end points are within a distance r of C. For a given r, the
+ deviation is maximum when the latitude of C is 45°, when endpoints are a
+ distance r away, and when their azimuths from the center are ± 45° or ± 135°. To
+ lowest order in r and the flattening f, the deviation is f (r/2a)^3 r.
+
+ CAUTION: The definition of this projection for a sphere is standard. However,
+ there is no standard for how it should be extended to an ellipsoid. The choices
+ are:
+  * Declare that the projection is undefined for an ellipsoid.
+  * Project to a tangent plane from the center of the ellipsoid. This causes great
+  ellipses to appear as straight lines in the projection; i.e., it generalizes the
+  spherical great circle to a great ellipse. This was proposed by independently by
+  Bowring and Williams in 1997.
+  * Project to the conformal sphere with the constant of integration chosen so
+  that the values of the latitude match for the center point and perform a central
+  projection onto the plane tangent to the conformal sphere at the center point.
+  This causes normal sections through the center point to appear as straight lines
+  in the projection; i.e., it generalizes the spherical great circle to a normal
+  section. This was proposed by I. G. Letoval'tsev, Generalization of the gnomonic
+  projection for a spheroid and the principal geodetic problems involved in the
+  alignment of surface routes, Geodesy and Aerophotography (5), 271–274 (1963).
+  * The projection given here. This causes geodesics close to the center point to
+  appear as straight lines in the projection; i.e., it generalizes the spherical
+  great circle to a geodesic.
+
+ The algorithms are described in
+  C. F. F. Karney, Algorithms for geodesics, J. Geodesy 87, 43-55 (2013)
+  Link: https://doi.org/10.1007/s00190-012-0578-z
+  Addenda: https://geographiclib.sourceforge.io/geod-addenda.html
+*/
+type Gnomonic struct {
+	Earth *Geodesic
+}
+
+/*
+ Forward performs the forward projection from geographic to gnomonic.
+
+  lat0: latitude of center point of projection (degrees). Should be in the range [-90°, 90°].
+  lon0: longitude of center point of projection (degrees). Should be in the range [−540°, 540°).
+  lat - latitude of point (degrees). Should be in the range [-90°, 90°].
+  lon - longitude of point (degrees). Should be in the range [−540°, 540°).
+
+ The scale of the projection is 1/rk2 in the "radial" direction, azi clockwise
+ from true north, and is 1/rk in the direction perpendicular to this. If the
+ point lies "over the horizon", i.e., if rk ≤ 0, then NaNs are returned for x and
+ y (the correct values are returned for azi and rk). A call to Forward followed
+ by a call to Reverse will return the original (lat, lon) (to within roundoff)
+ provided the point is not over the horizon.
+*/
+func (g *Gnomonic) Forward(lat0, lon0, lat, lon float64) GnomonicData {
+	fwd := newGnomonicData()
+	fwd.Lat0, fwd.Lon0, fwd.Lat, fwd.Lon = lat0, lon0, lat, lon
+
+	caps := capabilities.Azimuth | capabilities.GeodesicScale | capabilities.ReducedLength
+	inv := g.Earth.InverseWithCapabilities(lat0, lon0, lat, lon, caps)
+	fwd.Azi, fwd.Rk = inv.Azi2, inv.M12
+
+	if inv.M12 > 0 {
+		rho := inv.M12Reduced / inv.M12
+		x, y := sincosd(inv.Azi1)
+		fwd.X, fwd.Y = rho*x, rho*y
+	}
+	return fwd
+}
+
+/*
+ Reverse performs the reverse projection from gnomonic to geographic.
+
+  lat0: latitude of center point of projection (degrees). Should be in the range [-90°, 90°].
+  lon0: longitude of center point of projection (degrees). Should be in the range [−540°, 540°).
+  x: easting of point (meters).
+  y: northing of point (meters).
+
+ The returned values GnomonicData.Lat and GnomonicData.Lon are in the range [-180°, 180°].
+
+ The scale of the projection is 1/rk2 in the "radial" direction, azi clockwise
+ from true north, and is 1/rk in the direction perpendicular to this. Even though
+ all inputs should return a valid lat and lon, it's possible that the procedure
+ fails to converge for very large x or y; in this case NaNs are returned for all
+ the output arguments. A call to Reverse followed by a call to Forward will
+ return the original (x, y) (to roundoff).
+*/
+func (g *Gnomonic) Reverse(lat0, lon0, x, y float64) GnomonicData {
+	rev := newGnomonicData()
+	rev.Lat0, rev.Lon0, rev.X, rev.Y = lat0, lon0, x, y
+
+	azi0 := atan2d(x, y)
+	rho := math.Hypot(x, y)
+	a := g.Earth.EquatorialRadius()
+	s := a * math.Atan(rho/a)
+
+	little := rho <= a
+	if !little {
+		rho = 1 / rho
+	}
+
+	caps := capabilities.Latitude | capabilities.Longitude | capabilities.Azimuth | capabilities.DistanceIn |
+		capabilities.ReducedLength | capabilities.GeodesicScale
+	line := g.Earth.LineWithCapabilities(lat0, lon0, azi0, caps)
+
+	var pos Data
+	numIt := 10
+	trip := 0
+	tripEpsilon := 0.01 * math.Sqrt(epsilon)
+
+	for count := numIt; count > 0; count-- {
+		pos = line.PositionWithCapabilities(s, caps)
+		if trip > 0 {
+			break
+		}
+
+		var ds float64
+		if little {
+			ds = ((pos.M12Reduced / pos.M12) - rho) * pos.M12 * pos.M12
+		} else {
+			ds = (rho - (pos.M12 / pos.M12Reduced)) * pos.M12Reduced * pos.M12Reduced
+		}
+		s -= ds
+
+		if math.Abs(ds) < tripEpsilon*a {
+			trip++
+		}
+	}
+
+	if trip == 0 {
+		return rev
+	}
+
+	rev.Lat, rev.Lon, rev.Azi, rev.Rk = pos.Lat2, pos.Lon2, pos.Azi1, pos.M12
+	return rev
+}
+
+// GnomonicData describes the results of gnomonic projection. This is used to
+// return the results for a gnomonic projection of a point (lat, lon) given a
+// center point of projection (lat0, lon0). The returned GnomonicData structs
+// always include the parameters provided to Gnomonic.Forward and
+// Gnomonic.Reverse and it always includes the fields x, y, azi. and rk.
+type GnomonicData struct {
+	Lat0 float64
+	Lon0 float64
+	Lat  float64
+	Lon  float64
+	X    float64
+	Y    float64
+	Azi  float64
+	Rk   float64
+}
+
+func newGnomonicData() GnomonicData {
+	return GnomonicData{
+		Lat0: math.NaN(),
+		Lon0: math.NaN(),
+		Lat:  math.NaN(),
+		Lon:  math.NaN(),
+		X:    math.NaN(),
+		Y:    math.NaN(),
+		Azi:  math.NaN(),
+		Rk:   math.NaN(),
+	}
+}

geodesic/gnomonic_test.go

@@ -0,0 +1,29 @@
+package geodesic
+
+import (
+	"testing"
+
+	"github.com/stretchr/testify/assert"
+)
+
+func TestGnomonicForward(t *testing.T) {
+	t.Run("forward calc", func(t *testing.T) {
+		lat0, lon0 := 48+50/60.0, 2+20/60.0 // Paris
+		lat, lon := 50.9, 1.8               // Calais
+		g := Gnomonic{Earth: WGS84}
+		r := g.Forward(lat0, lon0, lat, lon)
+		assert.InDelta(t, -37543.7, r.X, 0.05)
+		assert.InDelta(t, 230103, r.Y, 0.25)
+	})
+}
+
+func TestGnomonicReverse(t *testing.T) {
+	t.Run("reverse calc", func(t *testing.T) {
+		lat0, lon0 := 48+50/60.0, 2+20/60.0 // Paris
+		x, y := -38e3, 230e3                // Calais
+		g := Gnomonic{Earth: WGS84}
+		r := g.Reverse(lat0, lon0, x, y)
+		assert.InDelta(t, 50.899, r.Lat, 0.0005)
+		assert.InDelta(t, 1.79353, r.Lon, 0.000005)
+	})
+}

geodesic/inverse.go

@@ -99,7 +99,7 @@ func (s *inverseSolver) genInverse(lat1, lon1, lat2, lon2 float64, caps capabili
 	// |bet2|. Alternatively (cbet1 >= -sbet1), abs(sbet2) + sbet1 is a better
 	// measure. This logic is used in assigning calp2 in Lambda12. Sometimes these
 	// quantities vanish and in that case we force bet2 = +/- bet1 exactly. An
-	// example where is is necessary is the inverse problem 48.522876735459 0
+	// example where this is necessary is the inverse problem 48.522876735459 0
 	// -48.52287673545898293 179.599720456223079643 which failed with Visual Studio
 	// 10 (Release and Debug)
 	if cbet1 < -sbet1 {

geodesic/test_cases.go

@@ -6,7 +6,6 @@ import (
 	"testing"
 
 	"github.com/pymaxion/geographiclib-go/geodesic/capabilities"
-
 	"github.com/stretchr/testify/assert"
 	"github.com/stretchr/testify/require"
 )

geodtest/geod_test.go

@@ -10,7 +10,6 @@ import (
 	"testing"
 
 	"github.com/pymaxion/geographiclib-go/geodesic"
-
 	"github.com/stretchr/testify/assert"
 )
 

go.mod

@@ -2,10 +2,10 @@ module github.com/pymaxion/geographiclib-go
 
 go 1.17
 
-require github.com/stretchr/testify v1.7.0
+require github.com/stretchr/testify v1.8.4
 
 require (
-	github.com/davecgh/go-spew v1.1.0 // indirect
+	github.com/davecgh/go-spew v1.1.1 // indirect
 	github.com/pmezard/go-difflib v1.0.0 // indirect
-	gopkg.in/yaml.v3 v3.0.0-20200313102051-9f266ea9e77c // indirect
+	gopkg.in/yaml.v3 v3.0.1 // indirect
 )

go.sum

@@ -1,11 +1,17 @@
-github.com/davecgh/go-spew v1.1.0 h1:ZDRjVQ15GmhC3fiQ8ni8+OwkZQO4DARzQgrnXU1Liz8=
 github.com/davecgh/go-spew v1.1.0/go.mod h1:J7Y8YcW2NihsgmVo/mv3lAwl/skON4iLHjSsI+c5H38=
+github.com/davecgh/go-spew v1.1.1 h1:vj9j/u1bqnvCEfJOwUhtlOARqs3+rkHYY13jYWTU97c=
+github.com/davecgh/go-spew v1.1.1/go.mod h1:J7Y8YcW2NihsgmVo/mv3lAwl/skON4iLHjSsI+c5H38=
 github.com/pmezard/go-difflib v1.0.0 h1:4DBwDE0NGyQoBHbLQYPwSUPoCMWR5BEzIk/f1lZbAQM=
 github.com/pmezard/go-difflib v1.0.0/go.mod h1:iKH77koFhYxTK1pcRnkKkqfTogsbg7gZNVY4sRDYZ/4=
 github.com/stretchr/objx v0.1.0/go.mod h1:HFkY916IF+rwdDfMAkV7OtwuqBVzrE8GR6GFx+wExME=
-github.com/stretchr/testify v1.7.0 h1:nwc3DEeHmmLAfoZucVR881uASk0Mfjw8xYJ99tb5CcY=
-github.com/stretchr/testify v1.7.0/go.mod h1:6Fq8oRcR53rry900zMqJjRRixrwX3KX962/h/Wwjteg=
+github.com/stretchr/objx v0.4.0/go.mod h1:YvHI0jy2hoMjB+UWwv71VJQ9isScKT/TqJzVSSt89Yw=
+github.com/stretchr/objx v0.5.0/go.mod h1:Yh+to48EsGEfYuaHDzXPcE3xhTkx73EhmCGUpEOglKo=
+github.com/stretchr/testify v1.7.1/go.mod h1:6Fq8oRcR53rry900zMqJjRRixrwX3KX962/h/Wwjteg=
+github.com/stretchr/testify v1.8.0/go.mod h1:yNjHg4UonilssWZ8iaSj1OCr/vHnekPRkoO+kdMU+MU=
+github.com/stretchr/testify v1.8.4 h1:CcVxjf3Q8PM0mHUKJCdn+eZZtm5yQwehR5yeSVQQcUk=
+github.com/stretchr/testify v1.8.4/go.mod h1:sz/lmYIOXD/1dqDmKjjqLyZ2RngseejIcXlSw2iwfAo=
 gopkg.in/check.v1 v0.0.0-20161208181325-20d25e280405 h1:yhCVgyC4o1eVCa2tZl7eS0r+SDo693bJlVdllGtEeKM=
 gopkg.in/check.v1 v0.0.0-20161208181325-20d25e280405/go.mod h1:Co6ibVJAznAaIkqp8huTwlJQCZ016jof/cbN4VW5Yz0=
-gopkg.in/yaml.v3 v3.0.0-20200313102051-9f266ea9e77c h1:dUUwHk2QECo/6vqA44rthZ8ie2QXMNeKRTHCNY2nXvo=
 gopkg.in/yaml.v3 v3.0.0-20200313102051-9f266ea9e77c/go.mod h1:K4uyk7z7BCEPqu6E+C64Yfv1cQ7kz7rIZviUmN+EgEM=
+gopkg.in/yaml.v3 v3.0.1 h1:fxVm/GzAzEWqLHuvctI91KS9hhNmmWOoWu0XTYJS7CA=
+gopkg.in/yaml.v3 v3.0.1/go.mod h1:K4uyk7z7BCEPqu6E+C64Yfv1cQ7kz7rIZviUmN+EgEM=

readme/readme_test.go

@@ -7,7 +7,6 @@ import (
 
 	"github.com/pymaxion/geographiclib-go/geodesic"
 	"github.com/pymaxion/geographiclib-go/geodesic/capabilities"
-
 	"github.com/stretchr/testify/assert"
 )