///|
let private_log10 : Double = 2.302585092994046

///|
fn safe_log(value : Double) -> Double {
  @math.ln(value.max(1.0e-300))
}

///|
fn safe_log10(value : Double) -> Double {
  safe_log(value) / private_log10
}

///|
fn safe_exp(value : Double) -> Double {
  @math.exp(value.max(-700.0).min(700.0))
}

///|
fn activity_model_name(model : String) -> String? {
  match model {
    "Davies" | "davies" => Some("Davies")
    "DebyeHuckel" | "debye_huckel" | "debye-huckel" => Some("DebyeHuckel")
    _ => None
  }
}

///|
/// Calculate Davies or limiting Debye-Hückel activity coefficients.  The
/// returned array aligns with the species array, and neutral species have
/// coefficient one.  The implementation deliberately rejects the range where
/// these dilute-solution models cease to be a defensible MVP approximation.
fn calculate_activity_coefficients(
  model : String,
  temperature_k : Double,
  ionic_strength_value : Double,
  species : Array[SpeciesDef],
) -> Result[Array[Double], SolveError] {
  let canonical = match activity_model_name(model) {
    Some(name) => name
    None =>
      return Err(
        ActivityModelOutOfRange("activity_model must be Davies or DebyeHuckel"),
      )
  }
  if !finite(temperature_k) || (temperature_k - 298.15).abs() > 1.0e-8 {
    return Err(
      ActivityModelOutOfRange(
        "activity models and builtin constants support 298.15 K only",
      ),
    )
  }
  if !finite(ionic_strength_value) || ionic_strength_value < 0.0 {
    return Err(ActivityModelOutOfRange("ionic strength cannot be negative"))
  }
  let limit = if canonical == "Davies" { 0.1 } else { 0.01 }
  if ionic_strength_value > limit {
    return Err(
      ActivityModelOutOfRange(
        "ionic strength exceeds the selected dilute-solution activity model",
      ),
    )
  }
  // A at 25 C. No unsupported temperature correction is applied.
  let a = 0.5092
  let root = ionic_strength_value.sqrt()
  let coefficients : Array[Double] = []
  for species_def in species {
    let z = species_def.charge.to_double()
    let z2 = z * z
    let log_gamma = if z2 == 0.0 {
      0.0
    } else if canonical == "Davies" {
      -a * z2 * (root / (1.0 + root) - 0.3 * ionic_strength_value)
    } else {
      -a * z2 * root
    }
    coefficients.push(safe_exp(log_gamma * private_log10))
  }
  Ok(coefficients)
}

///|
fn ionic_strength_of(
  concentrations : Array[Double],
  species : Array[SpeciesDef],
) -> Double {
  let mut total = 0.0
  let count = concentrations.length().min(species.length())
  for i = 0; i < count; i = i + 1 {
    let z = species[i].charge.to_double()
    total = total + concentrations[i] * z * z
  }
  0.5 * total
}