fix: align stop direction calculation with OBA Java semantics

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,15 +400,15 @@ 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
 	for _, v := range values {
 		diff := v - mean
 		sumSquares += diff * diff
 	}
-	return sumSquares / float64(len(values)-1) // Sample variance
+	return sumSquares / float64(len(values))
 }
 
 func median(values []float64) float64 {

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))
 	})