geom-plane.stanza

defpackage geom/plane :
  import core
  import math
  import geom/vec
  import geom/line
  import geom/polyline

public defstruct Plane :
  center : Vec3
  normal : Vec3
with :
  printer => true

public defn parallel-vector (p:Plane) -> Vec3 :
  val vs = to-tuple $ for v in [x3(1.0), y3(1.0), z3(1.0)] filter : v != normal(p)
  normalize(vs[0] % normal(p))

public defn offset (p:Plane, o:Double) -> Plane :
  Plane(center(p) + o * normal(p), normal(p))

public defn invert (p:Plane) -> Plane :
  Plane(center(p), -1.0 * normal(p))

public defn close-loop (il:Seqable<Vec3>) -> Tuple<Vec3> :
  val l = to-tuple(il)
  to-tuple $ cat(l, [l[0]])

public defn to-polyline (p:Plane, r:Double) -> PolyLine3 :
  val v1 = parallel-vector(p)
  val v2 = normalize(normal(p) % v1)
  val l  = rect(center(p), v1, v2, r)
  val res = PolyLine3 $ [close-loop(l), [l[0] l[2]], [l[1] l[3]]]
  res

public defn distance (p:Plane) -> Double :
  dot(normal(p), center(p))

public defn distance (p:Plane, v:Vec3) -> Double :
  dot(normal(p), v) - distance(p)

public defn project (p:Plane, v:Vec3) -> Vec3 :
  v - distance(p, v) * normal(p)

public defn intersect (p:Plane, l:Line3) -> False|Vec3 :
  val d = dot(dir(l), normal(p))
  if d == 0.0 :
    false
  else :
    ; val t = (dot(normal(p), pos(l)) - distance(p)) / d
    val t = dot(center(p) - pos(l), normal(p)) / d
    eval(l, t)

public defn intersect (p0:Plane, p1:Plane) -> False|Line3 :
  if normal(p0) != normal(p1) :
    val n01 = normalize(normal(p0) % normal(p1))
    val d0 = normalize(n01 % normal(p0))
    val d1 = normalize(n01 % normal(p1))
    val l0 = Line3(center(p0), d0)
    val l1 = Line3(center(p1), d1)
    val c01 = closest-point-to(l0, l1)
    Line3(c01, n01)