Fix GeodesicHaversineFormulas hemisphere inversion

src/GeoTools-Tests/GeodesicFormulasTest.class.st

@@ -55,9 +55,37 @@ GeodesicFormulasTest >> testAbsoluteCoordinatesFromDistanceInMetersAzimuthInRadi
 	coordinates := geodesicFormulas
 		               absoluteCoordinatesFrom: AbsoluteCoordinates zero
 		               distanceInMeters: 1000
-		               azimuthInRadians: 0.
+		               azimuthInRadians: (0 degreesToRadians).
+	self assert: coordinates class identicalTo: AbsoluteCoordinates.
 
-	self assert: coordinates class identicalTo: AbsoluteCoordinates
+	"limit cases"
+	"at the extreme south-east, when the requested point in on top (0°) need to return a new south-east point"
+	coordinates := geodesicFormulas
+		               absoluteCoordinatesFrom: (AbsoluteCoordinates latitudeInDegrees: -70 longitudeInDegrees: 180) 
+		               distanceInMeters: 500000
+		               azimuthInRadians: (0 degreesToRadians).
+	self assert: coordinates longitudeInDegrees equals: 180.
+
+	"at the extreme north-east, when the requested point in on bottom (180°) need to return a new north-east point"
+	coordinates := geodesicFormulas
+		               absoluteCoordinatesFrom: (AbsoluteCoordinates latitudeInDegrees: 70 longitudeInDegrees: 180) 
+		               distanceInMeters: 500000
+		               azimuthInRadians: (180 degreesToRadians).
+	self assert: coordinates longitudeInDegrees equals: 180.
+
+	"at the extreme north-west, when the requested point in on bottom (180°) need to return a new north-west point"
+	coordinates := geodesicFormulas
+		               absoluteCoordinatesFrom: (AbsoluteCoordinates latitudeInDegrees: 70 longitudeInDegrees: -180) 
+		               distanceInMeters: 500000
+		               azimuthInRadians: (180 degreesToRadians).
+	self assert: coordinates longitudeInDegrees equals: -180.
+
+	"at the extreme south-west, when the requested point in on bottom (180°) need to return a new south-west point"
+	coordinates := geodesicFormulas
+		               absoluteCoordinatesFrom: (AbsoluteCoordinates latitudeInDegrees: -70 longitudeInDegrees: -180) 
+		               distanceInMeters: 500000
+		               azimuthInRadians: (0 degreesToRadians).
+	self assert: coordinates longitudeInDegrees equals: -180.
 ]
 
 { #category : #tests }

src/GeoTools/GeodesicApproximativeFormulas.class.st

@@ -81,5 +81,7 @@ GeodesicApproximativeFormulas >> azimuthInRadiansFrom: aFromAbsoluteCoordinate t
 { #category : #tools }
 GeodesicApproximativeFormulas >> distanceInMetersFrom: aFromAbsoluteCoordinate to: anEndAbsoluteCoordinate [
 
-	^ self haversineDistanceFrom: aFromAbsoluteCoordinate to: anEndAbsoluteCoordinate
+	^ self
+		  fastDistanceInMetersFrom: aFromAbsoluteCoordinate
+		  to: anEndAbsoluteCoordinate
 ]

src/GeoTools/GeodesicFormulas.class.st

@@ -89,3 +89,19 @@ GeodesicFormulas >> haversineDistanceFrom: aFromAbsoluteCoordinate to: aToAbsolu
 
 	^ self fastDistanceInMetersFrom: aFromAbsoluteCoordinate to: aToAbsoluteCoordinate
 ]
+
+{ #category : #private }
+GeodesicFormulas >> normalizeLongitude: aLongitudeInRadians preservingHemisphereOf: originalLongitudeInRadians [
+	"Normalize longitude to [-π, π] while preserving antimeridian hemisphere"
+
+	| normalized |
+	normalized := (aLongitudeInRadians + (3 * Float pi) euclideanModulo:
+		               2 * Float pi) - Float pi.
+
+	"If both original and result are on/near antimeridian, preserve hemisphere"
+	((originalLongitudeInRadians abs - Float pi) abs < 1e-10 and: [
+		 (normalized abs - Float pi) abs < 1e-10 ]) ifTrue: [
+		^ originalLongitudeInRadians sign * Float pi ].
+
+	^ normalized
+]

src/GeoTools/GeodesicHaversineFormulas.class.st

@@ -39,9 +39,7 @@ GeodesicHaversineFormulas >> absoluteCoordinatesFrom: anAbsoluteCoordinates dist
 
 	lon2 := lon1 + ((anAzimuthInRadians sin * angularDistance sin * lat1 cos) 
 		arcTan: ((angularDistance cos) - (lat1 sin * lat2 sin))).
-
-	"Normalize longitude to [-π, π]"
-	lon2 := (((lon2 + (3 * Float pi)) euclideanModulo: (2 * Float pi)) - Float pi).
+	lon2 := self normalizeLongitude: lon2 preservingHemisphereOf: lon1.
 
 	^ AbsoluteCoordinates latitudeInRadians: lat2 longitudeInRadians: lon2
 ]

src/GeoTools/GeodesicVincentyFormulas.class.st

@@ -102,7 +102,9 @@ GeodesicVincentyFormulas >> absoluteCoordinatesFrom: anAbsoluteCoordinates dista
     lambda := (sinSigma * sinAlpha1) arcTan2: ((cosU1 * cosSigma) - (sinU1 * sinSigma * cosAlpha1)).
     C := ((f / 16) * cosSqAlpha * (4 + (f * (4 - (3 * cosSqAlpha))))).
     L := lambda - ((1 - C) * f * sinAlpha1 * (sigma + (C * sinSigma * (cos2SigmaM + (C * cosSigma * (-1 + (2 * ((cos2SigmaM) squared)))))))).
-    lon2 := (((lon1 + L + (3 * Float pi)) euclideanModulo: (2 * Float pi)) - Float pi).
+
+    lon2 := lon1 + L.
+	 lon2 := self normalizeLongitude: lon2 preservingHemisphereOf: lon1.
 
     ^ AbsoluteCoordinates latitudeInRadians: lat2 longitudeInRadians: lon2
 

src/GeoTools/WGS84.class.st

@@ -25,6 +25,12 @@ WGS84 class >> cleanUp [
 	poleCurveRadius := nil
 ]
 
+{ #category : #properties }
+WGS84 class >> eccentricity [
+
+	^ self eccentricitySquared sqrt
+]
+
 { #category : #properties }
 WGS84 class >> eccentricitySquared [
 	"~0.0066943799901413165"