Refactoring of fast formulas: create haversine class insteadof "appro…

src/GeoTools-Tests/GeodesicApproximativeFormulasTest.class.st

@@ -7,6 +7,12 @@ Class {
 	#category : #'GeoTools-Tests-Geodesic'
 }
 
+{ #category : #testing }
+GeodesicApproximativeFormulasTest class >> isDeprecated [ 
+
+	^ true
+]
+
 { #category : #accessing }
 GeodesicApproximativeFormulasTest >> geodesicFormulasClass [
 

src/GeoTools-Tests/GeodesicFormulasTest.class.st

@@ -81,3 +81,29 @@ GeodesicFormulasTest >> testDistanceInMetersFromTo [
 
 	self assert: distance equals: 0
 ]
+
+{ #category : #tests }
+GeodesicFormulasTest >> testFastModeFormulas [
+
+	self assert: geodesicFormulas fastModeFormulas class equals: geodesicFormulas class defaultFastModeFormulasClass.
+
+	geodesicFormulas fastModeFormulas: GeodesicVincentyFormulas new.
+	self assert: geodesicFormulas fastModeFormulas class equals: GeodesicVincentyFormulas.
+
+	geodesicFormulas fastModeFormulas: nil.
+	self assert: geodesicFormulas fastModeFormulas class equals: geodesicFormulas class defaultFastModeFormulasClass.
+
+]
+
+{ #category : #tests }
+GeodesicFormulasTest >> testFastModeThresholdInRadians [
+
+	self assert: geodesicFormulas fastModeThresholdInRadians equals: geodesicFormulas class defaultFastModeThresholdInRadians.
+
+	geodesicFormulas fastModeThresholdInRadians: 1e-12.
+	self assert: geodesicFormulas fastModeThresholdInRadians equals: 1e-12.
+
+	geodesicFormulas fastModeThresholdInRadians: nil.
+	self assert: geodesicFormulas fastModeThresholdInRadians equals: geodesicFormulas class defaultFastModeThresholdInRadians.
+
+]

src/GeoTools-Tests/GeodesicHaversineFormulasTest.class.st

@@ -0,0 +1,11 @@
+Class {
+	#name : #GeodesicHaversineFormulasTest,
+	#superclass : #GeodesicFormulasTest,
+	#category : #'GeoTools-Tests-Geodesic'
+}
+
+{ #category : #accessing }
+GeodesicHaversineFormulasTest >> geodesicFormulasClass [
+
+	^ GeodesicHaversineFormulas 
+]

src/GeoTools-Tests/GeodesicUtilsTest.class.st

@@ -15,6 +15,14 @@ GeodesicUtilsTest >> setUp [
 	WGS84 cleanUp
 ]
 
+{ #category : #running }
+GeodesicUtilsTest >> tearDown [ 
+
+	GeodesicUtils cleanUp.
+
+	super tearDown.
+]
+
 { #category : #tests }
 GeodesicUtilsTest >> testAbsoluteCoordinatesFromDistanceInMAzimuthInRadiansEast [
 
@@ -112,7 +120,7 @@ GeodesicUtilsTest >> testAbsoluteCoordinatesFromDistanceInMAzimuthInRadiansWest
 ]
 
 { #category : #tests }
-GeodesicUtilsTest >> testAbsoluteCoordinatesWithVincentyFromDistanceInMAzimuthInRadians [
+GeodesicUtilsTest >> testAbsoluteCoordinatesFromDistanceInMetersAzimuthInRadians [
 	"Déplacement vers l’Est sur 1000 m depuis l’équateur (0°, 0°) à azimut 90°"
 
     | start azimuth distanceInM coord |
@@ -121,65 +129,65 @@ GeodesicUtilsTest >> testAbsoluteCoordinatesWithVincentyFromDistanceInMAzimuthIn
     azimuth := 90 degreesToRadians.
     distanceInM := 1000.
 
-    coord := GeodesicUtils absoluteCoordinatesWithVincentyFrom: start distanceInM: distanceInM azimuthInRadians: azimuth.
+    coord := GeodesicUtils absoluteCoordinatesFrom: start distanceInMeters: distanceInM azimuthInRadians: azimuth.
     self assert: (coord latitudeInDegrees closeTo: 0).
     self assert: (coord longitudeInDegrees closeTo: 0.009013373 precision: 0.0001)
 ]
 
 { #category : #tests }
-GeodesicUtilsTest >> testAbsoluteCoordinatesWithVincentyFromDistanceInMAzimuthInRadians2 [
+GeodesicUtilsTest >> testAbsoluteCoordinatesFromDistanceInMetersAzimuthInRadians2 [
 	"Déplacement vers le Nord sur 1000 m depuis l’équateur (0°, 0°) à azimut 0°"
     | start azimuth distanceInM coord |
 
 	start := AbsoluteCoordinates latitudeInDegrees: 0 longitudeInDegrees: 0.
    azimuth := 0 degreesToRadians.
    distanceInM := 1000.
 
-   coord := GeodesicUtils absoluteCoordinatesWithVincentyFrom: start distanceInM: distanceInM azimuthInRadians: azimuth.
+   coord := GeodesicUtils absoluteCoordinatesFrom: start distanceInMeters: distanceInM azimuthInRadians: azimuth.
    self assert: (coord latitudeInDegrees closeTo: 0.008993216 precision: 1e-3).
    self assert: (coord longitudeInDegrees closeTo: 0.0). 
 ]
 
 { #category : #tests }
-GeodesicUtilsTest >> testAbsoluteCoordinatesWithVincentyFromDistanceInMAzimuthInRadians3 [
+GeodesicUtilsTest >> testAbsoluteCoordinatesFromDistanceInMetersAzimuthInRadians3 [
     "Déplacement vers l’Est sur 1000 m depuis l’équateur (0°, 0°) à azimut 90°"
     | start azimuth distanceInM coord |
 
     start := AbsoluteCoordinates latitudeInDegrees: 0 longitudeInDegrees: 0.
     azimuth := 90 degreesToRadians.
     distanceInM := 1000.
 
-    coord := GeodesicUtils absoluteCoordinatesWithVincentyFrom: start distanceInM: distanceInM azimuthInRadians: azimuth.
+    coord := GeodesicUtils absoluteCoordinatesFrom: start distanceInMeters: distanceInM azimuthInRadians: azimuth.
     self assert: (coord latitudeInDegrees closeTo: 0.0 precision: 1e-3).
     self assert: (coord longitudeInDegrees closeTo: 0.008993 precision: 1e-3).
 
 ]
 
 { #category : #tests }
-GeodesicUtilsTest >> testAbsoluteCoordinatesWithVincentyFromDistanceInMAzimuthInRadians4 [
+GeodesicUtilsTest >> testAbsoluteCoordinatesFromDistanceInMetersAzimuthInRadians4 [
     "Déplacement vers le Sud sur 1000 m depuis l’équateur (0°, 0°) à azimut 180°"
     | start azimuth distanceInM coord |
 
     start := AbsoluteCoordinates latitudeInDegrees: 0 longitudeInDegrees: 0.
     azimuth := 180 degreesToRadians.
     distanceInM := 1000.
 
-    coord := GeodesicUtils absoluteCoordinatesWithVincentyFrom: start distanceInM: distanceInM azimuthInRadians: azimuth.
+    coord := GeodesicUtils absoluteCoordinatesFrom: start distanceInMeters: distanceInM azimuthInRadians: azimuth.
     self assert: (coord latitudeInDegrees closeTo: (-0.008993) precision: 1e-3).
     self assert: (coord longitudeInDegrees closeTo: 0.0 precision: 1e-3).
 
 ]
 
 { #category : #tests }
-GeodesicUtilsTest >> testAbsoluteCoordinatesWithVincentyFromDistanceInMAzimuthInRadians5 [
+GeodesicUtilsTest >> testAbsoluteCoordinatesFromDistanceInMetersAzimuthInRadians5 [
     "Test d'un déplacement connu depuis Paris vers l'Est sur 1 km à azimut 90°"
 
     | start azimuth distance coord expectedLat expectedLon delta |
     start := AbsoluteCoordinates frParis.
     azimuth := 90 degreesToRadians.
     distance := 1000.
 
-    coord := GeodesicUtils absoluteCoordinatesWithVincentyFrom: start distanceInM: distance azimuthInRadians: azimuth.
+    coord := GeodesicUtils absoluteCoordinatesFrom: start distanceInMeters: distance azimuthInRadians: azimuth.
 
     expectedLat := 48.86.
     expectedLon := 2.34.

src/GeoTools-Tests/GeodesicVincentyFormulasTest.class.st

@@ -13,6 +13,56 @@ GeodesicVincentyFormulasTest >> geodesicFormulasClass [
 	^ GeodesicVincentyFormulas 
 ]
 
+{ #category : #tests }
+GeodesicVincentyFormulasTest >> testDistanceAtThresholdBoundary [
+	"Test behavior exactly at the threshold boundary"
+
+	| coord1 coord2 distance thresholdInMeters |
+
+	"Calculate distance corresponding to 1e-4 radians"
+	thresholdInMeters := 1e-4 * WGS84 semiMajorAxisInMeters. "~637 meters"
+
+	coord1 := AbsoluteCoordinates latitudeInDegrees: 48.8566 longitudeInDegrees: 2.3522 altitudeInMeters: 0.
+	coord2 := AbsoluteCoordinates latitudeInDegrees: (48.8566 + (1e-4 radiansToDegrees)) longitudeInDegrees: 2.3522 altitudeInMeters: 0.
+
+	distance := geodesicFormulas distanceInMetersFrom: coord1 to: coord2.
+
+	"Should be close to threshold distance"
+	self assert: (distance between: (thresholdInMeters * 0.9) and: (thresholdInMeters * 1.1))
+]
+
+{ #category : #tests }
+GeodesicVincentyFormulasTest >> testDistanceUsesHaversineForClosePoints [
+	"Verify that Haversine is used for points closer than the threshold"
+
+	| coord1 coord2 distance |
+
+	"Two points separated by ~500 meters (less than 1e-4 radians threshold)"
+	coord1 := AbsoluteCoordinates latitudeInDegrees: 48.8566 longitudeInDegrees: 2.3522 altitudeInMeters: 0. "Paris"
+	coord2 := AbsoluteCoordinates latitudeInDegrees: 48.8610 longitudeInDegrees: 2.3522 altitudeInMeters: 0. "~500m north"
+
+	distance := geodesicFormulas distanceInMetersFrom: coord1 to: coord2.
+
+	"Should be around 490 meters"
+	self assert: (distance between: 480 and: 500)
+]
+
+{ #category : #tests }
+GeodesicVincentyFormulasTest >> testDistanceUsesVincentyForDistantPoints [
+	"Verify that Vincenty is used for points farther than the threshold"
+
+	| coord1 coord2 distance |
+
+	"Two points separated by ~10 km (more than 1e-4 radians threshold)"
+	coord1 := AbsoluteCoordinates latitudeInDegrees: 48.8566 longitudeInDegrees: 2.3522 altitudeInMeters: 0. "Paris"
+	coord2 := AbsoluteCoordinates latitudeInDegrees: 48.9566 longitudeInDegrees: 2.3522 altitudeInMeters: 0. "~11 km north"
+
+	distance := geodesicFormulas distanceInMetersFrom: coord1 to: coord2.
+
+	"Should be around 11 km"
+	self assert: (distance between: 11000 and: 11200).
+]
+
 { #category : #tests }
 GeodesicVincentyFormulasTest >> testIterationsLimit [
 

src/GeoTools/GeodesicApproximativeFormulas.class.st

@@ -4,6 +4,12 @@ Class {
 	#category : #'GeoTools-Geodesic'
 }
 
+{ #category : #testing }
+GeodesicApproximativeFormulas class >> isDeprecated [ 
+
+	^ true
+]
+
 { #category : #tools }
 GeodesicApproximativeFormulas >> absoluteCoordinatesAlongGeodesicFrom: aFromAbsoluteCoordinate to: aToAbsoluteCoordinate atFraction: aFraction [
 

src/GeoTools/GeodesicFormulas.class.st

@@ -1,9 +1,27 @@
 Class {
 	#name : #GeodesicFormulas,
 	#superclass : #Object,
+	#instVars : [
+		'fastModeFormulas',
+		'fastModeThresholdInRadians'
+	],
 	#category : #'GeoTools-Geodesic'
 }
 
+{ #category : #accessing }
+GeodesicFormulas class >> defaultFastModeFormulasClass [
+
+	^ GeodesicHaversineFormulas 
+]
+
+{ #category : #accessing }
+GeodesicFormulas class >> defaultFastModeThresholdInRadians [
+	"Threshold below which Haversine formula is used instead of other formula (i.e. Vincenty) faster but less precise.
+	1e-4 radians ≈ 640 meters at equator"
+
+	^ 1e-4
+]
+
 { #category : #tools }
 GeodesicFormulas >> absoluteCoordinatesAlongGeodesicFrom: aFromAbsoluteCoordinate to: aToAbsoluteCoordinate atFraction: aFraction [
 
@@ -28,19 +46,46 @@ GeodesicFormulas >> distanceInMetersFrom: aFromAbsoluteCoordinate to: anEndAbsol
 	self subclassResponsibility 
 ]
 
-{ #category : #haversine }
+{ #category : #'tools - fast mode' }
+GeodesicFormulas >> fastDistanceInMetersFrom: aFromAbsoluteCoordinate to: aToAbsoluteCoordinate [
+
+	^ self fastModeFormulas distanceInMetersFrom: aFromAbsoluteCoordinate to: aToAbsoluteCoordinate
+]
+
+{ #category : #'fast mode' }
+GeodesicFormulas >> fastModeFormulas [
+
+	^ fastModeFormulas ifNil: [
+		  fastModeFormulas := self class defaultFastModeFormulasClass new ]
+]
+
+{ #category : #'fast mode' }
+GeodesicFormulas >> fastModeFormulas: anObject [
+
+	fastModeFormulas := anObject
+]
+
+{ #category : #'fast mode' }
+GeodesicFormulas >> fastModeThresholdInRadians [
+
+	^ fastModeThresholdInRadians ifNil: [
+			  fastModeThresholdInRadians := self class
+				                                defaultFastModeThresholdInRadians ]
+]
+
+{ #category : #'fast mode' }
+GeodesicFormulas >> fastModeThresholdInRadians: anObject [
+
+	fastModeThresholdInRadians := anObject
+]
+
+{ #category : #'fast mode' }
 GeodesicFormulas >> haversineDistanceFrom: aFromAbsoluteCoordinate to: aToAbsoluteCoordinate [
-	"Haversine formulas to compute distance from two positions of a big circle - efficient to compute distance around the equator"
-	"sources: https://fr.wikipedia.org/wiki/Formule_de_haversine"
 
-	| radius deltaLat deltaLon a c |	
-	radius := WGS84 semiMajorAxisInMeters.
-	deltaLat := aToAbsoluteCoordinate latitudeInRadians - aFromAbsoluteCoordinate latitudeInRadians.
-	deltaLon := aToAbsoluteCoordinate longitudeInRadians - aFromAbsoluteCoordinate longitudeInRadians.
+	self 
+		deprecated: 'Use #fastDistanceInMetersFrom:to: instead'
+		transformWith: '`@receiver haversineDistanceFrom: `@arg1 to: `@arg2' 
+		-> '`@receiver fastDistanceInMetersFrom: `@arg1 to: `@arg2'.
 
-	a := ((deltaLat / 2) sin squared) 
-			+ (aFromAbsoluteCoordinate latitudeInRadians cos * aToAbsoluteCoordinate latitudeInRadians cos * ((deltaLon / 2) sin squared)).
-	c := 2 * (a sqrt arcTan: ((1 - a) sqrt)).
-
-	^ radius * c
+	^ self fastDistanceInMetersFrom: aFromAbsoluteCoordinate to: aToAbsoluteCoordinate
 ]