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 }