///|
pub fn normal_gravity_at_height_mgal(
  latitude_deg : Double,
  height_m : Double,
) -> Double {
  let surface = normal_gravity_mgal(latitude_deg)
  let radius_ratio = MEAN_EARTH_RADIUS_M / (MEAN_EARTH_RADIUS_M + height_m)
  surface * radius_ratio * radius_ratio
}

///|
pub fn normal_gravity_vector(point : GeoPoint) -> Vector3 {
  {
    x: 0.0,
    y: 0.0,
    z: normal_gravity_at_height_mgal(point.lat, point.elevation_m),
  }
}

///|
pub fn gravity_components(
  observed_mgal : Double,
  point : GeoPoint,
  density_kg_m3? : Double = 2670.0,
  terrain_mgal? : Double = 0.0,
) -> GravityComponents {
  let normal = normal_gravity_at_height_mgal(point.lat, point.elevation_m)
  let free_air = observed_mgal -
    normal +
    free_air_correction_mgal(point.elevation_m)
  let bouguer = free_air -
    bouguer_slab_correction_mgal(point.elevation_m, density_kg_m3~)
  {
    normal_mgal: normal,
    free_air_mgal: free_air,
    bouguer_mgal: bouguer,
    terrain_mgal: bouguer + terrain_mgal,
  }
}

///|
pub fn spherical_bouguer_correction_mgal(
  height_m : Double,
  density_kg_m3? : Double = 2670.0,
) -> Double {
  let flat = bouguer_slab_correction_mgal(height_m, density_kg_m3~)
  let curvature = height_m.abs() / MEAN_EARTH_RADIUS_M
  flat * (1.0 - curvature / 2.0)
}

///|
pub fn terrain_prism_correction_mgal(
  density_kg_m3 : Double,
  horizontal_distance_m : Double,
  elevation_difference_m : Double,
) -> Double {
  if horizontal_distance_m <= 0.0 {
    0.0
  } else {
    0.00004193 *
    density_kg_m3 *
    elevation_difference_m *
    elevation_difference_m /
    (horizontal_distance_m + elevation_difference_m.abs())
  }
}

///|
pub fn terrain_correction_from_samples_mgal(
  station : GeoPoint,
  terrain : ArrayView[FieldSample],
  density_kg_m3? : Double = 2670.0,
) -> Double {
  let mut total = 0.0
  for sample in terrain {
    let distance = haversine_distance_m(station, sample.point)
    total += terrain_prism_correction_mgal(
      density_kg_m3,
      distance,
      sample.point.elevation_m - station.elevation_m,
    )
  }
  total
}

///|
pub fn gravity_anomaly_series(
  observed : ArrayView[FieldSample],
  density_kg_m3? : Double = 2670.0,
) -> FieldSeries {
  let points : Array[SeriesPoint] = []
  for sample in observed {
    let components = gravity_components(
      sample.value,
      sample.point,
      density_kg_m3~,
    )
    points.push({
      coordinate: sample.point.elevation_m,
      value: components.bouguer_mgal,
      flag: if is_finite(sample.value) {
        Good
      } else {
        Missing
      },
    })
  }
  { name: "bouguer_anomaly", points }
}

///|
pub fn gravity_anomaly_grid(
  observed : Grid,
  latitude_deg : Double,
  height_m : Double,
  density_kg_m3? : Double = 2670.0,
) -> Grid {
  let normal = normal_gravity_at_height_mgal(latitude_deg, height_m)
  observed.map(fn(value) {
    bouguer_anomaly_mgal(value, latitude_deg, height_m, density_kg_m3~) +
    normal -
    normal
  })
}

///|
pub fn normal_gravity_latitude_profile(
  start_latitude_deg : Double,
  finish_latitude_deg : Double,
  count : Int,
) -> FieldSeries {
  let points : Array[SeriesPoint] = []
  let safe = Int::max(0, count)
  if safe > 0 {
    for i in 0.. Grid {
  grid.map(fn(value) { free_air_correction_mgal(value) })
}

///|
pub fn bouguer_correction_grid(
  grid : Grid,
  density_kg_m3? : Double = 2670.0,
) -> Grid {
  grid.map(fn(value) { bouguer_slab_correction_mgal(value, density_kg_m3~) })
}

///|
pub fn gravity_residual_grid(observed : Grid, reference : Grid) -> Grid? {
  observed.combine(reference, fn(a, b) { a - b })
}

///|
pub fn gravity_uncertainty_mgal(
  height_uncertainty_m : Double,
  density_uncertainty_kg_m3 : Double,
  height_m : Double,
) -> Double {
  let height_part = 0.3086 * height_uncertainty_m.abs()
  let density_part = 0.00004193 *
    density_uncertainty_kg_m3.abs() *
    height_m.abs()
  (height_part * height_part + density_part * density_part).sqrt()
}

///|
pub fn gravity_normal_gradient_mgal_per_degree(latitude_deg : Double) -> Double {
  let epsilon = 0.001
  (
    normal_gravity_mgal(latitude_deg + epsilon) -
    normal_gravity_mgal(latitude_deg - epsilon)
  ) /
  (2.0 * epsilon)
}

///|
pub fn gravity_height_gradient_mgal_per_m(
  latitude_deg : Double,
  height_m : Double,
) -> Double {
  let epsilon = 0.1
  (
    normal_gravity_at_height_mgal(latitude_deg, height_m + epsilon) -
    normal_gravity_at_height_mgal(latitude_deg, height_m - epsilon)
  ) /
  (2.0 * epsilon)
}

///|
pub fn gravity_latitude_correction_mgal(
  observed_mgal : Double,
  from_latitude_deg : Double,
  to_latitude_deg : Double,
) -> Double {
  observed_mgal +
  normal_gravity_mgal(to_latitude_deg) -
  normal_gravity_mgal(from_latitude_deg)
}

///|
pub fn gravity_apply_reference(
  observed : ArrayView[FieldSample],
  reference_mgal : Double,
) -> Array[FieldSample] {
  let output : Array[FieldSample] = []
  for sample in observed {
    output.push({ ..sample, value: sample.value - reference_mgal })
  }
  output
}