///|
pub(all) struct BoundaryIntersection {
  line_fraction : Double
  boundary_fraction : Double
  edge_index : Int
  composition : @domain.TernaryComposition
} derive(Debug, Eq, ToJson)

///|
pub(all) struct BoundaryProjection {
  edge_index : Int
  edge_fraction : Double
  composition : @domain.TernaryComposition
  distance : Double
} derive(Debug, Eq, ToJson)

///|
/// Closed polygonal approximation of a ternary liquid-liquid binodal region.
///
/// Both branches are supplied from low to high solute fraction. The raffinate
/// branch is followed forward and the extract branch backward to form a closed
/// two-phase envelope.
pub(all) struct BinodalCurve {
  raffinate_branch : Array[@domain.TernaryComposition]
  extract_branch : Array[@domain.TernaryComposition]
  polygon : Array[@domain.TernaryComposition]
} derive(Debug, Eq)

///|
fn validate_branch(
  branch : ArrayView[@domain.TernaryComposition],
  name : String,
) -> Unit raise @domain.ExtractionError {
  if branch.length() < 2 {
    raise @domain.ExtractionError::InsufficientData(
      required=2,
      actual=branch.length(),
    )
  }
  for index in 1.. BinodalCurve raise @domain.ExtractionError {
  validate_branch(raffinate_branch, "raffinate branch")
  validate_branch(extract_branch, "extract branch")
  let polygon = raffinate_branch.to_owned()
  for reverse_index in 0.. Double {
  ax * by - ay * bx
}

///|
fn squared_distance(first : CartesianPoint, second : CartesianPoint) -> Double {
  let dx = first.x - second.x
  let dy = first.y - second.y
  dx * dx + dy * dy
}

///|
pub fn BinodalCurve::contains(
  self : BinodalCurve,
  composition : @domain.TernaryComposition,
  include_boundary? : Bool = true,
  tolerance? : Double = 1.0e-10,
) -> Bool {
  let point = to_cartesian(composition)
  let mut inside = false
  for index in 0..= -tolerance &&
      dot <= length_squared + tolerance {
      return include_boundary
    }
    let crosses_ray = (a.y > point.y) != (b.y > point.y)
    if crosses_ray {
      let intersection_x = a.x + (point.y - a.y) * (b.x - a.x) / (b.y - a.y)
      if intersection_x > point.x {
        inside = !inside
      }
    }
  }
  inside
}

///|
fn project_to_segment(
  point : CartesianPoint,
  a : CartesianPoint,
  b : CartesianPoint,
) -> (Double, CartesianPoint) {
  let abx = b.x - a.x
  let aby = b.y - a.y
  let length_squared = abx * abx + aby * aby
  let fraction = (if length_squared > 0.0 {
    ((point.x - a.x) * abx + (point.y - a.y) * aby) / length_squared
  } else {
    0.0
  }).clamp(min=0.0, max=1.0)
  (fraction, { x: a.x + fraction * abx, y: a.y + fraction * aby })
}

///|
pub fn BinodalCurve::nearest_boundary(
  self : BinodalCurve,
  composition : @domain.TernaryComposition,
) -> BoundaryProjection raise @domain.ExtractionError {
  let point = to_cartesian(composition)
  let mut best_edge = 0
  let mut best_fraction = 0.0
  let mut best_point = to_cartesian(self.polygon[0])
  let mut best_distance_squared = 1.0e300
  for index in 0.. (Double, Double)? {
  let rx = p2.x - p.x
  let ry = p2.y - p.y
  let sx = q2.x - q.x
  let sy = q2.y - q.y
  let denominator = cross_2d(rx, ry, sx, sy)
  if denominator.abs() <= tolerance {
    return None
  }
  let qpx = q.x - p.x
  let qpy = q.y - p.y
  let line_fraction = cross_2d(qpx, qpy, sx, sy) / denominator
  let boundary_fraction = cross_2d(qpx, qpy, rx, ry) / denominator
  if line_fraction >= -tolerance &&
    line_fraction <= 1.0 + tolerance &&
    boundary_fraction >= -tolerance &&
    boundary_fraction <= 1.0 + tolerance {
    Some(
      (
        line_fraction.clamp(min=0.0, max=1.0),
        boundary_fraction.clamp(min=0.0, max=1.0),
      ),
    )
  } else {
    None
  }
}

///|
/// Intersect a ternary mixing line with the polygonal two-phase boundary.
pub fn BinodalCurve::intersections(
  self : BinodalCurve,
  from : @domain.TernaryComposition,
  to : @domain.TernaryComposition,
  tolerance? : Double = 1.0e-10,
) -> Array[BoundaryIntersection] raise @domain.ExtractionError {
  @domain.require_positive("intersection tolerance", tolerance)
  let p = to_cartesian(from)
  let p2 = to_cartesian(to)
  if squared_distance(p, p2) <= tolerance * tolerance {
    raise @domain.ExtractionError::InvalidArgument(
      message="intersection line endpoints must be distinct",
    )
  }
  let results : Array[BoundaryIntersection] = []
  for index in 0.. {
        let duplicate = results.any(item => {
          (item.line_fraction - line_fraction).abs() <= tolerance
        })
        if !duplicate {
          results.push({
            line_fraction,
            boundary_fraction,
            edge_index: index,
            composition: from.mix(
              to,
              self_amount=1.0 - line_fraction,
              other_amount=line_fraction,
            ),
          })
        }
      }
      None => ()
    }
  }
  results.sort_by((left, right) => {
    left.line_fraction.compare(right.line_fraction)
  })
  results
}

///|
pub fn BinodalCurve::estimated_plait_point(
  self : BinodalCurve,
) -> @domain.TernaryComposition raise @domain.ExtractionError {
  let raffinate = self.raffinate_branch[self.raffinate_branch.length() - 1]
  let extract = self.extract_branch[self.extract_branch.length() - 1]
  raffinate.mix(extract, self_amount=1.0, other_amount=1.0)
}