///|
/// Validation diagnostics for a one-dimensional velocity model.
pub struct VelocityModelDiagnostics {
  valid : Bool
  layer_count : Int
  minimum_vp : Double
  maximum_vp : Double
  minimum_vs : Double
  maximum_vs : Double
  discontinuities : Int
  warnings : Array[String]
}

///|
/// One sample along a straight ray path.
pub struct RayPathSample {
  fraction : Double
  depth_km : Double
  distance_km : Double
  vp : Double
  vs : Double
  p_slowness : Double
  s_slowness : Double
}

///|
/// Travel-time table for a depth/distance grid.
pub struct TravelTimeTable {
  depths_km : Array[Double]
  distances_km : Array[Double]
  p_times : Array[Double]
  s_times : Array[Double]
  width : Int
}

///|
/// Return default depth samples used by diagnostics.
pub fn standard_depth_samples(
  maximum_depth : Double,
  count : Int,
) -> Array[Double] {
  let result = Array::new()
  if maximum_depth < 0.0 || count <= 0 {
    return result
  }
  if count == 1 {
    result.push(0.0)
    return result
  }
  for index = 0; index < count; index = index + 1 {
    result.push(maximum_depth * index.to_double() / (count - 1).to_double())
  }
  result
}

///|
/// Validate ordering and physical ranges of velocity layers.
pub fn validate_velocity_model_detailed(
  model : VelocityModel,
) -> VelocityModelDiagnostics {
  let warnings : Array[String] = Array::new()
  if model.layers.is_empty() {
    warnings.push("model has no layers")
    return {
      valid: false,
      layer_count: 0,
      minimum_vp: 0.0,
      maximum_vp: 0.0,
      minimum_vs: 0.0,
      maximum_vs: 0.0,
      discontinuities: 0,
      warnings,
    }
  }
  let mut min_vp = model.layers[0].vp
  let mut max_vp = model.layers[0].vp
  let mut min_vs = model.layers[0].vs
  let mut max_vs = model.layers[0].vs
  let mut discontinuities = 0
  let mut valid = true
  let mut previous_depth = -1.0e300
  for layer in model.layers {
    min_vp = min_vp.min(layer.vp)
    max_vp = max_vp.max(layer.vp)
    min_vs = min_vs.min(layer.vs)
    max_vs = max_vs.max(layer.vs)
    if layer.depth < previous_depth {
      valid = false
      warnings.push("layer depths are not sorted")
    }
    if layer.vp <= 0.0 || layer.vs <= 0.0 || layer.vs > layer.vp {
      valid = false
      warnings.push("layer contains invalid wave speeds")
    }
    if layer.depth > previous_depth && previous_depth > -1.0e200 {
      discontinuities += 1
    }
    previous_depth = layer.depth
  }
  {
    valid,
    layer_count: model.layers.length(),
    minimum_vp: min_vp,
    maximum_vp: max_vp,
    minimum_vs: min_vs,
    maximum_vs: max_vs,
    discontinuities,
    warnings,
  }
}

///|
/// Linearly interpolate a velocity between neighboring layers.
pub fn interpolated_velocity(
  model : VelocityModel,
  depth : Double,
  is_s_wave : Bool,
) -> Double {
  if model.layers.is_empty() {
    return 0.0
  }
  if model.layers.length() == 1 {
    return get_velocity_at_depth(model, depth, is_s_wave)
  }
  let first = model.layers[0]
  if depth <= first.depth {
    return if is_s_wave { first.vs } else { first.vp }
  }
  for index = 1; index < model.layers.length(); index = index + 1 {
    let left = model.layers[index - 1]
    let right = model.layers[index]
    if depth <= right.depth {
      let denominator = (right.depth - left.depth).max(1.0e-12)
      let fraction = (depth - left.depth) / denominator
      let left_velocity = if is_s_wave { left.vs } else { left.vp }
      let right_velocity = if is_s_wave { right.vs } else { right.vp }
      return left_velocity + fraction * (right_velocity - left_velocity)
    }
  }
  get_velocity_at_depth(model, depth, is_s_wave)
}

///|
/// Sample a straight source-to-receiver ray at a fixed resolution.
pub fn sample_ray_path(
  model : VelocityModel,
  source_depth : Double,
  distance : Double,
  samples : Int,
) -> Array[RayPathSample] {
  let result = Array::new()
  let count = samples.max(1)
  let path_length = (source_depth * source_depth + distance * distance).sqrt()
  for index = 0; index <= count; index = index + 1 {
    let fraction = index.to_double() / count.to_double()
    let depth = source_depth * (1.0 - fraction)
    let horizontal = distance * fraction
    let vp = interpolated_velocity(model, depth, false).max(1.0e-12)
    let vs = interpolated_velocity(model, depth, true).max(1.0e-12)
    result.push({
      fraction,
      depth_km: depth,
      distance_km: horizontal,
      vp,
      vs,
      p_slowness: 1.0 / vp,
      s_slowness: 1.0 / vs,
    })
  }
  ignore(path_length)
  result
}

///|
/// Integrate a sampled ray's slowness with the trapezoidal rule.
pub fn ray_path_travel_time(
  model : VelocityModel,
  source_depth : Double,
  distance : Double,
  is_s_wave : Bool,
  samples : Int,
) -> Double {
  let path_length = (source_depth * source_depth + distance * distance).sqrt()
  if path_length <= 0.0 {
    return 0.0
  }
  let path = sample_ray_path(model, source_depth, distance, samples)
  let count = path.length()
  if count < 2 {
    let velocity = if is_s_wave {
      get_velocity_at_depth(model, source_depth, true)
    } else {
      get_velocity_at_depth(model, source_depth, false)
    }
    return path_length / velocity.max(1.0e-12)
  }
  let mut total = 0.0
  for index = 1; index < count; index = index + 1 {
    let left = if is_s_wave {
      path[index - 1].s_slowness
    } else {
      path[index - 1].p_slowness
    }
    let right = if is_s_wave {
      path[index].s_slowness
    } else {
      path[index].p_slowness
    }
    total += 0.5 * (left + right)
  }
  total * path_length / (count - 1).to_double()
}

///|
/// Calculate a travel-time gradient with respect to source depth.
pub fn travel_time_depth_gradient(
  model : VelocityModel,
  source_depth : Double,
  distance : Double,
  is_s_wave : Bool,
) -> Double {
  let step = 0.01.max(source_depth.abs() * 0.001)
  let lower = travel_time(
    model,
    (source_depth - step).max(0.0),
    distance,
    is_s_wave,
  )
  let upper = travel_time(model, source_depth + step, distance, is_s_wave)
  (upper - lower) / (2.0 * step)
}

///|
/// Calculate a travel-time gradient with respect to epicentral distance.
pub fn travel_time_distance_gradient(
  model : VelocityModel,
  source_depth : Double,
  distance : Double,
  is_s_wave : Bool,
) -> Double {
  let step = 0.01.max(distance.abs() * 0.001)
  let lower = travel_time(
    model,
    source_depth,
    (distance - step).max(0.0),
    is_s_wave,
  )
  let upper = travel_time(model, source_depth, distance + step, is_s_wave)
  (upper - lower) / (2.0 * step)
}

///|
/// Build a rectangular P/S travel-time table.
pub fn make_travel_time_table(
  model : VelocityModel,
  depths : Array[Double],
  distances : Array[Double],
) -> TravelTimeTable {
  let p_times = Array::make(depths.length() * distances.length(), 0.0)
  let s_times = Array::make(p_times.length(), 0.0)
  for depth_index = 0
      depth_index < depths.length()
      depth_index = depth_index + 1 {
    for distance_index = 0
        distance_index < distances.length()
        distance_index = distance_index + 1 {
      let index = depth_index * distances.length() + distance_index
      p_times[index] = ray_path_travel_time(
        model,
        depths[depth_index],
        distances[distance_index],
        false,
        80,
      )
      s_times[index] = ray_path_travel_time(
        model,
        depths[depth_index],
        distances[distance_index],
        true,
        80,
      )
    }
  }
  {
    depths_km: depths.copy(),
    distances_km: distances.copy(),
    p_times,
    s_times,
    width: distances.length(),
  }
}

///|
/// Interpolate a value in a precomputed travel-time table.
pub fn lookup_travel_time(
  table : TravelTimeTable,
  depth_index : Int,
  distance_index : Int,
  is_s_wave : Bool,
) -> Double {
  if depth_index < 0 || distance_index < 0 || table.width <= 0 {
    return 0.0
  }
  let index = depth_index * table.width + distance_index
  let values = if is_s_wave { table.s_times } else { table.p_times }
  if index < 0 || index >= values.length() {
    0.0
  } else {
    values[index]
  }
}

///|
/// Return model diagnostics in a stable text format.
pub fn velocity_model_diagnostics_text(
  diagnostics : VelocityModelDiagnostics,
) -> String {
  "valid=" +
  diagnostics.valid.to_string() +
  " layers=" +
  diagnostics.layer_count.to_string() +
  " vp_min=" +
  diagnostics.minimum_vp.to_string() +
  " vp_max=" +
  diagnostics.maximum_vp.to_string() +
  " vs_min=" +
  diagnostics.minimum_vs.to_string() +
  " vs_max=" +
  diagnostics.maximum_vs.to_string() +
  " discontinuities=" +
  diagnostics.discontinuities.to_string() +
  " warnings=" +
  diagnostics.warnings.length().to_string()
}

///|
/// Return a CSV table of ray-path samples.
pub fn ray_path_to_csv(path : Array[RayPathSample]) -> String {
  let mut output = "fraction,depth_km,distance_km,vp,vs,p_slowness,s_slowness\n"
  for sample in path {
    output = output +
      sample.fraction.to_string() +
      "," +
      sample.depth_km.to_string() +
      "," +
      sample.distance_km.to_string() +
      "," +
      sample.vp.to_string() +
      "," +
      sample.vs.to_string() +
      "," +
      sample.p_slowness.to_string() +
      "," +
      sample.s_slowness.to_string() +
      "\n"
  }
  output
}