///|
pub(all) struct SensitivityPoint {
  parameter : String
  base_value : Double
  lower_value : Double
  upper_value : Double
  base_recovery : Double
  lower_recovery : Double
  upper_recovery : Double
  normalized_sensitivity : Double
} derive(Debug, Eq, ToJson)

///|
fn normalized_sensitivity(
  lower_input : Double,
  upper_input : Double,
  lower_output : Double,
  upper_output : Double,
  base_input : Double,
  base_output : Double,
) -> Double {
  let input_span = upper_input - lower_input
  let output_span = upper_output - lower_output
  if input_span == 0.0 || base_output == 0.0 {
    0.0
  } else {
    output_span / input_span * (base_input / base_output)
  }
}

///|
/// Central finite-difference sensitivity of single-stage recovery.
pub fn single_stage_sensitivity(
  carrier_flow~ : Double,
  feed_concentration~ : Double,
  solvent_flow~ : Double,
  distribution~ : Double,
  relative_step? : Double = 0.01,
) -> Array[SensitivityPoint] raise @domain.ExtractionError {
  @domain.require_positive("carrier flow", carrier_flow)
  @domain.require_positive("feed concentration", feed_concentration)
  @domain.require_positive("solvent flow", solvent_flow)
  @domain.require_positive("distribution coefficient", distribution)
  @domain.require_positive("relative step", relative_step)
  if relative_step >= 1.0 {
    raise @domain.ExtractionError::InvalidArgument(
      message="relative_step must be below one",
    )
  }
  let recovery = (flow, coefficient) => {
    @liquid.single_stage_constant_k(
      feed_flow=carrier_flow,
      feed_concentration~,
      solvent_flow=flow,
      distribution=coefficient,
    ).recovery
  }
  let base = recovery(solvent_flow, distribution)
  let flow_lower = solvent_flow * (1.0 - relative_step)
  let flow_upper = solvent_flow * (1.0 + relative_step)
  let flow_lower_recovery = recovery(flow_lower, distribution)
  let flow_upper_recovery = recovery(flow_upper, distribution)
  let distribution_lower = distribution * (1.0 - relative_step)
  let distribution_upper = distribution * (1.0 + relative_step)
  let distribution_lower_recovery = recovery(solvent_flow, distribution_lower)
  let distribution_upper_recovery = recovery(solvent_flow, distribution_upper)
  [
    {
      parameter: "solvent_flow",
      base_value: solvent_flow,
      lower_value: flow_lower,
      upper_value: flow_upper,
      base_recovery: base,
      lower_recovery: flow_lower_recovery,
      upper_recovery: flow_upper_recovery,
      normalized_sensitivity: normalized_sensitivity(
        flow_lower, flow_upper, flow_lower_recovery, flow_upper_recovery, solvent_flow,
        base,
      ),
    },
    {
      parameter: "distribution",
      base_value: distribution,
      lower_value: distribution_lower,
      upper_value: distribution_upper,
      base_recovery: base,
      lower_recovery: distribution_lower_recovery,
      upper_recovery: distribution_upper_recovery,
      normalized_sensitivity: normalized_sensitivity(
        distribution_lower, distribution_upper, distribution_lower_recovery, distribution_upper_recovery,
        distribution, base,
      ),
    },
  ]
}