///|
pub fn coefficient_table(
  reference_radius_m~ : Double,
  scale~ : Double,
  terms~ : Array[HarmonicTerm],
  source? : String = "custom",
) -> CoefficientTable {
  { reference_radius_m, scale, terms, source }
}

///|
pub fn coefficient_table_to_model(
  table : CoefficientTable,
) -> SphericalHarmonicModel {
  {
    reference_radius_m: table.reference_radius_m,
    scale: table.scale,
    terms: table.terms,
  }
}

///|
pub fn coefficient_table_is_valid(table : CoefficientTable) -> Bool {
  table.reference_radius_m > 0.0 &&
  is_finite(table.reference_radius_m) &&
  is_finite(table.scale) &&
  table.terms.all(valid_term)
}

///|
pub fn coefficient_terms_at_degree(
  table : CoefficientTable,
  degree : Int,
) -> Array[HarmonicTerm] {
  let output : Array[HarmonicTerm] = []
  for term in table.terms {
    if term.degree == degree {
      output.push(term)
    }
  }
  output
}

///|
pub fn coefficient_terms_at_order(
  table : CoefficientTable,
  order : Int,
) -> Array[HarmonicTerm] {
  let output : Array[HarmonicTerm] = []
  for term in table.terms {
    if term.order == order {
      output.push(term)
    }
  }
  output
}

///|
pub fn coefficient_max_degree(table : CoefficientTable) -> Int {
  let mut maximum = 0
  for term in table.terms {
    maximum = Int::max(maximum, term.degree)
  }
  maximum
}

///|
pub fn coefficient_degree_summaries(
  table : CoefficientTable,
) -> Array[DegreeSummary] {
  let output : Array[DegreeSummary] = []
  for degree in 0..<=coefficient_max_degree(table) {
    let mut terms = 0
    let mut cosine_energy = 0.0
    let mut sine_energy = 0.0
    for term in table.terms {
      if term.degree == degree {
        terms += 1
        cosine_energy += term.cosine * term.cosine
        sine_energy += term.sine * term.sine
      }
    }
    if terms > 0 {
      output.push({ degree, terms, cosine_energy, sine_energy })
    }
  }
  output
}

///|
pub fn coefficient_energy(table : CoefficientTable) -> Double {
  let mut total = 0.0
  for term in table.terms {
    total += term.cosine * term.cosine + term.sine * term.sine
  }
  total
}

///|
pub fn coefficient_truncate(
  table : CoefficientTable,
  maximum_degree : Int,
) -> CoefficientTable {
  let terms : Array[HarmonicTerm] = []
  for term in table.terms {
    if term.degree <= maximum_degree {
      terms.push(term)
    }
  }
  { ..table, terms, }
}

///|
pub fn coefficient_scale(
  table : CoefficientTable,
  factor : Double,
) -> CoefficientTable {
  let terms : Array[HarmonicTerm] = []
  for term in table.terms {
    terms.push({
      ..term,
      cosine: term.cosine * factor,
      sine: term.sine * factor,
    })
  }
  { ..table, scale: table.scale * factor, terms }
}

///|
pub fn coefficient_merge(
  left : CoefficientTable,
  right : CoefficientTable,
) -> CoefficientTable? {
  if left.reference_radius_m != right.reference_radius_m {
    None
  } else {
    let terms = left.terms.copy()
    for right_term in right.terms {
      let mut found = false
      for i in 0.. CoefficientTable {
  coefficient_truncate(
    { ..table, terms: table.terms.filter(fn(term) { term.degree != degree }) },
    coefficient_max_degree(table),
  )
}

///|
fn split_coefficient_line(line : String) -> Array[String] {
  let normalized = line
    .replace_all(old=";", new=" ")
    .replace_all(old=",", new=" ")
  let output : Array[String] = []
  for cell in normalized.split(" ") {
    let value = cell.trim()
    if !value.is_empty() {
      output.push(value.to_owned())
    }
  }
  output
}

///|
fn parse_int_cell(cell : StringView) -> Int raise {
  @strconv.from_str(cell.trim())
}

///|
/// Parse simple coefficient rows: degree order cosine sine.
pub fn parse_coefficient_table(
  text : String,
  reference_radius_m? : Double = MEAN_EARTH_RADIUS_M,
  scale? : Double = 1.0,
  source? : String = "text",
) -> CoefficientTable raise {
  let terms : Array[HarmonicTerm] = []
  let mut line_number = 0
  for raw in text.split("\n") {
    line_number += 1
    let line = raw.trim()
    if line.is_empty() || line.has_prefix("#") || line.has_prefix("degree") {
      ()
    } else {
      let cells = split_coefficient_line(line.to_owned())
      if cells.length() < 4 {
        fail(
          "coefficient line \{line_number}: expected degree order cosine sine",
        )
      }
      let degree = parse_int_cell(cells[0])
      let order = parse_int_cell(cells[1])
      let cosine = parse_double_cell(cells[2])
      let sine = parse_double_cell(cells[3])
      let term = { degree, order, cosine, sine }
      if !valid_term(term) {
        fail("coefficient line \{line_number}: invalid harmonic term")
      }
      terms.push(term)
    }
  }
  { reference_radius_m, scale, terms, source }
}

///|
pub fn coefficient_table_to_text(table : CoefficientTable) -> String {
  let rows : Array[String] = ["degree order cosine sine"]
  for term in table.terms {
    rows.push("\{term.degree} \{term.order} \{term.cosine} \{term.sine}")
  }
  rows.join("\n")
}

///|
pub fn evaluate_coefficient_table(
  table : CoefficientTable,
  point : GeoPoint,
) -> Double {
  evaluate_spherical_harmonic(coefficient_table_to_model(table), point)
}

///|
pub fn evaluate_coefficient_grid(table : CoefficientTable, grid : Grid) -> Grid {
  let values : Array[Double] = []
  for row in 0.. values.push(evaluate_coefficient_table(table, point))
        None => values.push(0.0 / 0.0)
      }
    }
  }
  { ..grid, values, }
}

///|
pub fn coefficient_roundtrip(
  table : CoefficientTable,
) -> CoefficientTable raise {
  parse_coefficient_table(
    coefficient_table_to_text(table),
    reference_radius_m=table.reference_radius_m,
    scale=table.scale,
    source=table.source,
  )
}

///|
fn find_term_at(
  table : CoefficientTable,
  degree : Int,
  order : Int,
) -> HarmonicTerm? {
  for term in table.terms {
    if term.degree == degree && term.order == order {
      return Some(term)
    }
  }
  None
}

///|
pub fn coefficient_difference(
  left : CoefficientTable,
  right : CoefficientTable,
) -> CoefficientTable? {
  if left.reference_radius_m != right.reference_radius_m {
    None
  } else {
    let terms : Array[HarmonicTerm] = []
    let max_degree = Int::max(
      coefficient_max_degree(left),
      coefficient_max_degree(right),
    )
    for degree in 0..<=max_degree {
      for order in 0..<=degree {
        let a = find_term_at(left, degree, order)
        let b = find_term_at(right, degree, order)
        match (a, b) {
          (Some(x), Some(y)) =>
            terms.push({
              degree,
              order,
              cosine: x.cosine - y.cosine,
              sine: x.sine - y.sine,
            })
          (Some(x), None) => terms.push(x)
          (None, Some(y)) =>
            terms.push({ ..y, cosine: -y.cosine, sine: -y.sine })
          (None, None) => ()
        }
      }
    }
    Some({
      reference_radius_m: left.reference_radius_m,
      scale: left.scale,
      terms,
      source: "difference",
    })
  }
}