From 3db3ecd23da48ad81275770cc481601566bda494 Mon Sep 17 00:00:00 2001
From: Adel Mohamed <adel.erorr159@gmail.com>
Date: Wed, 25 Feb 2026 14:13:04 +0200
Subject: [PATCH] fix: align stop direction calculation with OBA Java semantics

---
 .../gtfs/advanced_direction_calculator.go     | 42 ++++++++-----------
 .../advanced_direction_calculator_test.go     |  2 +-
 2 files changed, 19 insertions(+), 25 deletions(-)

diff --git a/internal/gtfs/advanced_direction_calculator.go b/internal/gtfs/advanced_direction_calculator.go
index 5ebbfe14..b238efb1 100644
--- a/internal/gtfs/advanced_direction_calculator.go
+++ b/internal/gtfs/advanced_direction_calculator.go
@@ -221,27 +221,23 @@ func (adc *AdvancedDirectionCalculator) computeFromShapes(ctx context.Context, s
 	// Calculate mean orientation vector
 	var xs, ys []float64
 	for _, orientation := range orientations {
-		x := math.Cos(orientation)
-		y := math.Sin(orientation)
-		xs = append(xs, x)
-		ys = append(ys, y)
+		xs = append(xs, math.Cos(orientation))
+		ys = append(ys, math.Sin(orientation))
 	}
 
 	xMu := mean(xs)
 	yMu := mean(ys)
 
-	// Check for opposite directions (mean vector is zero)
-	if xMu == 0.0 && yMu == 0.0 {
-		return "" // Ambiguous direction
+	// Check for ambiguous/opposite directions (mean vector near zero)
+	if math.Abs(xMu) < 1e-6 && math.Abs(yMu) < 1e-6 {
+		return ""
 	}
 
-	// Check variance threshold
+	// Directly compare Variance to varianceThreshold
 	xVariance := variance(xs, xMu)
 	yVariance := variance(ys, yMu)
-	xStdDev := math.Sqrt(xVariance)
-	yStdDev := math.Sqrt(yVariance)
 
-	if xStdDev > adc.varianceThreshold || yStdDev > adc.varianceThreshold {
+	if xVariance > adc.varianceThreshold || yVariance > adc.varianceThreshold {
 		return "" // Too much variance
 	}
 
@@ -335,21 +331,19 @@ func (adc *AdvancedDirectionCalculator) calculateOrientationAtStop(ctx context.C
 		indexFrom = 0
 	}
 	indexTo := closestIdx + shapePointWindow
-	if indexTo > len(shapePoints) {
-		indexTo = len(shapePoints)
+	if indexTo >= len(shapePoints) {
+		indexTo = len(shapePoints) - 1
 	}
 
-	// Calculate orientation from the window
-	// Use the bearing from the first point to the last point in the window
-	if indexTo > indexFrom+1 {
+	// Calculate orientation from the window using flat-earth approximation
+	if indexTo > indexFrom {
 		fromPoint := shapePoints[indexFrom]
-		toPoint := shapePoints[indexTo-1]
+		toPoint := shapePoints[indexTo]
+
+		dx := (toPoint.Lon - fromPoint.Lon) * math.Cos(fromPoint.Lat*math.Pi/180.0)
+		dy := toPoint.Lat - fromPoint.Lat
 
-		bearing := utils.BearingBetweenPoints(fromPoint.Lat, fromPoint.Lon, toPoint.Lat, toPoint.Lon)
-		// Convert bearing (0-360°, 0=North) to mathematical angle (radians, 0=East, counterclockwise)
-		// Bearing: 0°=N, 90°=E, 180°=S, 270°=W
-		// Math angle: 0=E, π/2=N, π=W, -π/2=S
-		orientation := (90.0 - bearing) * math.Pi / 180.0
+		orientation := math.Atan2(dy, dx)
 		return orientation, nil
 	}
 
@@ -406,7 +400,7 @@ func mean(values []float64) float64 {
 }
 
 func variance(values []float64, mean float64) float64 {
-	if len(values) <= 1 {
+	if len(values) == 0 {
 		return 0
 	}
 	sumSquares := 0.0
@@ -414,7 +408,7 @@ func variance(values []float64, mean float64) float64 {
 		diff := v - mean
 		sumSquares += diff * diff
 	}
-	return sumSquares / float64(len(values)-1) // Sample variance
+	return sumSquares / float64(len(values))
 }
 
 func median(values []float64) float64 {
diff --git a/internal/gtfs/advanced_direction_calculator_test.go b/internal/gtfs/advanced_direction_calculator_test.go
index 6c2ef17e..e2ac3c64 100644
--- a/internal/gtfs/advanced_direction_calculator_test.go
+++ b/internal/gtfs/advanced_direction_calculator_test.go
@@ -138,7 +138,7 @@ func TestStatisticalFunctions(t *testing.T) {
 		values := []float64{1, 2, 3, 4, 5}
 		m := mean(values)
 		v := variance(values, m)
-		assert.InDelta(t, 2.5, v, 0.001) // Sample variance of 1,2,3,4,5 is 2.5
+		assert.InDelta(t, 2.0, v, 0.001)
 
 		assert.Equal(t, 0.0, variance([]float64{5}, 5.0))
 	})