use super::funejson::Value; use super::math_atan::atan; use super::math_sin_cos::sin_cos; // The latitude/longitude box that contains every point within a distance of a // centre, by the method of J. P. Matuschek, "Finding Points Within a Distance // of a Latitude/Longitude Using Bounding Coordinates" // (http://janmatuschek.de/LatitudeLongitudeBoundingCoordinates). // // Sine, cosine and arctangent come from math.sin-cos and math.atan instead of // f64::sin and friends, because the platform maths library may differ in the // last bit between Rust, JavaScript and Python. The shared helpers use only // + - * /, and the rest here only + - * / sqrt floor ceil, all correctly // rounded, in the same order in all three. /// IUGG mean radius R1 of the GRS 80 ellipsoid (Moritz, Journal of Geodesy 74 (2000)). const EARTH_RADIUS_METRES: f64 = 6371008.8; /// Ten millionths of a degree, about a centimetre. const SCALE: f64 = 10000000.0; /// Within a millionth of a grid step of a 7-decimal value counts as that value, /// so a centre of 51.5074 stays 51.5074 rather than becoming 51.5073999. const SNAP: f64 = 0.000001; const PI: f64 = 3.141592653589793; const HALF_PI: f64 = PI / 2.0; const DEGREES: f64 = PI / 180.0; /// asin for 0 <= x, on math.atan: asin x = atan(x / sqrt(1 - x^2)), with /// (1 - x)(1 + x) to keep precision near 1. There is no shared asin; x >= 1 /// (reachable only by rounding, just short of a pole) is a quarter turn. fn asin_positive(x: f64) -> f64 { if x >= 1.0 { return HALF_PI; } atan(x / ((1.0 - x) * (1.0 + x)).sqrt()) } /// Outward rounding: a lower bound goes down and an upper bound goes up, so /// rounding can only grow the box, never cut the circle. fn round_down(x: f64) -> f64 { (x * SCALE + SNAP).floor() / SCALE + 0.0 } fn round_up(x: f64) -> f64 { (x * SCALE - SNAP).ceil() / SCALE + 0.0 } fn check_number(value: f64, what: &str) { if !value.is_finite() { panic!("{} must be a finite number", what); } } /// The box containing every point within `distance_metres` of (`lat`, `lng`). /// /// # Panics /// Panics if a coordinate is out of range or the distance is negative or not finite. pub fn bounding_box(lat: f64, lng: f64, distance_metres: f64) -> BoundingBox { check_number(lat, "latitude"); check_number(lng, "longitude"); check_number(distance_metres, "distance"); if lat < -90.0 || lat > 90.0 { panic!("latitude must be between -90 and 90 degrees"); } if lng < -180.0 || lng > 180.0 { panic!("longitude must be between -180 and 180 degrees"); } if distance_metres < 0.0 { panic!("distance must be 0 or more metres"); } let r = distance_metres / EARTH_RADIUS_METRES; let r_degrees = r / DEGREES; let mut min_lat = lat - r_degrees; let mut max_lat = lat + r_degrees; let mut min_lng; let mut max_lng; if min_lat > -90.0 && max_lat < 90.0 { let sin_r = sin_cos(r).sin; let cos_lat = sin_cos(lat * DEGREES).cos; let d_lng = asin_positive(sin_r / cos_lat) / DEGREES; min_lng = lng - d_lng; max_lng = lng + d_lng; // Past the antimeridian the bound comes round the other side, leaving // min_lng > max_lng: the box is the two strips either side of 180. if min_lng < -180.0 { min_lng += 360.0; } if max_lng > 180.0 { max_lng -= 360.0; } } else { // The circle reaches a pole, so it covers every longitude. if min_lat < -90.0 { min_lat = -90.0; } if max_lat > 90.0 { max_lat = 90.0; } min_lng = -180.0; max_lng = 180.0; } BoundingBox { min_lat: round_down(min_lat), min_lng: round_down(min_lng), max_lat: round_up(max_lat), max_lng: round_up(max_lng), } } pub fn bounding_box_to_value(b: &BoundingBox) -> Value { Value::obj(vec![ ("minLat", Value::Float(b.min_lat)), ("minLng", Value::Float(b.min_lng)), ("maxLat", Value::Float(b.max_lat)), ("maxLng", Value::Float(b.max_lng)), ]) } pub fn fune_vector(args: &[Value]) -> Value { // Refuse what the typed signature cannot hold, with the wording TypeScript // and Python use, rather than let the conversion below quietly change it. if !matches!(args[0], Value::Int(_) | Value::Float(_)) { panic!("latitude must be a finite number"); } if !matches!(args[1], Value::Int(_) | Value::Float(_)) { panic!("longitude must be a finite number"); } if !matches!(args[2], Value::Int(_) | Value::Float(_)) { panic!("distance must be a finite number"); } bounding_box_to_value(&bounding_box( args[0].as_f64(), args[1].as_f64(), args[2].as_f64(), )) }