///|
/// Reaction-kinetics and reactor-screening calculations.
pub struct ArrheniusModel {
  pre_exponential : Float
  activation_energy_kj : Float
  reference_temperature_k : Float
}

///|
pub fn arrhenius(
  pre_exponential : Float,
  activation_energy_kj : Float,
  reference_temperature_k : Float,
) -> ArrheniusModel {
  { pre_exponential, activation_energy_kj, reference_temperature_k }
}

///|
pub fn ArrheniusModel::rate_constant(
  self : ArrheniusModel,
  temperature_k : Float,
) -> Float {
  if temperature_k <= 0.0 {
    0.0
  } else {
    let ratio : Float = self.activation_energy_kj / 8.314 / temperature_k
    let factor : Float = if ratio >= 1.0 { 0.0 } else { 1.0 - ratio }
    self.pre_exponential * factor
  }
}

///|
pub fn ArrheniusModel::relative_rate(
  self : ArrheniusModel,
  temperature_k : Float,
) -> Float {
  if self.reference_temperature_k <= 0.0 {
    0.0
  } else {
    self.rate_constant(temperature_k) /
    self.rate_constant(self.reference_temperature_k)
  }
}

///|
pub struct ReactionOrder {
  order_a : Float
  order_b : Float
  stoich_a : Float
  stoich_b : Float
}

///|
pub fn reaction_order(
  order_a : Float,
  order_b : Float,
  stoich_a : Float,
  stoich_b : Float,
) -> ReactionOrder {
  { order_a, order_b, stoich_a, stoich_b }
}

///|
pub fn ReactionOrder::rate(
  self : ReactionOrder,
  constant : Float,
  concentration_a : Float,
  concentration_b : Float,
) -> Float {
  let a : Float = if concentration_a < 0.0 { 0.0 } else { concentration_a }
  let b : Float = if concentration_b < 0.0 { 0.0 } else { concentration_b }
  let order : Float = self.order_a + self.order_b
  let factor : Float = if order < 0.0 { 0.0 } else { order }
  constant * a * b * factor
}

///|
pub fn ReactionOrder::valid(self : ReactionOrder) -> Bool {
  self.order_a >= 0.0 &&
  self.order_b >= 0.0 &&
  self.stoich_a > 0.0 &&
  self.stoich_b >= 0.0
}

///|
pub struct CstrDesign {
  volume_m3 : Float
  flow_m3_h : Float
  rate_mol_m3_h : Float
  feed_mol_m3 : Float
}

///|
pub fn cstr_design(
  volume_m3 : Float,
  flow_m3_h : Float,
  rate_mol_m3_h : Float,
  feed_mol_m3 : Float,
) -> CstrDesign {
  { volume_m3, flow_m3_h, rate_mol_m3_h, feed_mol_m3 }
}

///|
pub fn CstrDesign::residence_hours(self : CstrDesign) -> Float {
  if self.flow_m3_h == 0.0 {
    0.0
  } else {
    self.volume_m3 / self.flow_m3_h
  }
}

///|
pub fn CstrDesign::conversion(self : CstrDesign) -> Float {
  if self.flow_m3_h == 0.0 || self.feed_mol_m3 == 0.0 {
    0.0
  } else {
    (self.rate_mol_m3_h * self.volume_m3 / self.flow_m3_h / self.feed_mol_m3)
    .min(1.0)
    .max(0.0)
  }
}

///|
pub fn CstrDesign::outlet(self : CstrDesign) -> Float {
  self.feed_mol_m3 * (1.0 - self.conversion())
}

///|
pub fn CstrDesign::valid(self : CstrDesign) -> Bool {
  self.volume_m3 > 0.0 &&
  self.flow_m3_h > 0.0 &&
  self.rate_mol_m3_h >= 0.0 &&
  self.feed_mol_m3 >= 0.0
}

///|
pub struct PfrDesign {
  volume_m3 : Float
  flow_m3_h : Float
  inlet_mol_m3 : Float
  rate_mol_m3_h : Float
}

///|
pub fn pfr_design(
  volume_m3 : Float,
  flow_m3_h : Float,
  inlet_mol_m3 : Float,
  rate_mol_m3_h : Float,
) -> PfrDesign {
  { volume_m3, flow_m3_h, inlet_mol_m3, rate_mol_m3_h }
}

///|
pub fn PfrDesign::space_time(self : PfrDesign) -> Float {
  if self.flow_m3_h == 0.0 {
    0.0
  } else {
    self.volume_m3 / self.flow_m3_h
  }
}

///|
pub fn PfrDesign::conversion(self : PfrDesign) -> Float {
  if self.inlet_mol_m3 == 0.0 || self.flow_m3_h == 0.0 {
    0.0
  } else {
    (self.rate_mol_m3_h * self.space_time() / self.inlet_mol_m3)
    .min(1.0)
    .max(0.0)
  }
}

///|
pub fn PfrDesign::outlet(self : PfrDesign) -> Float {
  self.inlet_mol_m3 * (1.0 - self.conversion())
}

///|
pub struct ReactionPath {
  name : String
  conversion : Float
  selectivity : Float
  yield_fraction : Float
}

///|
pub fn reaction_path(
  name : String,
  conversion : Float,
  selectivity : Float,
  yield_fraction : Float,
) -> ReactionPath {
  { name, conversion, selectivity, yield_fraction }
}

///|
pub fn ReactionPath::product_fraction(self : ReactionPath) -> Float {
  self.conversion * self.selectivity * self.yield_fraction
}

///|
pub fn ReactionPath::valid(self : ReactionPath) -> Bool {
  self.conversion >= 0.0 &&
  self.conversion <= 1.0 &&
  self.selectivity >= 0.0 &&
  self.selectivity <= 1.0 &&
  self.yield_fraction >= 0.0 &&
  self.yield_fraction <= 1.0
}

///|
pub fn reaction_path_table(paths : Array[ReactionPath]) -> ReportTable {
  let rows : Array[Array[String]] = []
  for p in paths {
    rows.push([
      p.name,
      "\{p.conversion}",
      "\{p.selectivity}",
      "\{p.yield_fraction}",
      "\{p.product_fraction()}",
    ])
  }
  table(
    ["path", "conversion", "selectivity", "yield", "product fraction"],
    rows,
  )
}

///|
pub fn first_order_conversion(
  rate_constant : Float,
  residence_hours : Float,
) -> Float {
  (rate_constant * residence_hours / (1.0 + rate_constant * residence_hours))
  .min(1.0)
  .max(0.0)
}

///|
pub fn first_order_half_life(rate_constant : Float) -> Float {
  if rate_constant <= 0.0 {
    0.0
  } else {
    0.693147 / rate_constant
  }
}

///|
pub fn reaction_heat(duty_mol_h : Float, heat_kj_mol : Float) -> Float {
  duty_mol_h * heat_kj_mol
}

///|
pub fn catalyst_productivity(
  product_kg : Float,
  catalyst_kg : Float,
  hours : Float,
) -> Float {
  if catalyst_kg == 0.0 || hours == 0.0 {
    0.0
  } else {
    product_kg / catalyst_kg / hours
  }
}

///|
pub fn reactor_gate(
  conversion : Float,
  minimum : Float,
  temperature_c : Float,
  maximum_c : Float,
) -> Bool {
  conversion >= minimum && temperature_c <= maximum_c
}