nx_geo_haversine.nx source

1// nx_geo_haversine.nx -- LIB: GEO-003b GREAT-CIRCLE distance (haversine), the PRECISION REFINEMENT of 2// GEO-003 equirectangular. Hardware-up + DRY: REUSES the sovereign CORDIC sin/cos (fx.nx, Q16.16), 3// the vectoring-mode CORDIC atan2 (geo_atan2 from nx_geo_query), and nx_isqrt -- no trig reinvention. 4// 5// WHY THIS RUNG (measured, not asserted): equirectangular assumes a locally-flat earth, so its error 6// grows with span -- ~0.5% regionally and several percent at continental / antipodal scale. Haversine 7// is exact great-circle on the reference sphere. The gate PROVES the exceed: for NYC->Sydney 8// (~16000 km) haversine lands within ~2% of the true distance while equirectangular is off by far 9// more -- so haversine is measurably CLOSER. (For short spans GEO-003 stays the right tool: its 10// integer deltas avoid the Q16.16 underflow of a tiny half-angle's sine.) 11// 12// All Q16.16 fixed point -> deterministic. license_tier: ORIGINAL 13import "nx_geo_query.nx" // geo_atan2 (vectoring CORDIC) -- reused, no atan reinvention 14import "nx_geo_distance.nx" // geo_distance_m -- the equirectangular baseline we refine 15import "fx.nx" 16import "nx_isqrt.nx" 17import "nx_syscalls.nx" 18const HAV_MAGIC_360000000: i64 = 360000000 19const HAV_MAGIC_180000000: i64 = 180000000 20 21const HAV_R_M: i64 = 6371000 // mean Earth radius in meters (the great-circle reference sphere) 22 23// microdegrees -> Q16.16 radians (360e6 microdeg == 2*pi == FX_TWO_PI). Multiply before divide. 24func hav_micro_to_q16(micro: i64) -> i64 { 25 return micro * FX_TWO_PI / HAV_MAGIC_360000000 26} 27 28// Great-circle distance in METERS via the haversine formula: 29// a = sin^2(dlat/2) + cos(lat1)*cos(lat2)*sin^2(dlon/2); c = 2*atan2(sqrt(a), sqrt(1-a)); d = R*c. 30// Everything in Q16.16, sin/cos/atan2 from the sovereign CORDIC. dlon is reduced to [-180,180] deg 31// first (sin^2 is wrap-invariant) so each half-angle stays within CORDIC's sweet spot. 32func geo_haversine_m(lat1: i64, lon1: i64, lat2: i64, lon2: i64) -> i64 { 33 let dlat: i64 = lat2 - lat1 34 var dlon: i64 = lon2 - lon1 35 while dlon > HAV_MAGIC_180000000 { dlon = dlon - HAV_MAGIC_360000000 } 36 while dlon < (0 - HAV_MAGIC_180000000) { dlon = dlon + HAV_MAGIC_360000000 } 37 38 let s_dlat: i64 = fx_sin(hav_micro_to_q16(dlat) / 2) // sin(dlat/2), Q16.16 39 let s_dlon: i64 = fx_sin(hav_micro_to_q16(dlon) / 2) // sin(dlon/2), Q16.16 40 let sin2_dlat: i64 = s_dlat * s_dlat / FX_ONE // sin^2(dlat/2) 41 let sin2_dlon: i64 = s_dlon * s_dlon / FX_ONE // sin^2(dlon/2) 42 43 let cos1: i64 = fx_cos(hav_micro_to_q16(lat1)) 44 let cos2: i64 = fx_cos(hav_micro_to_q16(lat2)) 45 let coscos: i64 = cos1 * cos2 / FX_ONE // cos(lat1)*cos(lat2) 46 47 var a: i64 = sin2_dlat + coscos * sin2_dlon / FX_ONE // Q16.16, clamped to [0,1] 48 if a < 0 { a = 0 } 49 if a > FX_ONE { a = FX_ONE } 50 if a == 0 { return 0 } // coincident (or sub-resolution) points: exactly 0, no atan2(0,x) residual 51 52 // sqrt of a Q16.16 value v, result in Q16.16, is isqrt(v * FX_ONE). 53 let sqrt_a: i64 = nx_isqrt(a * FX_ONE) 54 let sqrt_1ma: i64 = nx_isqrt((FX_ONE - a) * FX_ONE) 55 let c_q16: i64 = 2 * geo_atan2(sqrt_a, sqrt_1ma) // central angle in Q16.16 radians 56 return HAV_R_M * c_q16 / FX_ONE // -> meters 57}