///|
/// Configuration for three-component polarization analysis.
pub struct PolarizationConfig {
  window_size : Int
  step : Int
  iterations : Int
  noise_floor : Double
}

///|
/// Eigenvector and eigenvalue summary of a covariance window.
pub struct PolarizationAttributes {
  eigenvalue_major : Double
  eigenvalue_middle : Double
  eigenvalue_minor : Double
  azimuth_degrees : Double
  incidence_degrees : Double
  rectilinearity : Double
  planarity : Double
  ellipticity : Double
  energy : Double
  confidence : Double
}

///|
/// Polarization attributes associated with a physical time window.
pub struct PolarizationWindow {
  start_index : Int
  end_index : Int
  start_time : Double
  end_time : Double
  attributes : PolarizationAttributes
}

///|
/// A compact summary of a polarization track.
pub struct PolarizationSummary {
  windows : Int
  mean_azimuth : Double
  mean_incidence : Double
  mean_rectilinearity : Double
  mean_planarity : Double
  maximum_energy : Double
  most_linear_window : Int
}

///|
/// Create a practical default for local P/S phase analysis.
pub fn default_polarization_config() -> PolarizationConfig {
  { window_size: 64, step: 16, iterations: 24, noise_floor: 1.0e-12 }
}

///|
/// Return a normalized covariance matrix in row-major order.
pub fn triplet_covariance_matrix(triplet : ComponentTriplet) -> Array[Double] {
  let terms = component_covariance(triplet)
  [
    terms[0],
    terms[3],
    terms[4],
    terms[3],
    terms[1],
    terms[5],
    terms[4],
    terms[5],
    terms[2],
  ]
}

///|
/// Calculate a 3 by 3 matrix-vector product.
fn matrix_vector(
  matrix : Array[Double],
  vector : Array[Double],
) -> Array[Double] {
  [
    matrix[0] * vector[0] + matrix[1] * vector[1] + matrix[2] * vector[2],
    matrix[3] * vector[0] + matrix[4] * vector[1] + matrix[5] * vector[2],
    matrix[6] * vector[0] + matrix[7] * vector[1] + matrix[8] * vector[2],
  ]
}

///|
/// Normalize a three-dimensional vector, returning zero for a zero vector.
fn normalize_vector(vector : Array[Double]) -> Array[Double] {
  let norm = (vector[0] * vector[0] +
  vector[1] * vector[1] +
  vector[2] * vector[2]).sqrt()
  if norm <= 0.0 {
    [0.0, 0.0, 0.0]
  } else {
    [vector[0] / norm, vector[1] / norm, vector[2] / norm]
  }
}

///|
/// Estimate the dominant eigenvector by power iteration.
pub fn dominant_covariance_vector(
  matrix : Array[Double],
  iterations : Int,
) -> Array[Double] {
  if matrix.length() < 9 {
    return [0.0, 0.0, 0.0]
  }
  let count = if iterations <= 0 { 1 } else { iterations }
  let mut vector = normalize_vector([1.0, 0.5, 0.25])
  for _ in 0.. Double {
  let projected = matrix_vector(matrix, vector)
  vector[0] * projected[0] + vector[1] * projected[1] + vector[2] * projected[2]
}

///|
/// Remove a dominant eigenvector from a covariance matrix by deflation.
fn deflate_covariance(
  matrix : Array[Double],
  vector : Array[Double],
  eigenvalue : Double,
) -> Array[Double] {
  let result = matrix.copy()
  for row = 0; row < 3; row = row + 1 {
    for column = 0; column < 3; column = column + 1 {
      result[row * 3 + column] -= eigenvalue * vector[row] * vector[column]
    }
  }
  result
}

///|
/// Estimate the three non-negative covariance eigenvalues.
pub fn covariance_eigenvalues(
  matrix : Array[Double],
  iterations : Int,
) -> Array[Double] {
  if matrix.length() < 9 {
    return [0.0, 0.0, 0.0]
  }
  let major_vector = dominant_covariance_vector(matrix, iterations)
  let major = covariance_rayleigh(matrix, major_vector).max(0.0)
  let residual = deflate_covariance(matrix, major_vector, major)
  let middle_vector = dominant_covariance_vector(residual, iterations)
  let middle = covariance_rayleigh(residual, middle_vector).max(0.0)
  let minor = (matrix[0] + matrix[4] + matrix[8] - major - middle).max(0.0)
  let values = [major, middle, minor]
  values.sort()
  values.rev_in_place()
  values
}

///|
/// Convert a dominant component vector to azimuth and incidence angles.
pub fn polarization_angles(vector : Array[Double]) -> (Double, Double) {
  if vector.length() < 3 {
    return (0.0, 0.0)
  }
  let azimuth = @math.atan2(vector[1], vector[0]) * 180.0 / @math.PI
  let horizontal = (vector[0] * vector[0] + vector[1] * vector[1]).sqrt()
  let incidence = @math.atan2(horizontal, vector[2].abs()) * 180.0 / @math.PI
  let normalized_azimuth = if azimuth < 0.0 { azimuth + 360.0 } else { azimuth }
  (normalized_azimuth, incidence)
}

///|
/// Compute standard polarization attributes from one component window.
pub fn polarization_attributes(
  triplet : ComponentTriplet,
  iterations : Int,
  noise_floor : Double,
) -> PolarizationAttributes {
  let matrix = triplet_covariance_matrix(triplet)
  let vector = dominant_covariance_vector(matrix, iterations)
  let values = covariance_eigenvalues(matrix, iterations)
  let major = values[0]
  let middle = values[1]
  let minor = values[2]
  let (azimuth, incidence) = polarization_angles(vector)
  let energy = matrix[0] + matrix[4] + matrix[8]
  let safe_major = major.max(noise_floor.max(1.0e-20))
  let rectilinearity = (1.0 - middle / safe_major).max(0.0).min(1.0)
  let planarity = if major <= 0.0 { 0.0 } else { 1.0 - minor / major }
  let ellipticity = if major <= 0.0 { 0.0 } else { (middle / major).sqrt() }
  let confidence = if energy <= noise_floor.max(1.0e-20) {
    0.0
  } else {
    (major / energy).max(0.0).min(1.0)
  }
  {
    eigenvalue_major: major,
    eigenvalue_middle: middle,
    eigenvalue_minor: minor,
    azimuth_degrees: azimuth,
    incidence_degrees: incidence,
    rectilinearity,
    planarity: planarity.max(0.0).min(1.0),
    ellipticity,
    energy,
    confidence,
  }
}

///|
/// Build a component window from a triplet and a half-open sample range.
pub fn triplet_window(
  triplet : ComponentTriplet,
  start : Int,
  end : Int,
) -> ComponentTriplet {
  {
    station: triplet.station,
    sampling_rate: triplet.sampling_rate,
    start_time: triplet.start_time +
    start.max(0).to_double() / triplet.sampling_rate.max(1.0e-12),
    north: slice_samples(triplet.north, start, end),
    east: slice_samples(triplet.east, start, end),
    vertical: slice_samples(triplet.vertical, start, end),
  }
}

///|
/// Analyze every complete component window.
pub fn analyze_polarization(
  triplet : ComponentTriplet,
  config : PolarizationConfig,
) -> Array[PolarizationWindow] {
  let result = Array::new()
  let count = triplet.north
    .length()
    .min(triplet.east.length())
    .min(triplet.vertical.length())
  let size = config.window_size
  let step = if config.step <= 0 { size } else { config.step }
  if count == 0 || size <= 0 || step <= 0 {
    return result
  }
  let mut start = 0
  while start + size <= count {
    let end = start + size
    let window = triplet_window(triplet, start, end)
    let attributes = polarization_attributes(
      window,
      config.iterations,
      config.noise_floor,
    )
    result.push({
      start_index: start,
      end_index: end,
      start_time: triplet.start_time +
      start.to_double() / triplet.sampling_rate.max(1.0e-12),
      end_time: triplet.start_time +
      (end - 1).to_double() / triplet.sampling_rate.max(1.0e-12),
      attributes,
    })
    start += step
  }
  result
}

///|
/// Find the window most consistent with a rectilinear P-wave motion.
pub fn best_linear_polarization(
  windows : Array[PolarizationWindow],
) -> PolarizationWindow? {
  if windows.is_empty() {
    return None
  }
  let mut best = 0
  for index = 1; index < windows.length(); index = index + 1 {
    let current = windows[index].attributes.rectilinearity *
      windows[index].attributes.confidence
    let previous = windows[best].attributes.rectilinearity *
      windows[best].attributes.confidence
    if current > previous {
      best = index
    }
  }
  Some(windows[best])
}

///|
/// Compute a robust circular mean of azimuths weighted by confidence.
pub fn weighted_azimuth_mean(windows : Array[PolarizationWindow]) -> Double {
  let mut x = 0.0
  let mut y = 0.0
  for window in windows {
    let angle = window.attributes.azimuth_degrees * @math.PI / 180.0
    let weight = window.attributes.confidence.max(0.0)
    x += @math.cos(angle) * weight
    y += @math.sin(angle) * weight
  }
  if x.abs() + y.abs() <= 1.0e-20 {
    0.0
  } else {
    let angle = @math.atan2(y, x) * 180.0 / @math.PI
    if angle < 0.0 {
      angle + 360.0
    } else {
      angle
    }
  }
}

///|
/// Summarize a polarization track for event reports.
pub fn summarize_polarization(
  windows : Array[PolarizationWindow],
) -> PolarizationSummary {
  let mut incidence = 0.0
  let mut rectilinearity = 0.0
  let mut planarity = 0.0
  let mut energy = 0.0
  let mut best_score = -1.0
  let mut best_index = -1
  for index = 0; index < windows.length(); index = index + 1 {
    let attributes = windows[index].attributes
    incidence += attributes.incidence_degrees
    rectilinearity += attributes.rectilinearity
    planarity += attributes.planarity
    if attributes.energy > energy {
      energy = attributes.energy
    }
    let score = attributes.rectilinearity * attributes.confidence
    if score > best_score {
      best_score = score
      best_index = index
    }
  }
  let divisor = windows.length().to_double().max(1.0)
  {
    windows: windows.length(),
    mean_azimuth: weighted_azimuth_mean(windows),
    mean_incidence: incidence / divisor,
    mean_rectilinearity: rectilinearity / divisor,
    mean_planarity: planarity / divisor,
    maximum_energy: energy,
    most_linear_window: best_index,
  }
}

///|
/// Compute the covariance trace of one component window.
pub fn polarization_energy(triplet : ComponentTriplet) -> Double {
  let matrix = triplet_covariance_matrix(triplet)
  matrix[0] + matrix[4] + matrix[8]
}

///|
/// Return a scalar score for a P-wave candidate.
pub fn p_wave_polarization_score(attributes : PolarizationAttributes) -> Double {
  attributes.rectilinearity *
  attributes.confidence *
  (1.0 - attributes.incidence_degrees / 180.0).max(0.0)
}

///|
/// Return a scalar score for an S-wave candidate.
pub fn s_wave_polarization_score(attributes : PolarizationAttributes) -> Double {
  attributes.planarity * attributes.ellipticity.max(0.0) * attributes.confidence
}

///|
/// Convert a window track into a stable CSV table.
pub fn polarization_windows_to_csv(
  windows : Array[PolarizationWindow],
) -> String {
  let mut output = "start_index,end_index,start_time,end_time,azimuth,incidence,rectilinearity,planarity,ellipticity,energy,confidence\n"
  for window in windows {
    let a = window.attributes
    output = output +
      window.start_index.to_string() +
      "," +
      window.end_index.to_string() +
      "," +
      window.start_time.to_string() +
      "," +
      window.end_time.to_string() +
      "," +
      a.azimuth_degrees.to_string() +
      "," +
      a.incidence_degrees.to_string() +
      "," +
      a.rectilinearity.to_string() +
      "," +
      a.planarity.to_string() +
      "," +
      a.ellipticity.to_string() +
      "," +
      a.energy.to_string() +
      "," +
      a.confidence.to_string() +
      "\n"
  }
  output
}

///|
/// Return whether a polarization window is suitable for automatic picking.
pub fn polarization_is_confident(
  attributes : PolarizationAttributes,
  minimum_confidence : Double,
  minimum_energy : Double,
) -> Bool {
  attributes.confidence >= minimum_confidence &&
  attributes.energy >= minimum_energy
}

///|
/// Blend two component windows with a confidence-based weight.
pub fn blend_component_windows(
  left : ComponentTriplet,
  right : ComponentTriplet,
  right_weight : Double,
) -> ComponentTriplet {
  let count = left.north
    .length()
    .min(right.north.length())
    .min(left.east.length())
    .min(right.east.length())
    .min(left.vertical.length())
    .min(right.vertical.length())
  let weight = right_weight.max(0.0).min(1.0)
  let north = Array::make(count, 0.0)
  let east = Array::make(count, 0.0)
  let vertical = Array::make(count, 0.0)
  for index = 0; index < count; index = index + 1 {
    north[index] = left.north[index] * (1.0 - weight) +
      right.north[index] * weight
    east[index] = left.east[index] * (1.0 - weight) + right.east[index] * weight
    vertical[index] = left.vertical[index] * (1.0 - weight) +
      right.vertical[index] * weight
  }
  {
    station: left.station,
    sampling_rate: left.sampling_rate,
    start_time: left.start_time,
    north,
    east,
    vertical,
  }
}

///|
/// Format a polarization summary for an operator-facing report.
pub fn polarization_summary_text(summary : PolarizationSummary) -> String {
  "windows=" +
  summary.windows.to_string() +
  " mean_azimuth=" +
  summary.mean_azimuth.to_string() +
  " mean_incidence=" +
  summary.mean_incidence.to_string() +
  " rectilinearity=" +
  summary.mean_rectilinearity.to_string() +
  " planarity=" +
  summary.mean_planarity.to_string() +
  " maximum_energy=" +
  summary.maximum_energy.to_string() +
  " best_window=" +
  summary.most_linear_window.to_string()
}