///|
/// A sampled instrument response curve for calibration diagnostics.
pub struct ResponseCurvePoint {
  frequency : Double
  gain : Double
  normalized_gain : Double
  relative_db : Double
}

///|
/// A frequency response curve with metadata.
pub struct ResponseCurve {
  points : Array[ResponseCurvePoint]
  minimum_gain : Double
  maximum_gain : Double
  corner_frequency : Double
  bandwidth_hz : Double
}

///|
/// Diagnostics from applying or removing a response.
pub struct CalibrationDiagnostics {
  input_quantity : WaveQuantity
  output_quantity : WaveQuantity
  sensitivity : Double
  frequency : Double
  gain : Double
  input_peak : Double
  output_peak : Double
  scale_factor : Double
  stable : Bool
}

///|
/// Return a readable quantity name.
pub fn wave_quantity_name(quantity : WaveQuantity) -> String {
  match quantity {
    Counts => "counts"
    Displacement => "displacement"
    Velocity => "velocity"
    Acceleration => "acceleration"
  }
}

///|
/// Return a logarithmically spaced frequency axis.
pub fn response_frequency_axis(
  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 lower = @math.ln(minimum)
  let upper = @math.ln(maximum)
  for index = 0; index < points; index = index + 1 {
    let fraction = index.to_double() / (points - 1).to_double()
    result.push(@math.exp(lower + fraction * (upper - lower)))
  }
  result
}

///|
/// Compute a sampled response curve and normalize it at the reference frequency.
pub fn sample_response_curve(
  response : InstrumentResponse,
  frequencies : Array[Double],
) -> ResponseCurve {
  let points = Array::new(capacity=frequencies.length())
  let reference_gain = response_gain(response, response.reference_frequency)
    .abs()
    .max(1.0e-20)
  let mut minimum_gain = 1.0e300
  let mut maximum_gain = 0.0
  for frequency in frequencies {
    let gain = response_gain(response, frequency).abs()
    let normalized = gain / reference_gain
    let relative_db = 20.0 * @math.log10(normalized.max(1.0e-20))
    minimum_gain = minimum_gain.min(gain)
    maximum_gain = maximum_gain.max(gain)
    points.push({ frequency, gain, normalized_gain: normalized, relative_db })
  }
  let mut corner = response.reference_frequency
  for value in frequencies {
    if response_gain(response, value).abs() >= reference_gain / 2.0 {
      corner = value
      break
    }
  }
  {
    points,
    minimum_gain: if points.is_empty() {
      0.0
    } else {
      minimum_gain
    },
    maximum_gain,
    corner_frequency: corner,
    bandwidth_hz: if frequencies.is_empty() {
      0.0
    } else {
      maximum_sample(frequencies) - minimum_sample(frequencies)
    },
  }
}

///|
/// Interpolate a sampled response curve in frequency.
pub fn interpolate_response_curve(
  curve : ResponseCurve,
  frequency : Double,
) -> Double {
  if curve.points.is_empty() {
    return 0.0
  }
  if frequency <= curve.points[0].frequency {
    return curve.points[0].gain
  }
  for index = 1; index < curve.points.length(); index = index + 1 {
    let left = curve.points[index - 1]
    let right = curve.points[index]
    if frequency <= right.frequency {
      let fraction = (frequency - left.frequency) /
        (right.frequency - left.frequency).max(1.0e-12)
      return left.gain + fraction * (right.gain - left.gain)
    }
  }
  curve.points[curve.points.length() - 1].gain
}

///|
/// Apply a frequency-independent response scale to a waveform.
pub fn apply_response_scale(
  trace : Trace,
  response : InstrumentResponse,
  frequency : Double,
) -> Trace {
  let gain = response_gain(response, frequency)
  if gain == 0.0 {
    return trace.clone()
  }
  trace_with_data(trace, trace.data.map(value => value * gain))
}

///|
/// Remove a frequency-independent response scale safely.
pub fn remove_response_scale(
  trace : Trace,
  response : InstrumentResponse,
  frequency : Double,
) -> Trace {
  let gain = response_gain(response, frequency)
  if gain.abs() <= 1.0e-20 {
    return trace.clone()
  }
  trace_with_data(trace, trace.data.map(value => value / gain))
}

///|
/// Return diagnostics for removing an instrument response.
pub fn response_calibration_diagnostics(
  trace : Trace,
  response : InstrumentResponse,
  frequency : Double,
) -> CalibrationDiagnostics {
  let output = remove_response_scale(trace, response, frequency)
  let gain = response_gain(response, frequency)
  {
    input_quantity: response.input_quantity,
    output_quantity: response.output_quantity,
    sensitivity: response.sensitivity,
    frequency,
    gain,
    input_peak: absolute_peak(trace.data),
    output_peak: absolute_peak(output.data),
    scale_factor: if gain.abs() <= 1.0e-20 {
      0.0
    } else {
      1.0 / gain
    },
    stable: gain.abs() > 1.0e-20 && trace_is_well_formed(output),
  }
}

///|
/// Calculate a response-correction scale over a frequency band.
pub fn response_band_scale(
  response : InstrumentResponse,
  low_frequency : Double,
  high_frequency : Double,
  samples : Int,
) -> Double {
  let frequencies = response_frequency_axis(
    low_frequency, high_frequency, samples,
  )
  if frequencies.is_empty() {
    1.0
  } else {
    let gains = frequencies.map(frequency => {
      response_gain(response, frequency).abs()
    })
    median_samples(gains)
  }
}

///|
/// Return whether a response is safe to invert over a band.
pub fn response_is_invertible(
  response : InstrumentResponse,
  low_frequency : Double,
  high_frequency : Double,
  samples : Int,
  minimum_gain : Double,
) -> Bool {
  response_frequency_axis(low_frequency, high_frequency, samples).all(frequency => {
    response_gain(response, frequency).abs() >= minimum_gain.abs()
  })
}

///|
/// Compute a normalized response curve text line.
pub fn response_curve_summary(curve : ResponseCurve) -> String {
  "points=" +
  curve.points.length().to_string() +
  " min_gain=" +
  curve.minimum_gain.to_string() +
  " max_gain=" +
  curve.maximum_gain.to_string() +
  " corner=" +
  curve.corner_frequency.to_string() +
  " bandwidth=" +
  curve.bandwidth_hz.to_string()
}

///|
/// Export response curve points for plotting.
pub fn response_curve_to_csv(curve : ResponseCurve) -> String {
  let mut output = "frequency,gain,normalized_gain,relative_db\n"
  for point in curve.points {
    output = output +
      point.frequency.to_string() +
      "," +
      point.gain.to_string() +
      "," +
      point.normalized_gain.to_string() +
      "," +
      point.relative_db.to_string() +
      "\n"
  }
  output
}

///|
/// Format response calibration diagnostics.
pub fn calibration_diagnostics_text(
  diagnostics : CalibrationDiagnostics,
) -> String {
  "input=" +
  wave_quantity_name(diagnostics.input_quantity) +
  " output=" +
  wave_quantity_name(diagnostics.output_quantity) +
  " frequency=" +
  diagnostics.frequency.to_string() +
  " gain=" +
  diagnostics.gain.to_string() +
  " input_peak=" +
  diagnostics.input_peak.to_string() +
  " output_peak=" +
  diagnostics.output_peak.to_string() +
  " stable=" +
  diagnostics.stable.to_string()
}