(* rh takes an embedding and a facet and returns the pair (g, dep).
dep is a linear function of n space. (g +. (dep p)) is the distance
from p to the plane of the facet. *)
let rh embd (vl, (o, _)) =
   let va = Array.of_list vl in
   let n = Array.length va in
   let norg = Linear.vneg embd.(va.(n-1)) in
   let nm = Array.init (n-1) (fun j -> Linear.vadd norg embd.(va.(j))) in
   let un = let z = cross nm in if o then z else Linear.vneg z in
   let nx = Linear.ip un in
   (nx norg, nx);;

(* Does a ray pierce a facet? *)
type situ = Ingress of float | Egress of float | Miss;;
let rpf embd fac ((pos, dir ) : ray) = let (g, f) = rh embd fac and n = Array.length embd.(0) in 
let side = (f pos) and app = (f dir) in let off = side -. g  in let time = off /. app in
if let rec inside sc = sc = 0 || (Linear.determ
  (Array.init n (fun j -> if j < n-2 then Linear.vadd () )) > 0. && inside (sc - 1)) in inside n
then (if off < 0. then Ingress time else Egress time) else Miss;;