///|
pub fn geo_to_ecef(point : GeoPoint) -> Vector3 {
  let latitude = deg_to_rad(point.lat)
  let longitude = deg_to_rad(point.lon)
  let e2 = WGS84_F * (2.0 - WGS84_F)
  let radius = WGS84_A /
    (1.0 - e2 * @math.sin(latitude) * @math.sin(latitude)).sqrt()
  {
    x: (radius + point.elevation_m) * @math.cos(latitude) * @math.cos(longitude),
    y: (radius + point.elevation_m) * @math.cos(latitude) * @math.sin(longitude),
    z: (radius * (1.0 - e2) + point.elevation_m) * @math.sin(latitude),
  }
}

///|
pub fn ecef_to_geo(vector : Vector3) -> GeoPoint {
  let longitude = @math.atan2(vector.y, vector.x)
  let horizontal = (vector.x * vector.x + vector.y * vector.y).sqrt()
  let mut latitude = @math.atan2(vector.z, horizontal * (1.0 - WGS84_F))
  for _ in 0..<8 {
    let sin_latitude = @math.sin(latitude)
    let radius = WGS84_A /
      (1.0 - WGS84_F * (2.0 - WGS84_F) * sin_latitude * sin_latitude).sqrt()
    latitude = @math.atan2(
      vector.z + WGS84_F * (2.0 - WGS84_F) * radius * sin_latitude,
      horizontal,
    )
  }
  let sin_latitude = @math.sin(latitude)
  let radius = WGS84_A /
    (1.0 - WGS84_F * (2.0 - WGS84_F) * sin_latitude * sin_latitude).sqrt()
  let elevation = horizontal / @math.cos(latitude) - radius
  {
    lat: rad_to_deg(latitude),
    lon: normalize_longitude(rad_to_deg(longitude)),
    elevation_m: elevation,
  }
}

///|
pub fn ecef_difference(a : GeoPoint, b : GeoPoint) -> Vector3 {
  vector_subtract(geo_to_ecef(b), geo_to_ecef(a))
}

///|
pub fn ecef_distance_m(a : GeoPoint, b : GeoPoint) -> Double {
  vector_magnitude(ecef_difference(a, b))
}

///|
pub fn enu_basis(origin : GeoPoint) -> (Vector3, Vector3, Vector3) {
  let lat = deg_to_rad(origin.lat)
  let lon = deg_to_rad(origin.lon)
  let east = { x: -@math.sin(lon), y: @math.cos(lon), z: 0.0 }
  let north = {
    x: -@math.sin(lat) * @math.cos(lon),
    y: -@math.sin(lat) * @math.sin(lon),
    z: @math.cos(lat),
  }
  let up = {
    x: @math.cos(lat) * @math.cos(lon),
    y: @math.cos(lat) * @math.sin(lon),
    z: @math.sin(lat),
  }
  (east, north, up)
}

///|
pub fn ecef_to_enu(origin : GeoPoint, point : GeoPoint) -> Vector3 {
  let delta = ecef_difference(origin, point)
  let (east, north, up) = enu_basis(origin)
  { x: dot(delta, east), y: dot(delta, north), z: dot(delta, up) }
}

///|
pub fn enu_to_ecef(origin : GeoPoint, delta : Vector3) -> Vector3 {
  let origin_ecef = geo_to_ecef(origin)
  let (east, north, up) = enu_basis(origin)
  vector_add(
    origin_ecef,
    vector_add(
      vector_scale(east, delta.x),
      vector_add(vector_scale(north, delta.y), vector_scale(up, delta.z)),
    ),
  )
}

///|
pub fn enu_to_geo(origin : GeoPoint, delta : Vector3) -> GeoPoint {
  ecef_to_geo(enu_to_ecef(origin, delta))
}

///|
pub fn geodetic_height(point : GeoPoint) -> Double {
  vector_magnitude(geo_to_ecef(point)) - WGS84_A
}

///|
pub fn geocentric_latitude(point : GeoPoint) -> Double {
  let vector = geo_to_ecef(point)
  rad_to_deg(
    @math.atan2(vector.z, (vector.x * vector.x + vector.y * vector.y).sqrt()),
  )
}

///|
pub fn geocentric_radius_m(point : GeoPoint) -> Double {
  vector_magnitude(geo_to_ecef(point))
}

///|
pub fn local_north_east(
  point : GeoPoint,
  bearing_deg : Double,
  distance_m : Double,
) -> GeoPoint {
  let bearing = deg_to_rad(bearing_deg)
  let north = distance_m * @math.cos(bearing)
  let east = distance_m * @math.sin(bearing)
  enu_to_geo(point, { x: east, y: north, z: 0.0 })
}

///|
pub fn split_antimeridian(
  start : GeoPoint,
  finish : GeoPoint,
) -> Array[GeoPoint] {
  let delta = finish.lon - start.lon
  if delta.abs() <= 180.0 {
    [start, finish]
  } else {
    let adjusted_finish = {
      ..finish,
      lon: if delta > 0.0 {
        finish.lon - 360.0
      } else {
        finish.lon + 360.0
      },
    }
    [start, adjusted_finish]
  }
}

///|
pub fn unwrap_point_sequence(points : ArrayView[GeoPoint]) -> Array[GeoPoint] {
  unwrap_longitudes(points)
}

///|
pub fn angular_distance_deg(a : GeoPoint, b : GeoPoint) -> Double {
  rad_to_deg(haversine_distance_m(a, b) / MEAN_EARTH_RADIUS_M)
}

///|
pub fn angular_distance_rad(a : GeoPoint, b : GeoPoint) -> Double {
  haversine_distance_m(a, b) / MEAN_EARTH_RADIUS_M
}

///|
pub fn central_angle_from_chord(chord_m : Double) -> Double {
  2.0 *
  @math.asin((chord_m / (2.0 * MEAN_EARTH_RADIUS_M)).clamp(min=-1.0, max=1.0))
}

///|
pub fn chord_distance_m(angle_rad : Double) -> Double {
  2.0 * MEAN_EARTH_RADIUS_M * @math.sin(angle_rad / 2.0)
}

///|
pub fn geodesic_area_triangle(
  a : GeoPoint,
  b : GeoPoint,
  c : GeoPoint,
) -> Double {
  let side_a = angular_distance_rad(b, c)
  let side_b = angular_distance_rad(a, c)
  let side_c = angular_distance_rad(a, b)
  let semiperimeter = (side_a + side_b + side_c) / 2.0
  let excess = 4.0 *
    @math.atan(
      (@math.tan(semiperimeter / 2.0) *
      @math.tan((semiperimeter - side_a) / 2.0) *
      @math.tan((semiperimeter - side_b) / 2.0) *
      @math.tan((semiperimeter - side_c) / 2.0)).sqrt(),
    )
  excess * MEAN_EARTH_RADIUS_M * MEAN_EARTH_RADIUS_M
}

///|
pub fn ecef_midpoint(a : GeoPoint, b : GeoPoint) -> GeoPoint {
  ecef_to_geo(vector_scale(vector_add(geo_to_ecef(a), geo_to_ecef(b)), 0.5))
}

///|
pub fn interpolate_geodesic(
  a : GeoPoint,
  b : GeoPoint,
  fraction : Double,
) -> GeoPoint {
  let f = fraction.clamp(min=0.0, max=1.0)
  let distance = haversine_distance_m(a, b)
  destination_point(a, initial_bearing_deg(a, b), distance * f)
}

///|
pub fn geodesic_polyline(
  start : GeoPoint,
  finish : GeoPoint,
  spacing_m : Double,
) -> Array[GeoPoint] {
  let distance = haversine_distance_m(start, finish)
  let count = if spacing_m <= 0.0 {
    2
  } else {
    Int::max(2, (distance / spacing_m).ceil().to_int() + 1)
  }
  profile_points(start, finish, count)
}

///|
pub fn geodesic_polyline_length(points : ArrayView[GeoPoint]) -> Double {
  polygon_perimeter_m(points, false)
}

///|
pub fn local_scale_factors(latitude_deg : Double) -> Vector3 {
  {
    x: local_degree_longitude_m(latitude_deg),
    y: local_degree_latitude_m(latitude_deg),
    z: 1.0,
  }
}

///|
pub fn local_distance_components(a : GeoPoint, b : GeoPoint) -> Vector3 {
  local_enu_delta_m(a, b)
}

///|
pub fn local_azimuth_deg(a : GeoPoint, b : GeoPoint) -> Double {
  initial_bearing_deg(a, b)
}

///|
pub fn local_slope_percent(a : GeoPoint, b : GeoPoint) -> Double {
  let horizontal = haversine_distance_m(a, b)
  if horizontal == 0.0 {
    0.0
  } else {
    100.0 * (b.elevation_m - a.elevation_m) / horizontal
  }
}

///|
pub fn local_slope_angle_deg(a : GeoPoint, b : GeoPoint) -> Double {
  rad_to_deg(@math.atan(local_slope_percent(a, b) / 100.0))
}