///|
pub struct MaterialProperty {
  name : String
  molecular_weight : Float
  density : Float
  heat_capacity : Float
  boiling_point : Float
  safety_limit : Float
} derive(Debug, Eq)

///|
pub fn material_property(
  name : String,
  molecular_weight : Float,
  density : Float,
  heat_capacity : Float,
  boiling_point : Float,
  safety_limit : Float,
) -> MaterialProperty {
  {
    name,
    molecular_weight,
    density,
    heat_capacity,
    boiling_point,
    safety_limit,
  }
}

///|
pub fn MaterialProperty::is_valid(self : MaterialProperty) -> Bool {
  self.name.trim().length() > 0 &&
  self.molecular_weight > 0.0 &&
  self.density > 0.0 &&
  self.heat_capacity > 0.0
}

///|
pub fn MaterialProperty::thermal_capacity(
  self : MaterialProperty,
  mass : Float,
) -> Float {
  mass * self.heat_capacity
}

///|
pub struct ReactorCase {
  name : String
  feed : Float
  product : Float
  byproduct : Float
  residence_time : Float
  temperature : Float
  pressure : Float
} derive(Debug, Eq)

///|
pub fn reactor_case(
  name : String,
  feed : Float,
  product : Float,
  byproduct : Float,
  residence_time : Float,
  temperature : Float,
  pressure : Float,
) -> ReactorCase {
  { name, feed, product, byproduct, residence_time, temperature, pressure }
}

///|
pub fn ReactorCase::conversion(self : ReactorCase) -> Float {
  if self.feed <= 0.0 {
    0.0
  } else {
    clamp((self.feed - self.product) / self.feed, 0.0, 1.0)
  }
}

///|
pub fn ReactorCase::product_yield(self : ReactorCase) -> Float {
  if self.feed <= 0.0 {
    0.0
  } else {
    clamp(self.product / self.feed, 0.0, 1.0)
  }
}

///|
pub fn ReactorCase::selectivity(self : ReactorCase) -> Float {
  if self.byproduct <= 0.0 {
    self.product
  } else {
    self.product / self.byproduct
  }
}

///|
pub fn ReactorCase::space_time_yield(self : ReactorCase) -> Float {
  if self.residence_time <= 0.0 {
    0.0
  } else {
    self.product / self.residence_time
  }
}

///|
pub fn ReactorCase::is_safe(
  self : ReactorCase,
  temperature_limit : Float,
  pressure_limit : Float,
) -> Bool {
  self.temperature >= 0.0 &&
  self.temperature <= temperature_limit &&
  self.pressure >= 0.0 &&
  self.pressure <= pressure_limit
}

///|
pub struct DistillationCase {
  feed_rate : Float
  feed_fraction : Float
  distillate_rate : Float
  distillate_fraction : Float
  bottoms_rate : Float
  bottoms_fraction : Float
  reflux_ratio : Float
  stages : Int
} derive(Debug, Eq)

///|
pub fn distillation_case(
  feed_rate : Float,
  feed_fraction : Float,
  distillate_rate : Float,
  distillate_fraction : Float,
  bottoms_rate : Float,
  bottoms_fraction : Float,
  reflux_ratio : Float,
  stages : Int,
) -> DistillationCase {
  {
    feed_rate,
    feed_fraction,
    distillate_rate,
    distillate_fraction,
    bottoms_rate,
    bottoms_fraction,
    reflux_ratio,
    stages,
  }
}

///|
pub fn DistillationCase::component_balance(self : DistillationCase) -> Float {
  self.feed_rate * self.feed_fraction -
  self.distillate_rate * self.distillate_fraction -
  self.bottoms_rate * self.bottoms_fraction
}

///|
pub fn DistillationCase::total_balance(self : DistillationCase) -> Float {
  self.feed_rate - self.distillate_rate - self.bottoms_rate
}

///|
pub fn DistillationCase::is_closed(
  self : DistillationCase,
  tolerance? : Float = 0.001,
) -> Bool {
  self.component_balance().abs() <= tolerance &&
  self.total_balance().abs() <= tolerance
}

///|
pub fn DistillationCase::recovery(self : DistillationCase) -> Float {
  if self.feed_rate * self.feed_fraction == 0.0 {
    0.0
  } else {
    self.distillate_rate *
    self.distillate_fraction /
    (self.feed_rate * self.feed_fraction)
  }
}

///|
pub fn DistillationCase::minimum_stages(self : DistillationCase) -> Int {
  if self.reflux_ratio <= 0.0 {
    0
  } else {
    (Float::from_int(self.stages) *
    (1.0 + self.reflux_ratio) /
    (2.0 + self.reflux_ratio))
    .ceil()
    .to_int()
  }
}

///|
pub struct HeatExchangerCase {
  hot_inlet : Float
  hot_outlet : Float
  cold_inlet : Float
  cold_outlet : Float
  hot_capacity_rate : Float
  cold_capacity_rate : Float
  area : Float
  overall_u : Float
} derive(Debug, Eq)

///|
pub fn heat_exchanger_case(
  hot_inlet : Float,
  hot_outlet : Float,
  cold_inlet : Float,
  cold_outlet : Float,
  hot_capacity_rate : Float,
  cold_capacity_rate : Float,
  area : Float,
  overall_u : Float,
) -> HeatExchangerCase {
  {
    hot_inlet,
    hot_outlet,
    cold_inlet,
    cold_outlet,
    hot_capacity_rate,
    cold_capacity_rate,
    area,
    overall_u,
  }
}

///|
pub fn HeatExchangerCase::hot_duty(self : HeatExchangerCase) -> Float {
  self.hot_capacity_rate * (self.hot_inlet - self.hot_outlet)
}

///|
pub fn HeatExchangerCase::cold_duty(self : HeatExchangerCase) -> Float {
  self.cold_capacity_rate * (self.cold_outlet - self.cold_inlet)
}

///|
pub fn HeatExchangerCase::duty_residual(self : HeatExchangerCase) -> Float {
  self.hot_duty() - self.cold_duty()
}

///|
pub fn HeatExchangerCase::effectiveness(self : HeatExchangerCase) -> Float {
  let cold_capacity = if self.cold_capacity_rate < self.hot_capacity_rate {
    self.cold_capacity_rate
  } else {
    self.hot_capacity_rate
  }
  if cold_capacity <= 0.0 || self.hot_inlet <= self.cold_inlet {
    0.0
  } else {
    clamp(
      self.cold_duty() / (cold_capacity * (self.hot_inlet - self.cold_inlet)),
      0.0,
      1.0,
    )
  }
}

///|
pub fn HeatExchangerCase::estimated_area(self : HeatExchangerCase) -> Float {
  if self.overall_u <= 0.0 {
    0.0
  } else {
    self.hot_duty() / self.overall_u
  }
}

///|
pub fn HeatExchangerCase::area_margin(self : HeatExchangerCase) -> Float {
  self.area - self.estimated_area()
}

///|
pub struct PumpCase {
  flow_rate : Float
  differential_pressure : Float
  efficiency : Float
  fluid_density : Float
  installed_power : Float
} derive(Debug, Eq)

///|
pub fn pump_case(
  flow_rate : Float,
  differential_pressure : Float,
  efficiency : Float,
  fluid_density : Float,
  installed_power : Float,
) -> PumpCase {
  {
    flow_rate,
    differential_pressure,
    efficiency,
    fluid_density,
    installed_power,
  }
}

///|
pub fn PumpCase::hydraulic_power(self : PumpCase) -> Float {
  self.flow_rate * self.differential_pressure
}

///|
pub fn PumpCase::shaft_power(self : PumpCase) -> Float {
  if self.efficiency <= 0.0 {
    0.0
  } else {
    self.hydraulic_power() / self.efficiency
  }
}

///|
pub fn PumpCase::power_margin(self : PumpCase) -> Float {
  self.installed_power - self.shaft_power()
}

///|
pub fn PumpCase::specific_speed(self : PumpCase, head : Float) -> Float {
  if head <= 0.0 {
    0.0
  } else {
    self.flow_rate.sqrt() / (head * head.sqrt())
  }
}

///|
pub fn PumpCase::is_overloaded(self : PumpCase) -> Bool {
  self.power_margin() < 0.0
}

///|
pub struct PressureDropCase {
  length : Float
  diameter : Float
  velocity : Float
  density : Float
  viscosity : Float
  roughness : Float
  fittings_k : Float
} derive(Debug, Eq)

///|
pub fn pressure_drop_case(
  length : Float,
  diameter : Float,
  velocity : Float,
  density : Float,
  viscosity : Float,
  roughness : Float,
  fittings_k : Float,
) -> PressureDropCase {
  { length, diameter, velocity, density, viscosity, roughness, fittings_k }
}

///|
pub fn PressureDropCase::reynolds(self : PressureDropCase) -> Float {
  if self.viscosity <= 0.0 {
    0.0
  } else {
    self.density * self.velocity * self.diameter / self.viscosity
  }
}

///|
pub fn PressureDropCase::friction_factor(self : PressureDropCase) -> Float {
  let reynolds = self.reynolds()
  if reynolds <= 0.0 {
    0.0
  } else if reynolds < 2300.0 {
    64.0 / reynolds
  } else {
    let relative_roughness : Float = if self.diameter <= 0.0 {
      0.0
    } else {
      self.roughness / self.diameter
    }
    let denominator = relative_roughness / 3.7 + 5.74 / reynolds
    0.25 / (denominator * denominator)
  }
}

///|
pub fn PressureDropCase::major_loss(self : PressureDropCase) -> Float {
  if self.diameter <= 0.0 {
    0.0
  } else {
    self.friction_factor() *
    self.length /
    self.diameter *
    self.density *
    self.velocity *
    self.velocity /
    2.0
  }
}

///|
pub fn PressureDropCase::minor_loss(self : PressureDropCase) -> Float {
  self.fittings_k * self.density * self.velocity * self.velocity / 2.0
}

///|
pub fn PressureDropCase::total_loss(self : PressureDropCase) -> Float {
  self.major_loss() + self.minor_loss()
}

///|
pub struct SafetyReview {
  name : String
  measured : Float
  design_limit : Float
  warning_fraction : Float
  unit : String
} derive(Debug, Eq)

///|
pub fn safety_review(
  name : String,
  measured : Float,
  design_limit : Float,
  warning_fraction : Float,
  unit : String,
) -> SafetyReview {
  { name, measured, design_limit, warning_fraction, unit }
}

///|
pub fn SafetyReview::utilization(self : SafetyReview) -> Float {
  if self.design_limit <= 0.0 {
    0.0
  } else {
    self.measured / self.design_limit
  }
}

///|
pub fn SafetyReview::is_warning(self : SafetyReview) -> Bool {
  self.utilization() >= self.warning_fraction
}

///|
pub fn SafetyReview::is_exceeded(self : SafetyReview) -> Bool {
  self.utilization() > 1.0
}

///|
pub fn SafetyReview::status(self : SafetyReview) -> String {
  if self.is_exceeded() {
    "exceeded"
  } else if self.is_warning() {
    "warning"
  } else {
    "normal"
  }
}

///|
pub fn thermal_duty(
  mass_flow : Float,
  heat_capacity : Float,
  inlet : Float,
  outlet : Float,
) -> Float {
  mass_flow * heat_capacity * (outlet - inlet)
}

///|
pub fn latent_duty(
  mass_flow : Float,
  latent_heat : Float,
  vapor_fraction : Float,
) -> Float {
  mass_flow * latent_heat * vapor_fraction
}

///|
pub fn total_heat_duty(
  sensible : Float,
  latent : Float,
  reaction : Float,
) -> Float {
  sensible + latent + reaction
}

///|
pub fn log_mean_temperature_difference(
  delta_one : Float,
  delta_two : Float,
) -> Float {
  if delta_one <= 0.0 || delta_two <= 0.0 {
    0.0
  } else if (delta_one - delta_two).abs() < 0.0000001 {
    delta_one
  } else {
    2.0 * delta_one * delta_two / (delta_one + delta_two)
  }
}

///|
pub fn reynolds_number(
  density : Float,
  velocity : Float,
  diameter : Float,
  viscosity : Float,
) -> Float {
  if viscosity <= 0.0 {
    0.0
  } else {
    density * velocity * diameter / viscosity
  }
}

///|
pub fn froude_number(
  velocity : Float,
  length : Float,
  gravity : Float,
) -> Float {
  if length <= 0.0 || gravity <= 0.0 {
    0.0
  } else {
    velocity / (gravity * length).sqrt()
  }
}

///|
pub fn residence_time(volume : Float, flow_rate : Float) -> Float {
  if flow_rate <= 0.0 {
    0.0
  } else {
    volume / flow_rate
  }
}

///|
pub fn turnover_rate(volume : Float, flow_rate : Float) -> Float {
  if volume <= 0.0 {
    0.0
  } else {
    flow_rate / volume
  }
}

///|
pub fn recovery(inlet : Float, recovered : Float) -> Float {
  if inlet <= 0.0 {
    0.0
  } else {
    clamp(recovered / inlet, 0.0, 1.0)
  }
}

///|
pub fn loss_fraction(inlet : Float, recovered : Float) -> Float {
  1.0 - recovery(inlet, recovered)
}

///|
pub fn safety_factor(allowable : Float, applied : Float) -> Float {
  if applied <= 0.0 {
    0.0
  } else {
    allowable / applied
  }
}

///|
pub fn annualized_energy(power : Float, operating_hours : Float) -> Float {
  power * operating_hours
}

///|
pub fn carbon_intensity(energy : Float, emission_factor : Float) -> Float {
  energy * emission_factor
}

///|
pub fn normalized_error(observed : Float, expected : Float) -> Float {
  if expected == 0.0 {
    observed.abs()
  } else {
    (observed - expected).abs() / expected.abs()
  }
}

///|
pub fn percent_change(before : Float, after : Float) -> Float {
  if before == 0.0 {
    0.0
  } else {
    (after - before) / before
  }
}

///|
pub fn interpolate_linear(
  x : Float,
  x_one : Float,
  y_one : Float,
  x_two : Float,
  y_two : Float,
) -> Float {
  if x_two == x_one {
    y_one
  } else {
    y_one + (x - x_one) / (x_two - x_one) * (y_two - y_one)
  }
}

///|
pub fn weighted_pressure_drop(
  lengths : Array[Float],
  drops : Array[Float],
) -> Float {
  let mut total : Float = 0.0
  for index, length in lengths {
    if index < drops.length() {
      total = total + length * drops[index]
    }
  }
  total
}

///|
pub fn sum_positive(values : Array[Float]) -> Float {
  let mut total : Float = 0.0
  for value in values {
    if value > 0.0 {
      total = total + value
    }
  }
  total
}

///|
pub fn sum_negative(values : Array[Float]) -> Float {
  let mut total : Float = 0.0
  for value in values {
    if value < 0.0 {
      total = total + value
    }
  }
  total
}

///|
pub fn balance_residual(inlet : Float, outlet : Float) -> Float {
  inlet - outlet
}

///|
pub fn balance_error_percent(inlet : Float, outlet : Float) -> Float {
  if inlet == 0.0 {
    0.0
  } else {
    100.0 * balance_residual(inlet, outlet).abs() / inlet.abs()
  }
}