///|
/// Common strong-motion measures derived from one acceleration trace.
pub struct GroundMotionMetrics {
  pga : Double
  pgv : Double
  pgd : Double
  rms_acceleration : Double
  arias_intensity : Double
  cav : Double
  significant_duration : Double
  bracketed_duration : Double
  dominant_period : Double
}

///|
/// A response-spectrum ordinate at one oscillator period.
pub struct ResponseSpectrumPoint {
  period_seconds : Double
  damping : Double
  spectral_acceleration : Double
  spectral_velocity : Double
  spectral_displacement : Double
}

///|
/// A response spectrum with stable period ordering.
pub struct ResponseSpectrum {
  points : Array[ResponseSpectrumPoint]
  sample_rate : Double
  damping : Double
  peak_period : Double
  peak_acceleration : Double
}

///|
/// Duration measures based on cumulative Arias intensity.
pub struct DurationMeasures {
  start_time : Double
  end_time : Double
  bracketed_duration : Double
  significant_start : Double
  significant_end : Double
  significant_duration : Double
}

///|
/// A simple macroseismic intensity proxy from acceleration and distance.
pub struct IntensityEstimate {
  pga_g : Double
  hypocentral_distance_km : Double
  intensity : Double
  label : String
}

///|
/// Return the largest absolute amplitude.
pub fn absolute_peak(data : Array[Double]) -> Double {
  maximum_sample(data).abs().max(minimum_sample(data).abs())
}

///|
/// Integrate acceleration with a configurable baseline correction.
pub fn integrate_acceleration(
  acceleration : Array[Double],
  sample_rate : Double,
  remove_drift : Bool,
) -> Array[Double] {
  let corrected = if remove_drift {
    detrend_copy(acceleration)
  } else {
    acceleration.copy()
  }
  cumulative_integral(corrected, sample_rate)
}

///|
/// Integrate a velocity trace into displacement.
pub fn integrate_velocity(
  velocity : Array[Double],
  sample_rate : Double,
  remove_drift : Bool,
) -> Array[Double] {
  let corrected = if remove_drift {
    detrend_copy(velocity)
  } else {
    velocity.copy()
  }
  cumulative_integral(corrected, sample_rate)
}

///|
/// Compute cumulative absolute velocity.
pub fn cumulative_absolute_velocity(
  acceleration : Array[Double],
  sample_rate : Double,
) -> Double {
  if acceleration.is_empty() || sample_rate <= 0.0 {
    0.0
  } else {
    let dt = 1.0 / sample_rate
    let mut total = 0.0
    for value in acceleration {
      total += value.abs() * dt
    }
    total
  }
}

///|
/// Compute Arias intensity using SI acceleration samples.
pub fn arias_intensity(
  acceleration : Array[Double],
  sample_rate : Double,
) -> Double {
  if acceleration.is_empty() || sample_rate <= 0.0 {
    return 0.0
  }
  let gravity = 9.80665
  let dt = 1.0 / sample_rate
  let mut total = 0.0
  for value in acceleration {
    total += value * value * dt
  }
  @math.PI / (2.0 * gravity) * total
}

///|
/// Return a cumulative Arias intensity curve.
pub fn cumulative_arias_intensity(
  acceleration : Array[Double],
  sample_rate : Double,
) -> Array[Double] {
  let result = Array::make(acceleration.length(), 0.0)
  if sample_rate <= 0.0 {
    return result
  }
  let gravity = 9.80665
  let coefficient = @math.PI / (2.0 * gravity * sample_rate)
  let mut total = 0.0
  for index = 0; index < acceleration.length(); index = index + 1 {
    total += acceleration[index] * acceleration[index] * coefficient
    result[index] = total
  }
  result
}

///|
/// Interpolate the first index where a cumulative curve reaches a fraction.
pub fn cumulative_fraction_index(
  curve : Array[Double],
  fraction : Double,
) -> Int {
  if curve.is_empty() {
    return -1
  }
  let total = curve[curve.length() - 1]
  if total <= 0.0 {
    return 0
  }
  let target = total * fraction.max(0.0).min(1.0)
  for index = 0; index < curve.length(); index = index + 1 {
    if curve[index] >= target {
      return index
    }
  }
  curve.length() - 1
}

///|
/// Compute bracketed and significant Arias duration measures.
pub fn arias_duration_measures(
  acceleration : Array[Double],
  sample_rate : Double,
  start_time : Double,
) -> DurationMeasures {
  let curve = cumulative_arias_intensity(acceleration, sample_rate)
  if curve.is_empty() || sample_rate <= 0.0 {
    return {
      start_time,
      end_time: start_time,
      bracketed_duration: 0.0,
      significant_start: start_time,
      significant_end: start_time,
      significant_duration: 0.0,
    }
  }
  let i5 = cumulative_fraction_index(curve, 0.05)
  let i95 = cumulative_fraction_index(curve, 0.95)
  let i5s = cumulative_fraction_index(curve, 0.05)
  let i75 = cumulative_fraction_index(curve, 0.75)
  let i95s = cumulative_fraction_index(curve, 0.95)
  let t5 = start_time + i5.to_double() / sample_rate
  let t95 = start_time + i95.to_double() / sample_rate
  let t5s = start_time + i5s.to_double() / sample_rate
  let t95s = start_time + i95s.to_double() / sample_rate
  ignore(i75)
  {
    start_time: t5,
    end_time: t95,
    bracketed_duration: (t95 - t5).max(0.0),
    significant_start: t5s,
    significant_end: t95s,
    significant_duration: (t95s - t5s).max(0.0),
  }
}

///|
/// Compute peak acceleration, velocity, displacement, and duration measures.
pub fn ground_motion_metrics(
  acceleration : Array[Double],
  sample_rate : Double,
  start_time : Double,
) -> GroundMotionMetrics {
  let velocity = integrate_acceleration(acceleration, sample_rate, true)
  let displacement = integrate_velocity(velocity, sample_rate, true)
  let duration = arias_duration_measures(acceleration, sample_rate, start_time)
  {
    pga: absolute_peak(acceleration),
    pgv: absolute_peak(velocity),
    pgd: absolute_peak(displacement),
    rms_acceleration: rms_samples(acceleration),
    arias_intensity: arias_intensity(acceleration, sample_rate),
    cav: cumulative_absolute_velocity(acceleration, sample_rate),
    significant_duration: duration.significant_duration,
    bracketed_duration: duration.bracketed_duration,
    dominant_period: 0.0,
  }
}

///|
/// Compute a damped single-degree-of-freedom response ordinate.
pub fn response_ordinate(
  acceleration : Array[Double],
  sample_rate : Double,
  period_seconds : Double,
  damping : Double,
) -> ResponseSpectrumPoint {
  if acceleration.is_empty() || sample_rate <= 0.0 || period_seconds <= 0.0 {
    return {
      period_seconds,
      damping,
      spectral_acceleration: 0.0,
      spectral_velocity: 0.0,
      spectral_displacement: 0.0,
    }
  }
  let dt = 1.0 / sample_rate
  let omega = 2.0 * @math.PI / period_seconds
  let stiffness = omega * omega
  let damping_term = 2.0 * damping.max(0.0) * omega
  let mut displacement = 0.0
  let mut velocity = 0.0
  let mut previous_displacement = 0.0
  let mut peak_displacement = 0.0
  let mut peak_velocity = 0.0
  let mut peak_relative_acceleration = 0.0
  for input in acceleration {
    let restoring = stiffness * displacement
    let dissipative = damping_term * velocity
    let relative_acceleration = input - restoring - dissipative
    previous_displacement = displacement
    displacement += velocity * dt + 0.5 * relative_acceleration * dt * dt
    velocity += relative_acceleration * dt
    let actual_velocity = (displacement - previous_displacement) / dt
    peak_displacement = peak_displacement.max(displacement.abs())
    peak_velocity = peak_velocity.max(actual_velocity.abs())
    peak_relative_acceleration = peak_relative_acceleration.max(
      relative_acceleration.abs(),
    )
  }
  {
    period_seconds,
    damping,
    spectral_acceleration: peak_relative_acceleration.max(
      stiffness * peak_displacement,
    ),
    spectral_velocity: peak_velocity,
    spectral_displacement: peak_displacement,
  }
}

///|
/// Compute a response spectrum over a supplied period array.
pub fn response_spectrum(
  acceleration : Array[Double],
  sample_rate : Double,
  periods : Array[Double],
  damping : Double,
) -> ResponseSpectrum {
  let points = Array::new(capacity=periods.length())
  let mut peak_period = 0.0
  let mut peak_acceleration = 0.0
  for period in periods {
    let point = response_ordinate(acceleration, sample_rate, period, damping)
    points.push(point)
    if point.spectral_acceleration > peak_acceleration {
      peak_acceleration = point.spectral_acceleration
      peak_period = period
    }
  }
  { points, sample_rate, damping, peak_period, peak_acceleration }
}

///|
/// Return common logarithmically spaced engineering periods.
pub fn standard_response_periods(
  minimum : Double,
  maximum : Double,
  points : Int,
) -> Array[Double] {
  let result = Array::new()
  if minimum <= 0.0 || maximum < minimum || points <= 0 {
    return result
  }
  if points == 1 {
    result.push(minimum)
    return result
  }
  let log_min = @math.ln(minimum)
  let log_max = @math.ln(maximum)
  for index = 0; index < points; index = index + 1 {
    let fraction = index.to_double() / (points - 1).to_double()
    result.push(@math.exp(log_min + fraction * (log_max - log_min)))
  }
  result
}

///|
/// Estimate an oscillator period from a spectrum's dominant frequency.
pub fn dominant_period_from_spectrum(spectrum : Spectrum) -> Double {
  let frequency = dominant_frequency(spectrum)
  if frequency <= 0.0 {
    0.0
  } else {
    1.0 / frequency
  }
}

///|
/// Estimate an intensity label from PGA in g and hypocentral distance.
pub fn estimate_intensity(
  pga_g : Double,
  hypocentral_distance_km : Double,
) -> IntensityEstimate {
  let distance = hypocentral_distance_km.max(1.0)
  let attenuation = pga_g.abs() / (1.0 + distance / 100.0)
  let intensity = (3.0 * @math.log10(attenuation.max(1.0e-8)) + 7.0)
    .max(0.0)
    .min(12.0)
  let rounded = intensity.round().to_int()
  let label = if rounded >= 10 {
    "extreme"
  } else if rounded >= 7 {
    "strong"
  } else if rounded >= 4 {
    "moderate"
  } else {
    "light"
  }
  { pga_g, hypocentral_distance_km: distance, intensity, label }
}

///|
/// Convert SI acceleration to units of gravity.
pub fn acceleration_to_g(acceleration : Double) -> Double {
  acceleration / 9.80665
}

///|
/// Build strong-motion measures directly from a trace.
pub fn trace_ground_motion_metrics(trace : Trace) -> GroundMotionMetrics {
  ground_motion_metrics(trace.data, trace.sampling_rate, trace.start_time)
}

///|
/// Return a stable response-spectrum CSV table.
pub fn response_spectrum_to_csv(spectrum : ResponseSpectrum) -> String {
  let mut output = "period,damping,spectral_acceleration,spectral_velocity,spectral_displacement\n"
  for point in spectrum.points {
    output = output +
      point.period_seconds.to_string() +
      "," +
      point.damping.to_string() +
      "," +
      point.spectral_acceleration.to_string() +
      "," +
      point.spectral_velocity.to_string() +
      "," +
      point.spectral_displacement.to_string() +
      "\n"
  }
  output
}

///|
/// Format strong-motion metrics for an operator report.
pub fn ground_motion_metrics_text(metrics : GroundMotionMetrics) -> String {
  "pga=" +
  metrics.pga.to_string() +
  " pgv=" +
  metrics.pgv.to_string() +
  " pgd=" +
  metrics.pgd.to_string() +
  " rms=" +
  metrics.rms_acceleration.to_string() +
  " arias=" +
  metrics.arias_intensity.to_string() +
  " cav=" +
  metrics.cav.to_string() +
  " significant_duration=" +
  metrics.significant_duration.to_string() +
  " bracketed_duration=" +
  metrics.bracketed_duration.to_string() +
  " dominant_period=" +
  metrics.dominant_period.to_string()
}