=== modified file 'src/point.ml'
--- src/point.ml	2011-09-06 15:20:02 +0000
+++ src/point.ml	2012-10-12 01:15:32 +0000
@@ -27,7 +27,7 @@
   | p -> p
 
 let to_spherical = function
-  | Cartesian (x, y, z) -> Spherical (acos z, atan (y /. x))
+  | Cartesian (x, y, z) -> Spherical (acos z, atan2 y x)
   | p -> p
 
 let rec ( =* ) p p' = match (p, p') with
@@ -55,3 +55,54 @@
   | _ -> assert false) (Array.map to_cartesian points);
   Cartesian (!sx *. factor, !sy *. factor, !sz *. factor)
 
+let rec plus p p' = match p, p' with
+  | Cartesian (x, y, z), Cartesian (x', y', z') ->
+    Cartesian (x +. x', y +. y', z +. z')
+  | _ -> plus (to_cartesian p) (to_cartesian p')
+
+let rec times s = function
+  | Cartesian (x, y, z) -> Cartesian (s *. x, s *. y, s *. z)
+  | (Spherical _) as p -> times s (to_cartesian p)
+
+let rec cross p p' = match p, p' with
+  | Cartesian (x, y, z), Cartesian (x', y', z') ->
+    Cartesian (y *. z' -. z *. y', z *. x' -. x *. z', x *. y' -. y *. x')
+  | _ -> cross (to_cartesian p) (to_cartesian p')
+
+let magnitude = function
+  | Cartesian (x, y, z) -> sqrt (x *. x +. y *. y +. z *. z)
+  | Spherical _ -> 1.0
+
+let norm = function
+  | (Cartesian _) as p -> times (1.0 /. magnitude p) p
+  | (Spherical _) as p -> p
+
+let midpoint p p' =
+  (* Midpoint of arc from p to p'
+   *)
+  norm (times 0.5 (plus p p'))
+
+let incenter p1 p2 p3 =
+  (* Calculate spherical incenter by intersecting angle bisectors.
+   *
+   * a12 = aux point at fixed dist along p1->p2
+   * a13 = aux point at fixed dist along p1->p3
+   * m1  = midpoint of a12->a13
+   * p1->m1 bisects angle at p1
+   *
+   * Similarly a21, a23, m2, p2->m2 are used to bisect angle at p2.
+   *
+   * Intersect great circles of p1->m1 and p2->m2 to get antipodal
+   * incenter candidates.  Pick the one that's closest to the centroid
+   * of p1, p2, p3.
+   *)
+  let a12 = norm (plus p1 (norm (cross (cross p1 p2) p1))) in
+  let a13 = norm (plus p1 (norm (cross (cross p1 p3) p1))) in
+  let m1 = midpoint a12 a13 in
+  let a21 = norm (plus p2 (norm (cross (cross p2 p1) p2))) in
+  let a23 = norm (plus p2 (norm (cross (cross p2 p3) p2))) in
+  let m2 = midpoint a21 a23 in
+  let cand1 = norm (cross (cross p1 m1) (cross p2 m2)) in
+  let cand2 = times (-1.0) cand1 in
+  let dc = distance (center [| p1; p2; p3 |]) in
+  if dc cand1 < dc cand2 then cand1 else cand2

=== modified file 'src/point.mli'
--- src/point.mli	2011-09-06 15:20:02 +0000
+++ src/point.mli	2012-10-11 20:57:03 +0000
@@ -29,3 +29,5 @@
 val to_string : ?after:string -> point -> string
 
 val center : point array -> point
+
+val incenter : point -> point -> point -> point

=== modified file 'src/wythoff.ml'
--- src/wythoff.ml	2011-09-14 10:20:57 +0000
+++ src/wythoff.ml	2012-10-11 23:50:11 +0000
@@ -76,7 +76,7 @@
   let x = match barpos with
     | One -> point_p
     | Bissected -> P.Spherical (p0q /. ( 1. +. q0r /. p0r), 0.)
-    | Incenter -> failwith "Not yet implemented"
+    | Incenter -> P.incenter point_p point_q point_r
     | Snub -> failwith "Not yet implemented"
     | Custom x -> x
   in