///|
pub fn magnetic_field_from_components(
  north_nt : Double,
  east_nt : Double,
  down_nt : Double,
) -> MagneticField {
  let horizontal = (north_nt * north_nt + east_nt * east_nt).sqrt()
  let total = (horizontal * horizontal + down_nt * down_nt).sqrt()
  {
    north_nt,
    east_nt,
    down_nt,
    total_nt: total,
    declination_deg: normalize_longitude(
      rad_to_deg(@math.atan2(east_nt, north_nt)),
    ),
    inclination_deg: rad_to_deg(@math.atan2(down_nt, horizontal)),
  }
}

///|
pub fn magnetic_field_from_direction(
  total_nt : Double,
  declination_deg : Double,
  inclination_deg : Double,
) -> MagneticField {
  let (north, east, down) = magnetic_components(
    total_nt, declination_deg, inclination_deg,
  )
  magnetic_field_from_components(north, east, down)
}

///|
pub fn magnetic_field_add(
  a : MagneticField,
  b : MagneticField,
) -> MagneticField {
  magnetic_field_from_components(
    a.north_nt + b.north_nt,
    a.east_nt + b.east_nt,
    a.down_nt + b.down_nt,
  )
}

///|
pub fn magnetic_field_scale(
  field : MagneticField,
  factor : Double,
) -> MagneticField {
  magnetic_field_from_components(
    field.north_nt * factor,
    field.east_nt * factor,
    field.down_nt * factor,
  )
}

///|
pub fn magnetic_pole_field(
  pole : MagneticPole,
  point : GeoPoint,
) -> MagneticField {
  let distance = Double::max(
    1.0,
    haversine_distance_m(
      { lat: pole.latitude_deg, lon: pole.longitude_deg, elevation_m: 0.0 },
      point,
    ),
  )
  let direction = initial_bearing_deg(point, {
    lat: pole.latitude_deg,
    lon: pole.longitude_deg,
    elevation_m: 0.0,
  })
  let inclination = magnetic_inclination_deg(point.lat)
  magnetic_field_from_direction(
    pole.strength_nt /
    (distance / MEAN_EARTH_RADIUS_M) /
    (distance / MEAN_EARTH_RADIUS_M),
    direction,
    inclination,
  )
}

///|
pub fn magnetic_field_at(
  point : GeoPoint,
  pole? : MagneticPole = {
    latitude_deg: DIPOLE_NORTH_LAT,
    longitude_deg: DIPOLE_NORTH_LON,
    strength_nt: 50000.0,
  },
) -> MagneticField {
  magnetic_pole_field(pole, point)
}

///|
pub fn magnetic_declination_grid(grid : Grid) -> Grid {
  grid.map_cells(fn(cell) {
    magnetic_declination_deg(cell.point.lat, cell.point.lon)
  })
}

///|
pub fn magnetic_inclination_grid(grid : Grid) -> Grid {
  grid.map_cells(fn(cell) { magnetic_inclination_deg(cell.point.lat) })
}

///|
pub fn magnetic_intensity_grid(
  grid : Grid,
  pole? : MagneticPole = {
    latitude_deg: DIPOLE_NORTH_LAT,
    longitude_deg: DIPOLE_NORTH_LON,
    strength_nt: 50000.0,
  },
) -> Grid {
  grid.map_cells(fn(cell) { magnetic_field_at(cell.point, pole~).total_nt })
}

///|
pub fn magnetic_profile(
  points : ArrayView[GeoPoint],
  pole? : MagneticPole = {
    latitude_deg: DIPOLE_NORTH_LAT,
    longitude_deg: DIPOLE_NORTH_LON,
    strength_nt: 50000.0,
  },
) -> FieldSeries {
  let series : Array[SeriesPoint] = []
  let distances = profile_distances(points)
  for i in 0.. Double {
  signed_angle_difference_deg(
    magnetic_declination_deg(a.lat, a.lon),
    magnetic_declination_deg(b.lat, b.lon),
  )
}

///|
pub fn magnetic_vector_error(
  expected : MagneticField,
  actual : MagneticField,
) -> Double {
  let dn = actual.north_nt - expected.north_nt
  let de = actual.east_nt - expected.east_nt
  let dd = actual.down_nt - expected.down_nt
  (dn * dn + de * de + dd * dd).sqrt()
}

///|
pub fn magnetic_components_to_samples(
  points : ArrayView[GeoPoint],
  total_nt : Double,
  declination_deg : Double,
  inclination_deg : Double,
) -> Array[FieldSample] {
  let field = magnetic_field_from_direction(
    total_nt, declination_deg, inclination_deg,
  )
  let output : Array[FieldSample] = []
  for point in points {
    output.push({ point, value: field.total_nt, name: "magnetic_total_nt" })
  }
  output
}