diff --git a/src/GeoTools-Tests/GeodesicFormulasTest.class.st b/src/GeoTools-Tests/GeodesicFormulasTest.class.st index 6d3d2b2..b86773e 100644 --- a/src/GeoTools-Tests/GeodesicFormulasTest.class.st +++ b/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 } diff --git a/src/GeoTools/GeodesicApproximativeFormulas.class.st b/src/GeoTools/GeodesicApproximativeFormulas.class.st index 078e3f0..69e8b5e 100644 --- a/src/GeoTools/GeodesicApproximativeFormulas.class.st +++ b/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 ] diff --git a/src/GeoTools/GeodesicFormulas.class.st b/src/GeoTools/GeodesicFormulas.class.st index 51e5eb4..300a26f 100644 --- a/src/GeoTools/GeodesicFormulas.class.st +++ b/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 +] diff --git a/src/GeoTools/GeodesicHaversineFormulas.class.st b/src/GeoTools/GeodesicHaversineFormulas.class.st index 9a13c85..4bfb734 100644 --- a/src/GeoTools/GeodesicHaversineFormulas.class.st +++ b/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 ] diff --git a/src/GeoTools/GeodesicVincentyFormulas.class.st b/src/GeoTools/GeodesicVincentyFormulas.class.st index 8e5f0d1..61e50c0 100644 --- a/src/GeoTools/GeodesicVincentyFormulas.class.st +++ b/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 diff --git a/src/GeoTools/WGS84.class.st b/src/GeoTools/WGS84.class.st index db0fbff..67335b5 100644 --- a/src/GeoTools/WGS84.class.st +++ b/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"