///|
/// CAD_LENGTH_EPS is the linear coincidence tolerance for boundary checks.
const CAD_LENGTH_EPS : Double = 1.0e-9

///|
/// CAD_AREA_EPS is the area/perimeter degeneracy tolerance.
const CAD_AREA_EPS : Double = 1.0e-12

///|
/// CAD_ANGLE_EPS is the angular/radius agreement tolerance in radians.
const CAD_ANGLE_EPS : Double = 1.0e-6

///|
/// CAD_DEGENERATE_EPS is the threshold below which a length, area, or volume
/// is treated as zero.
const CAD_DEGENERATE_EPS : Double = 1.0e-12

///|
pub fn polygon_area_2d(pts : Array[Vec2]) -> Double {
  let n = pts.length()
  if n < 3 {
    return 0.0
  }
  let mut area = 0.0
  for i in 0.. Array[Vec2] {
  if polygon_area_2d(pts) < 0.0 {
    pts.rev()
  } else {
    pts
  }
}

///|
pub fn force_cw(pts : Array[Vec2]) -> Array[Vec2] {
  if polygon_area_2d(pts) > 0.0 {
    pts.rev()
  } else {
    pts
  }
}

///|
pub fn is_point_in_polygon(p : Vec2, poly : Array[Vec2]) -> Bool {
  let n = poly.length()
  if n < 3 {
    return false
  }
  let mut inside = false
  let (px, py) = p
  for i in 0.. py) != (yj > py) &&
      px < (xj - xi) * (py - yi) / (yj - yi) + xi
    if intersect {
      inside = !inside
    }
  }
  inside
}

///|
pub fn cycle_area_3d(pts : Array[Vec3], cycle : Array[Int32]) -> Vec3 {
  let n = cycle.length()
  if n < 3 {
    return (0.0, 0.0, 0.0)
  }
  let mut ax = 0.0
  let mut ay = 0.0
  let mut az = 0.0
  for i in 0.. Vec3 {
  let (ax, ay, az) = a
  let (bx, by, bz) = b
  (ax - bx, ay - by, az - bz)
}

///|
fn v_cross(a : Vec3, b : Vec3) -> Vec3 {
  let (ax, ay, az) = a
  let (bx, by, bz) = b
  (ay * bz - az * by, az * bx - ax * bz, ax * by - ay * bx)
}

///|
fn v_dot(a : Vec3, b : Vec3) -> Double {
  let (ax, ay, az) = a
  let (bx, by, bz) = b
  ax * bx + ay * by + az * bz
}

///|
fn v_len(a : Vec3) -> Double {
  let (x, y, z) = a
  Double::sqrt(x * x + y * y + z * z)
}

///|
fn v_norm(a : Vec3) -> Vec3 {
  let l = v_len(a)
  if l == 0.0 {
    (1.0, 0.0, 0.0)
  } else {
    let (x, y, z) = a
    (x / l, y / l, z / l)
  }
}

///|
fn reverse_cycle(verts : Array[Int32]) -> Array[Int32] {
  verts.rev()
}

///|
pub fn rotate_point(p : Vec3, rx : Double, ry : Double, rz : Double) -> Vec3 {
  let m = rotation_matrix(rx, ry, rz)
  apply_matrix(m, p)
}

///|
#derive(Debug, Eq)
pub(all) struct Matrix3x3 {
  m00 : Double
  m01 : Double
  m02 : Double
  m10 : Double
  m11 : Double
  m12 : Double
  m20 : Double
  m21 : Double
  m22 : Double
}

///|
pub extend Matrix3x3 with @debug.Debug::{to_repr}

///|
pub extend Matrix3x3 with Eq::{not_equal, equal}

///|
fn rotation_matrix(rx : Double, ry : Double, rz : Double) -> Matrix3x3 {
  let rad_x = rx * @math.PI / 180.0
  let rad_y = ry * @math.PI / 180.0
  let rad_z = rz * @math.PI / 180.0
  let sx = @math.sin(rad_x)
  let cx = @math.cos(rad_x)
  let sy = @math.sin(rad_y)
  let cy = @math.cos(rad_y)
  let sz = @math.sin(rad_z)
  let cz = @math.cos(rad_z)

  // Combined matrix Rz * Ry * Rx
  Matrix3x3::{
    m00: cy * cz,
    m01: sx * sy * cz - cx * sz,
    m02: cx * sy * cz + sx * sz,
    m10: cy * sz,
    m11: sx * sy * sz + cx * cz,
    m12: cx * sy * sz - sx * cz,
    m20: -sy,
    m21: sx * cy,
    m22: cx * cy,
  }
}

///|
fn apply_matrix(m : Matrix3x3, p : Vec3) -> Vec3 {
  let (x, y, z) = p
  (
    m.m00 * x + m.m01 * y + m.m02 * z,
    m.m10 * x + m.m11 * y + m.m12 * z,
    m.m20 * x + m.m21 * y + m.m22 * z,
  )
}

///|
fn multiply_matrices(a : Matrix3x3, b : Matrix3x3) -> Matrix3x3 {
  Matrix3x3::{
    m00: a.m00 * b.m00 + a.m01 * b.m10 + a.m02 * b.m20,
    m01: a.m00 * b.m01 + a.m01 * b.m11 + a.m02 * b.m21,
    m02: a.m00 * b.m02 + a.m01 * b.m12 + a.m02 * b.m22,
    m10: a.m10 * b.m00 + a.m11 * b.m10 + a.m12 * b.m20,
    m11: a.m10 * b.m01 + a.m11 * b.m11 + a.m12 * b.m21,
    m12: a.m10 * b.m02 + a.m11 * b.m12 + a.m12 * b.m22,
    m20: a.m20 * b.m00 + a.m21 * b.m10 + a.m22 * b.m20,
    m21: a.m20 * b.m01 + a.m21 * b.m11 + a.m22 * b.m21,
    m22: a.m20 * b.m02 + a.m21 * b.m12 + a.m22 * b.m22,
  }
}

///|
pub fn matrix_to_euler(m : Matrix3x3) -> Vec3 {
  // Extract Euler angles from rotation matrix assuming order Rz * Ry * Rx
  let ry = -@math.asin(m.m20)
  if Double::abs(@math.cos(ry)) > 0.0001 {
    let rx = @math.atan2(m.m21, m.m22)
    let rz = @math.atan2(m.m10, m.m00)
    (rx * 180.0 / @math.PI, ry * 180.0 / @math.PI, rz * 180.0 / @math.PI)
  } else {
    // Gimbal lock
    let rx = 0.0
    let rz = @math.atan2(-m.m01, m.m11)
    (rx, ry * 180.0 / @math.PI, rz * 180.0 / @math.PI)
  }
}

///|
/// add_plane_face_with_edge_loops emits a planar face whose outer boundary is a
/// polygon cycle and whose holes are already-built edge loops (for example
/// analytic circles). The outer cycle is oriented so its area vector points
/// along `expected_out`; the hole loops are used exactly as supplied.
fn add_plane_face_with_edge_loops(
  ctx : StepContext,
  edges : Array[EdgeInfo],
  pts : Array[Vec3],
  cp_refs : Array[@step.Ref],
  outer_cycle_in : Array[Int32],
  hole_loops_in : Array[@step.Ref],
  expected_out : Vec3,
  reverse~ : Bool = false,
) -> @step.Ref {
  let area_vec = cycle_area_3d(pts, outer_cycle_in)
  let cycle = if v_dot(area_vec, expected_out) < 0.0 {
    reverse_cycle(outer_cycle_in)
  } else {
    outer_cycle_in
  }
  let b0 = cycle[0]
  let axis_raw = v_norm(expected_out)
  let refd_raw = v_norm(v_sub(pts[cycle[1]], pts[cycle[0]]))
  let a2p_ref = ctx.axis2_placement_3d(
    cp_refs[b0],
    ctx.direction(axis_raw.0, axis_raw.1, axis_raw.2),
    ctx.direction(refd_raw.0, refd_raw.1, refd_raw.2),
  )
  let plane_ref = ctx.add(
    @repository.any_entity_plane(@geometry.Plane::new("", a2p_ref)),
  )
  let outer_loop_ref = face_loop(ctx, edges, cycle)
  let bounds : Array[@step.Ref] = [ctx.face_outer_bound(outer_loop_ref, true)]
  for loop_ref in hole_loops_in {
    bounds.push(ctx.face_bound(loop_ref, true))
  }
  ctx.advanced_face(bounds, plane_ref, !reverse)
}

///|
/// add_plane_face_from_loops emits a planar face from an already-built outer
/// edge loop (for example an analytic circle) and zero or more hole loops. The
/// caller supplies the loop orientations; this helper only attaches them to the
/// plane. The plane normal is `axis` when `reverse` is false.
fn add_plane_face_from_loops(
  ctx : StepContext,
  origin_cp : @step.Ref,
  axis : Vec3,
  refd : Vec3,
  outer_loop : @step.Ref,
  hole_loops : Array[@step.Ref],
  reverse~ : Bool = false,
) -> @step.Ref {
  let axis_raw = v_norm(axis)
  let ref_raw = v_norm(refd)
  let a2p_ref = ctx.axis2_placement_3d(
    origin_cp,
    ctx.direction(axis_raw.0, axis_raw.1, axis_raw.2),
    ctx.direction(ref_raw.0, ref_raw.1, ref_raw.2),
  )
  let plane_ref = ctx.add(
    @repository.any_entity_plane(@geometry.Plane::new("", a2p_ref)),
  )
  let bounds : Array[@step.Ref] = [ctx.face_outer_bound(outer_loop, true)]
  for loop_ref in hole_loops {
    bounds.push(ctx.face_bound(loop_ref, true))
  }
  ctx.advanced_face(bounds, plane_ref, !reverse)
}

///|
pub(all) struct EdgeInfo {
  start : Int32
  end : Int32
  edge : @step.Ref
}

///|
fn find_edge(
  edges : Array[EdgeInfo],
  a : Int32,
  b : Int32,
) -> (@step.Ref, Bool) {
  for e in edges {
    if e.start == a && e.end == b {
      return (e.edge, true)
    }
    if e.start == b && e.end == a {
      return (e.edge, false)
    }
  }
  abort("edge not found between points \{a} and \{b}")
}

///|
fn face_loop(
  ctx : StepContext,
  edges : Array[EdgeInfo],
  verts : Array[Int32],
) -> @step.Ref {
  let oe_refs : Array[@step.Ref] = []
  let n = verts.length()
  for i in 0.. @step.Ref {
  let a2p_ref = ctx.axis2_placement_3d(location_pt, axis_dir, ref_dir)
  let plane = @geometry.Plane::new("", a2p_ref)
  let plane_ref = ctx.add(@repository.any_entity_plane(plane))
  let loop_ref = face_loop(ctx, edges, verts_cycle)
  let fob_ref = ctx.face_outer_bound(loop_ref, true)
  ctx.advanced_face([fob_ref], plane_ref, !reverse)
}

///|
fn add_plane_face_oriented(
  ctx : StepContext,
  edges : Array[EdgeInfo],
  pts : Array[Vec3],
  cp_refs : Array[@step.Ref],
  verts_cycle : Array[Int32],
  expected_out : Vec3,
  reverse~ : Bool = false,
) -> @step.Ref {
  let area_vec = cycle_area_3d(pts, verts_cycle)
  let cycle = if v_dot(area_vec, expected_out) < 0.0 {
    reverse_cycle(verts_cycle)
  } else {
    verts_cycle
  }
  let b0 = cycle[0]
  let axis_raw = v_norm(expected_out)
  let refd_raw = v_norm(v_sub(pts[cycle[1]], pts[cycle[0]]))
  let axis_ref = ctx.direction(axis_raw.0, axis_raw.1, axis_raw.2)
  let ref_ref = ctx.direction(refd_raw.0, refd_raw.1, refd_raw.2)
  add_plane_face(ctx, edges, cycle, cp_refs[b0], axis_ref, ref_ref, reverse~)
}

///|
fn add_plane_face_with_holes_oriented(
  ctx : StepContext,
  edges : Array[EdgeInfo],
  pts : Array[Vec3],
  cp_refs : Array[@step.Ref],
  outer_cycle_in : Array[Int32],
  hole_cycles_in : Array[Array[Int32]],
  expected_out : Vec3,
  reverse~ : Bool = false,
) -> @step.Ref {
  let area_vec = cycle_area_3d(pts, outer_cycle_in)
  let outer_cycle = if v_dot(area_vec, expected_out) < 0.0 {
    reverse_cycle(outer_cycle_in)
  } else {
    outer_cycle_in
  }
  // Use the oriented outer loop's area vector. After the reversal above it
  // always points along expected_out; using the original area_vec instead
  // inverted every hole whenever the outer loop had to be reversed, so the
  // cap traversed its hole edges in the same direction as the adjacent
  // (twisted) side faces. That leaves the shell with inconsistently oriented
  // edges, which some STEP tessellators report as mesh gaps.
  let outer_area = v_dot(cycle_area_3d(pts, outer_cycle), expected_out)

  // Holes must be oriented opposite to the outer boundary for STEP faces.
  let hole_cycles = hole_cycles_in.map(fn(hc) {
    let ha_vec = cycle_area_3d(pts, hc)
    // If hole has same orientation as outer (relative to expected_out), reverse it.
    if v_dot(ha_vec, expected_out) * outer_area > 0.0 {
      reverse_cycle(hc)
    } else {
      hc
    }
  })

  //
  let b0 = outer_cycle[0]
  let axis_raw = v_norm(expected_out)
  let refd_raw = v_norm(v_sub(pts[outer_cycle[1]], pts[outer_cycle[0]]))
  let axis_ref = ctx.direction(axis_raw.0, axis_raw.1, axis_raw.2)
  let ref_ref = ctx.direction(refd_raw.0, refd_raw.1, refd_raw.2)

  //
  let a2p_ref = ctx.axis2_placement_3d(cp_refs[b0], axis_ref, ref_ref)
  let plane = @geometry.Plane::new("", a2p_ref)
  let plane_ref = ctx.add(@repository.any_entity_plane(plane))

  //
  let outer_loop_ref = face_loop(ctx, edges, outer_cycle)
  let fob_ref = ctx.face_outer_bound(outer_loop_ref, true)

  //
  let bounds : Array[@step.Ref] = [fob_ref]
  bounds.push_iter(
    hole_cycles
    .iter()
    .map(fn(loop_pts) {
      let loop_ref = face_loop(ctx, edges, loop_pts)
      ctx.face_bound(loop_ref, true)
    }),
  )
  ctx.advanced_face(bounds, plane_ref, !reverse)
}