import math from .math_atan import atan from .math_round_float import round_float from .math_sin_cos import sin_cos # Sine, cosine and arctangent come from math.sin-cos and math.atan, not # math.sin or math.atan2: those come from the platform maths library, which may # differ in the last bit between Python, JavaScript and Rust, and a value next # to a rounding boundary can then round differently. With the shared helpers # the unrounded distance is the same double everywhere, and math.round-float # rounds it the same way everywhere. # IUGG mean radius R1 = (2a + b) / 3 of the GRS 80 ellipsoid (Moritz, # "Geodetic Reference System 1980", Journal of Geodesy 74 (2000) 128-133). EARTH_RADIUS_METRES = 6371008.8 PI = 3.141592653589793 HALF_PI = PI / 2 DEGREES = PI / 180 def _atan2_positive(y: float, x: float) -> float: # atan2 for y >= 0 and x >= 0. x is 0 only for exactly antipodal points, # where y / x would divide by zero. if x == 0: return 0.0 if y == 0 else HALF_PI return atan(y / x) def _is_number(value: object) -> bool: return isinstance(value, (int, float)) and not isinstance(value, bool) and math.isfinite(value) def _check_latitude(value: float) -> None: if not _is_number(value): raise TypeError("latitude must be a finite number of degrees") if value < -90 or value > 90: raise ValueError("latitude must be between -90 and 90 degrees") def _check_longitude(value: float) -> None: if not _is_number(value): raise TypeError("longitude must be a finite number of degrees") if value < -180 or value > 180: raise ValueError("longitude must be between -180 and 180 degrees") def distance(from_lat: float, from_lng: float, to_lat: float, to_lng: float) -> float: """Great-circle distance in metres (haversine, IUGG mean radius), to the millimetre.""" _check_latitude(from_lat) _check_longitude(from_lng) _check_latitude(to_lat) _check_longitude(to_lng) from_lat, from_lng, to_lat, to_lng = float(from_lat), float(from_lng), float(to_lat), float(to_lng) # A longitude difference of 359 degrees is a 1 degree step across the # antimeridian; sin^2 of the half-angle is the same either way. sin_half_lat = sin_cos(((to_lat - from_lat) * DEGREES) / 2).sin sin_half_lng = sin_cos(((to_lng - from_lng) * DEGREES) / 2).sin cos1 = sin_cos(from_lat * DEGREES).cos cos2 = sin_cos(to_lat * DEGREES).cos a = sin_half_lat * sin_half_lat + cos1 * cos2 * sin_half_lng * sin_half_lng if a < 0: a = 0.0 if a > 1: a = 1.0 # atan2 rather than asin(sqrt(a)): stays accurate for antipodal points too. c = 2 * _atan2_positive(math.sqrt(a), math.sqrt(1 - a)) # Half away from zero on the exact value of the double (math.round-float). return round_float(EARTH_RADIUS_METRES * c, 3)