turash/bugulma/backend/internal/geospatial/distance_calculator.go

package geospatial

import (
	"fmt"
	"math"
)

// DistanceCalculatorImpl implements DistanceCalculator interface
type DistanceCalculatorImpl struct {
	config *Config
}

// NewDistanceCalculator creates a new distance calculator
func NewDistanceCalculator(config *Config) DistanceCalculator {
	return &DistanceCalculatorImpl{
		config: config,
	}
}

// HaversineDistance calculates the great-circle distance between two points
// using the Haversine formula. This is accurate for most use cases and is fast.
func (dc *DistanceCalculatorImpl) HaversineDistance(p1, p2 Point) (float64, error) {
	if err := validatePoint(p1); err != nil {
		return 0, err
	}
	if err := validatePoint(p2); err != nil {
		return 0, err
	}

	// Convert degrees to radians
	lat1Rad := p1.Latitude * math.Pi / 180
	lon1Rad := p1.Longitude * math.Pi / 180
	lat2Rad := p2.Latitude * math.Pi / 180
	lon2Rad := p2.Longitude * math.Pi / 180

	// Calculate differences
	deltaLat := lat2Rad - lat1Rad
	deltaLon := lon2Rad - lon1Rad

	// Haversine formula
	a := math.Sin(deltaLat/2)*math.Sin(deltaLat/2) +
		math.Cos(lat1Rad)*math.Cos(lat2Rad)*
			math.Sin(deltaLon/2)*math.Sin(deltaLon/2)
	c := 2 * math.Atan2(math.Sqrt(a), math.Sqrt(1-a))

	// Distance in kilometers
	distance := dc.config.EarthRadiusKm * c

	return distance, nil
}

// VincentyDistance calculates the distance between two points using Vincenty's formulae.
// This is more accurate than Haversine for longer distances as it accounts for Earth's ellipsoidal shape.
// However, it's more computationally expensive.
func (dc *DistanceCalculatorImpl) VincentyDistance(p1, p2 Point) (float64, error) {
	if err := validatePoint(p1); err != nil {
		return 0, err
	}
	if err := validatePoint(p2); err != nil {
		return 0, err
	}

	if !dc.config.UseEllipsoidalModel {
		// Fall back to Haversine if ellipsoidal model is not enabled
		return dc.HaversineDistance(p1, p2)
	}

	// WGS84 ellipsoid parameters
	const (
		a = 6378137.0         // Semi-major axis (meters)
		f = 1 / 298.257223563 // Flattening
	)

	// Convert to radians
	lat1 := p1.Latitude * math.Pi / 180
	lon1 := p1.Longitude * math.Pi / 180
	lat2 := p2.Latitude * math.Pi / 180
	lon2 := p2.Longitude * math.Pi / 180

	L := lon2 - lon1
	U1 := math.Atan((1 - f) * math.Tan(lat1))
	U2 := math.Atan((1 - f) * math.Tan(lat2))
	sinU1 := math.Sin(U1)
	cosU1 := math.Cos(U1)
	sinU2 := math.Sin(U2)
	cosU2 := math.Cos(U2)

	lambda := L
	lambdaP := 2 * math.Pi
	iterLimit := 100

	var sinLambda, cosLambda, sinSigma, cosSigma, sigma, sinAlpha, cos2Alpha, cos2SigmaM, C float64

	for math.Abs(lambda-lambdaP) > 1e-12 && iterLimit > 0 {
		sinLambda = math.Sin(lambda)
		cosLambda = math.Cos(lambda)
		sinSigma = math.Sqrt((cosU2*sinLambda)*(cosU2*sinLambda) +
			(cosU1*sinU2-sinU1*cosU2*cosLambda)*(cosU1*sinU2-sinU1*cosU2*cosLambda))
		if sinSigma == 0 {
			return 0, nil // Co-incident points
		}
		cosSigma = sinU1*sinU2 + cosU1*cosU2*cosLambda
		sigma = math.Atan2(sinSigma, cosSigma)
		sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma
		cos2Alpha = 1 - sinAlpha*sinAlpha
		cos2SigmaM = cosSigma - 2*sinU1*sinU2/cos2Alpha
		if math.IsNaN(cos2SigmaM) {
			cos2SigmaM = 0 // Equatorial line
		}
		C = f / 16 * cos2Alpha * (4 + f*(4-3*cos2Alpha))
		lambdaP = lambda
		lambda = L + (1-C)*f*sinAlpha*(sigma+C*sinSigma*(cos2SigmaM+C*cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)))
		iterLimit--
	}

	if iterLimit == 0 {
		// Formula failed to converge, fall back to Haversine
		return dc.HaversineDistance(p1, p2)
	}

	u2 := cos2Alpha * (a*a - 6378137.0*6378137.0) / (6378137.0 * 6378137.0)
	A := 1 + u2/16384*(4096+u2*(-768+u2*(320-175*u2)))
	B := u2 / 1024 * (256 + u2*(-128+u2*(74-47*u2)))
	deltaSigma := B * sinSigma * (cos2SigmaM + B/4*(cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)-
		B/6*cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM*cos2SigmaM)))
	s := 6378137.0 * A * (sigma - deltaSigma)

	// Convert from meters to kilometers
	return s / 1000, nil
}

// CalculateDistance selects the appropriate distance calculation method
func (dc *DistanceCalculatorImpl) CalculateDistance(p1, p2 Point) (float64, error) {
	switch dc.config.DefaultDistanceMethod {
	case "vincenty":
		return dc.VincentyDistance(p1, p2)
	case "haversine":
		fallthrough
	default:
		return dc.HaversineDistance(p1, p2)
	}
}

// validatePoint validates a geographic point
func validatePoint(p Point) error {
	if p.Latitude < -90 || p.Latitude > 90 {
		return fmt.Errorf("%w: got %f", ErrInvalidLatitude, p.Latitude)
	}
	if p.Longitude < -180 || p.Longitude > 180 {
		return fmt.Errorf("%w: got %f", ErrInvalidLongitude, p.Longitude)
	}
	return nil
}