///|
/// A station waveform prepared for array beamforming.
pub struct ArrayTrace {
  station : String
  position : StationPosition
  data : Array[Double]
  sample_rate : Double
  start_time : Double
}

///|
/// Local east/north coordinates relative to an array reference.
pub struct ArrayCoordinate {
  station : String
  east_km : Double
  north_km : Double
  elevation_km : Double
}

///|
/// Grid and beamforming parameters.
pub struct BeamConfiguration {
  slowness_min : Double
  slowness_max : Double
  slowness_step : Double
  window_size : Int
  hop_size : Int
  reference_velocity : Double
}

///|
/// One beamforming solution for a slowness vector.
pub struct BeamResult {
  slowness_east : Double
  slowness_north : Double
  back_azimuth : Double
  apparent_velocity : Double
  semblance : Double
  energy : Double
  peak_time : Double
  traces_used : Int
}

///|
/// A time series produced by stacking aligned station traces.
pub struct BeamTrace {
  start_time : Double
  sample_rate : Double
  data : Array[Double]
  delays : Array[Int]
  stations : Array[String]
}

///|
/// A rectangular slowness map with one score per grid point.
pub struct SlownessMap {
  east_values : Array[Double]
  north_values : Array[Double]
  semblance : Array[Double]
  energy : Array[Double]
  best_index : Int
}

///|
/// Default configuration for local microseismic arrays.
pub fn default_beam_configuration() -> BeamConfiguration {
  {
    slowness_min: -0.5,
    slowness_max: 0.5,
    slowness_step: 0.05,
    window_size: 128,
    hop_size: 32,
    reference_velocity: 4.0,
  }
}

///|
/// Convert latitude/longitude to a local east/north approximation.
pub fn local_array_coordinates(
  reference : StationPosition,
  stations : Array[StationPosition],
) -> Array[ArrayCoordinate] {
  let result = Array::new(capacity=stations.length())
  let latitude_scale = 111.32
  let longitude_scale = 111.32 *
    @math.cos(reference.latitude * @math.PI / 180.0)
  for station in stations {
    result.push({
      station: station.station,
      east_km: (station.longitude - reference.longitude) * longitude_scale,
      north_km: (station.latitude - reference.latitude) * latitude_scale,
      elevation_km: station.elevation_km - reference.elevation_km,
    })
  }
  result
}

///|
/// Create a coordinate list using the first station as the reference.
pub fn coordinates_from_stations(
  stations : Array[StationPosition],
) -> Array[ArrayCoordinate] {
  if stations.is_empty() {
    return []
  }
  local_array_coordinates(stations[0], stations)
}

///|
/// Return the station coordinate matching an identifier.
pub fn find_array_coordinate(
  coordinates : Array[ArrayCoordinate],
  station : String,
) -> ArrayCoordinate? {
  for coordinate in coordinates {
    if coordinate.station == station {
      return Some(coordinate)
    }
  }
  None
}

///|
/// Calculate integer sample delays for a slowness vector.
pub fn beam_delays(
  coordinates : Array[ArrayCoordinate],
  slowness_east : Double,
  slowness_north : Double,
  sample_rate : Double,
) -> Array[Int] {
  let result = Array::new(capacity=coordinates.length())
  let rate = sample_rate.max(1.0e-12)
  for coordinate in coordinates {
    let seconds = coordinate.east_km * slowness_east +
      coordinate.north_km * slowness_north
    result.push((seconds * rate).round().to_int())
  }
  result
}

///|
/// Stack traces after applying integer station delays.
pub fn stack_delayed_traces(
  traces : Array[Array[Double]],
  delays : Array[Int],
) -> Array[Double] {
  if traces.is_empty() {
    return []
  }
  let mut count = traces[0].length()
  for data in traces {
    if data.length() < count {
      count = data.length()
    }
  }
  let result = Array::make(count, 0.0)
  let used = Array::make(count, 0)
  let station_count = traces.length().min(delays.length())
  for station = 0; station < station_count; station = station + 1 {
    let data = traces[station]
    let delay = delays[station]
    for index = 0; index < count; index = index + 1 {
      let source = index - delay
      if source >= 0 && source < data.length() {
        result[index] += data[source]
        used[index] += 1
      }
    }
  }
  for index = 0; index < count; index = index + 1 {
    if used[index] > 0 {
      result[index] /= used[index].to_double()
    }
  }
  result
}

///|
/// Calculate a semblance score for aligned station samples.
pub fn beam_semblance(
  traces : Array[Array[Double]],
  delays : Array[Int],
) -> Double {
  if traces.is_empty() {
    return 0.0
  }
  let mut count = traces[0].length()
  for data in traces {
    if data.length() < count {
      count = data.length()
    }
  }
  let stations = traces.length().min(delays.length())
  if count == 0 || stations == 0 {
    return 0.0
  }
  let mut numerator = 0.0
  let mut denominator = 0.0
  for index = 0; index < count; index = index + 1 {
    let mut sum = 0.0
    let mut energy = 0.0
    let mut used = 0
    for station = 0; station < stations; station = station + 1 {
      let source = index - delays[station]
      if source >= 0 && source < traces[station].length() {
        let value = traces[station][source]
        sum += value
        energy += value * value
        used += 1
      }
    }
    numerator += sum * sum
    if used > 0 {
      denominator += used.to_double() * energy
    }
  }
  if denominator <= 0.0 {
    0.0
  } else {
    numerator / denominator
  }
}

///|
/// Return the time of the largest absolute sample in a beam.
pub fn beam_peak_time(
  data : Array[Double],
  sample_rate : Double,
  start_time : Double,
) -> Double {
  if data.is_empty() || sample_rate <= 0.0 {
    return start_time
  }
  let mut best = 0
  for index = 1; index < data.length(); index = index + 1 {
    if data[index].abs() > data[best].abs() {
      best = index
    }
  }
  start_time + best.to_double() / sample_rate
}

///|
/// Build a beam result from a slowness vector and station traces.
pub fn evaluate_beam(
  traces : Array[ArrayTrace],
  coordinates : Array[ArrayCoordinate],
  slowness_east : Double,
  slowness_north : Double,
) -> (BeamResult, BeamTrace) {
  let data = traces.map(trace => trace.data)
  let sample_rate = if traces.is_empty() { 0.0 } else { traces[0].sample_rate }
  let start_time = if traces.is_empty() { 0.0 } else { traces[0].start_time }
  let delays = beam_delays(
    coordinates, slowness_east, slowness_north, sample_rate,
  )
  let stacked = stack_delayed_traces(data, delays)
  let semblance = beam_semblance(data, delays)
  let magnitude = (slowness_east * slowness_east +
  slowness_north * slowness_north).sqrt()
  let velocity = if magnitude <= 1.0e-12 { 1.0e300 } else { 1.0 / magnitude }
  let azimuth = @math.atan2(slowness_east, slowness_north) * 180.0 / @math.PI
  let normalized_azimuth = if azimuth < 0.0 { azimuth + 360.0 } else { azimuth }
  let result = {
    slowness_east,
    slowness_north,
    back_azimuth: normalized_azimuth,
    apparent_velocity: velocity,
    semblance,
    energy: signal_energy(stacked),
    peak_time: beam_peak_time(stacked, sample_rate, start_time),
    traces_used: traces.length().min(coordinates.length()),
  }
  let station_names = traces.map(trace => trace.station)
  (
    result,
    { start_time, sample_rate, data: stacked, delays, stations: station_names },
  )
}

///|
/// Return a linearly spaced slowness axis including both endpoints.
pub fn slowness_axis(
  minimum : Double,
  maximum : Double,
  step : Double,
) -> Array[Double] {
  let result = Array::new()
  if step <= 0.0 || maximum < minimum {
    return result
  }
  let mut value = minimum
  while value <= maximum + step * 0.001 {
    result.push(value)
    value += step
  }
  result
}

///|
/// Scan a rectangular slowness grid and keep its best solution.
pub fn scan_slowness_grid(
  traces : Array[ArrayTrace],
  coordinates : Array[ArrayCoordinate],
  config : BeamConfiguration,
) -> SlownessMap {
  let east_values = slowness_axis(
    config.slowness_min,
    config.slowness_max,
    config.slowness_step,
  )
  let north_values = east_values.copy()
  let semblance = Array::make(east_values.length() * north_values.length(), 0.0)
  let energy = Array::make(semblance.length(), 0.0)
  let mut best_index = -1
  let mut best_score = -1.0
  for north_index = 0
      north_index < north_values.length()
      north_index = north_index + 1 {
    for east_index = 0
        east_index < east_values.length()
        east_index = east_index + 1 {
      let index = north_index * east_values.length() + east_index
      let (beam, _) = evaluate_beam(
        traces,
        coordinates,
        east_values[east_index],
        north_values[north_index],
      )
      semblance[index] = beam.semblance
      energy[index] = beam.energy
      if beam.semblance > best_score {
        best_score = beam.semblance
        best_index = index
      }
    }
  }
  { east_values, north_values, semblance, energy, best_index }
}

///|
/// Convert the best map cell into a beam result.
pub fn best_map_beam(
  map : SlownessMap,
  traces : Array[ArrayTrace],
  coordinates : Array[ArrayCoordinate],
) -> BeamResult? {
  if map.best_index < 0 ||
    map.east_values.is_empty() ||
    map.north_values.is_empty() {
    return None
  }
  let width = map.east_values.length()
  let north_index = map.best_index / width
  let east_index = map.best_index % width
  if north_index >= map.north_values.length() {
    return None
  }
  let (result, _) = evaluate_beam(
    traces,
    coordinates,
    map.east_values[east_index],
    map.north_values[north_index],
  )
  Some(result)
}

///|
/// Produce a sliding beam time series for a fixed slowness vector.
pub fn beam_time_series(
  traces : Array[ArrayTrace],
  coordinates : Array[ArrayCoordinate],
  slowness_east : Double,
  slowness_north : Double,
  window_size : Int,
  hop_size : Int,
) -> Array[BeamResult] {
  let result = Array::new()
  if traces.is_empty() || window_size <= 0 || hop_size <= 0 {
    return result
  }
  let mut count = traces[0].data.length()
  for trace in traces {
    if trace.data.length() < count {
      count = trace.data.length()
    }
  }
  let mut start = 0
  while start + window_size <= count {
    let windowed = traces.map(trace => {
      station: trace.station,
      position: trace.position,
      data: slice_samples(trace.data, start, start + window_size),
      sample_rate: trace.sample_rate,
      start_time: trace.start_time +
      start.to_double() / trace.sample_rate.max(1.0e-12),
    })
    let (beam, _) = evaluate_beam(
      windowed, coordinates, slowness_east, slowness_north,
    )
    result.push(beam)
    start += hop_size
  }
  result
}

///|
/// Compute a detection score that combines semblance and beam energy.
pub fn beam_detection_score(
  result : BeamResult,
  noise_energy : Double,
) -> Double {
  let noise = noise_energy.max(1.0e-20)
  result.semblance.max(0.0) * @math.log10(1.0 + result.energy / noise)
}

///|
/// Return a map as rows suitable for contour or heatmap rendering.
pub fn slowness_map_to_csv(map : SlownessMap) -> String {
  let mut output = "east_slowness,north_slowness,semblance,energy\n"
  let width = map.east_values.length()
  for index = 0; index < map.semblance.length(); index = index + 1 {
    let north_index = if width == 0 { 0 } else { index / width }
    let east_index = if width == 0 { 0 } else { index % width }
    if north_index < map.north_values.length() && east_index < width {
      output = output +
        map.east_values[east_index].to_string() +
        "," +
        map.north_values[north_index].to_string() +
        "," +
        map.semblance[index].to_string() +
        "," +
        map.energy[index].to_string() +
        "\n"
    }
  }
  output
}

///|
/// Validate that traces can participate in a common array stack.
pub fn validate_array_traces(
  traces : Array[ArrayTrace],
) -> Result[Unit, String] {
  if traces.is_empty() {
    return Err("Array requires at least one station trace")
  }
  let reference_rate = traces[0].sample_rate
  let reference_time = traces[0].start_time
  for trace in traces {
    if trace.sample_rate <= 0.0 {
      return Err("Array trace has a non-positive sample rate")
    }
    if (trace.sample_rate - reference_rate).abs() > 1.0e-9 {
      return Err("Array traces use different sample rates")
    }
    if (trace.start_time - reference_time).abs() > 1.0 / reference_rate {
      return Err("Array traces are not approximately aligned")
    }
  }
  Ok(())
}

///|
/// Estimate a back-azimuth from a best beam result.
pub fn beam_back_azimuth(result : BeamResult) -> Double {
  result.back_azimuth
}

///|
/// Format a beam result for a structured processing log.
pub fn beam_result_text(result : BeamResult) -> String {
  "slow_east=" +
  result.slowness_east.to_string() +
  " slow_north=" +
  result.slowness_north.to_string() +
  " back_azimuth=" +
  result.back_azimuth.to_string() +
  " velocity=" +
  result.apparent_velocity.to_string() +
  " semblance=" +
  result.semblance.to_string() +
  " energy=" +
  result.energy.to_string() +
  " peak_time=" +
  result.peak_time.to_string() +
  " traces=" +
  result.traces_used.to_string()
}