diff --git a/arrays.scad b/arrays.scad
index 36841fe..835fd65 100644
--- a/arrays.scad
+++ b/arrays.scad
@@ -18,20 +18,6 @@
 
 // Section: List Query Operations
 
-// Function: is_simple_list()
-// Description:
-//   Returns true just when all elements of `list` are simple values.
-// Usage:
-//   is_simple_list(list)
-// Arguments:
-//   list = The list to check.
-// Example:
-//   a = is_simple_list([3,4,5,6,7,8,9]);  Returns: true
-//   b = is_simple_list([3,4,5,[6],7,8]);  Returns: false
-function is_simple_list(list) =
-		is_list(list)
-		&& []==[for(e=list) if(is_list(e)) 0];
-
 
 // Function: select()
 // Description:
@@ -73,9 +59,6 @@ function select(list, start, end=undef) =
             :   concat([for (i = [s:1:l-1]) list[i]], [for (i = [0:1:e]) list[i]]) ;
 
 
-
-
-
 // Function: slice()
 // Description:
 //   Returns a slice of a list.  The first item is index 0.
@@ -101,24 +84,23 @@ function slice(list,start,end) =
         ) [for (i=[s:1:e-1]) if (e>s) list[i]];
 
 
-
-
 // Function: in_list()
 // Description: Returns true if value `val` is in list `list`. When `val==NAN` the answer will be false for any list.
 // Arguments:
 //   val = The simple value to search for.
 //   list = The list to search.
-//   idx = If given, searches the given subindexes for matches for `val`.
+//   idx = If given, searches the given subindex for matches for `val`.
 // Example:
 //   in_list("bar", ["foo", "bar", "baz"]);  // Returns true.
 //   in_list("bee", ["foo", "bar", "baz"]);  // Returns false.
 //   in_list("bar", [[2,"foo"], [4,"bar"], [3,"baz"]], idx=1);  // Returns true.
 function in_list(val,list,idx=undef) = 
+    assert( is_list(list) && (is_undef(idx) || is_finite(idx)),
+            "Invalid input." )
     let( s = search([val], list, num_returns_per_match=1, index_col_num=idx)[0] )
-    s==[] || s[0]==[] ? false
+    s==[] || s==[[]] ? false
     : is_undef(idx) ? val==list[s] 
     : val==list[s][idx];
-    
 
 
 // Function: min_index()
@@ -209,7 +191,6 @@ function repeat(val, n, i=0) =
     [for (j=[1:1:n[i]]) repeat(val, n, i+1)];
 
 
-
 // Function: list_range()
 // Usage:
 //   list_range(n, [s], [e])
@@ -249,7 +230,6 @@ function list_range(n=undef, s=0, e=undef, step=undef) =
     
 
 
-
 // Section: List Manipulation
 
 // Function: reverse()
@@ -315,8 +295,6 @@ function deduplicate(list, closed=false, eps=EPSILON) =
     : [for (i=[0:1:l-1]) if (i==end || !approx(list[i], list[(i+1)%l], eps)) list[i]];
 
 
-
-
 // Function: deduplicate_indexed()
 // Usage:
 //   new_idxs = deduplicate_indexed(list, indices, [closed], [eps]);
@@ -351,8 +329,6 @@ function deduplicate_indexed(list, indices, closed=false, eps=EPSILON) =
     ];
 
 
-
-
 // Function: repeat_entries()
 // Usage:
 //   newlist = repeat_entries(list, N)
@@ -390,8 +366,6 @@ function repeat_entries(list, N, exact = true) =
                : [for (val=reps_guess) round(val)]
     )
     [for(i=[0:length-1]) each repeat(list[i],reps[i])];
-    
-
 
 
 // Function: list_set()
@@ -431,7 +405,6 @@ function list_set(list=[],indices,values,dflt=0,minlen=0) =
                 dflt ,
           each repeat(dflt, minlen-max(indices))
         ];
-      
 
 
 // Function: list_insert()
@@ -465,8 +438,6 @@ function list_insert(list, indices, values, _i=0) =
         ];
 
 
-
-
 // Function: list_remove()
 // Usage:
 //   list_remove(list, indices)
@@ -494,8 +465,6 @@ function list_remove(list, indices) =
         ]; 
 
 
-
-
 // Function: list_remove_values()
 // Usage:
 //   list_remove_values(list,values,all=false) =
@@ -565,8 +534,6 @@ function list_bset(indexset, valuelist, dflt=0) =
     );
 
 
-
-
 // Section: List Length Manipulation
 
 // Function: list_shortest()
@@ -579,7 +546,6 @@ function list_shortest(array) =
     min([for (v = array) len(v)]);
 
 
-
 // Function: list_longest()
 // Description:
 //   Returns the length of the longest sublist in a list of lists.
@@ -629,7 +595,6 @@ function list_fit(array, length, fill) =
               : list_pad(array,length,fill);
 
 
-
 // Section: List Shuffling and Sorting
 
 // Function: shuffle()
@@ -684,6 +649,7 @@ function _sort_vectors2(arr) =
     ) 
     concat( _sort_vectors2(lesser), equal, _sort_vectors2(greater) );
 
+
 // Sort a vector of vectors based on the first three entries of each vector
 // Lexicographic order, remaining entries of vector ignored
 function _sort_vectors3(arr) =
@@ -711,7 +677,6 @@ function _sort_vectors3(arr) =
     ) concat( _sort_vectors3(lesser), equal, _sort_vectors3(greater) );
 
 
-
 // Sort a vector of vectors based on the first four entries of each vector
 // Lexicographic order, remaining entries of vector ignored
 function _sort_vectors4(arr) =
@@ -742,45 +707,38 @@ function _sort_vectors4(arr) =
                                                     && y[3]>pivot[3] )))))) 
                     y ]
     ) concat( _sort_vectors4(lesser), equal, _sort_vectors4(greater) );
-    
 
 
+// when idx==undef, returns the sorted array
+// otherwise, returns the indices of the sorted array
 function _sort_general(arr, idx=undef) =
     (len(arr)<=1) ? arr :
+    is_undef(idx) 
+    ? _sort_scalar(arr)
+    : let( arrind=[for(k=[0:len(arr)-1], ark=[arr[k]]) [ k, [for (i=idx) ark[i]] ] ] )
+      _indexed_sort(arrind);
+      
+// given a list of pairs, return the first element of each pair of the list sorted by the second element of the pair
+// the sorting is done using compare_vals()
+function _indexed_sort(arrind) = 
+    arrind==[] ? [] : len(arrind)==1? [arrind[0][0]] :
+    let( pivot = arrind[floor(len(arrind)/2)][1] )
     let(
-        pivot = arr[floor(len(arr)/2)],
-        pivotval = idx==undef? pivot : [for (i=idx) pivot[i]],
-        compare = 
-            is_undef(idx) ? [for(entry=arr) compare_vals(entry, pivotval) ] :
-            [ for (entry = arr) 
-                  let( val = [for (i=idx) entry[i] ] )
-                  compare_vals(val, pivotval) ] ,
-        lesser  = [ for (i = [0:1:len(arr)-1]) if (compare[i] < 0) arr[i] ],
-        equal   = [ for (i = [0:1:len(arr)-1]) if (compare[i] ==0) arr[i] ],
-        greater = [ for (i = [0:1:len(arr)-1]) if (compare[i] > 0) arr[i] ]
-    )
-    concat(_sort_general(lesser,idx), equal, _sort_general(greater,idx));
+        lesser  = [ for (entry=arrind) if (compare_vals(entry[1], pivot) <0 ) entry ],
+        equal   = [ for (entry=arrind) if (compare_vals(entry[1], pivot)==0 ) entry[0] ],
+        greater = [ for (entry=arrind) if (compare_vals(entry[1], pivot) >0 ) entry ]
+      )
+    concat(_indexed_sort(lesser), equal, _indexed_sort(greater));
+
+
+// returns true for valid index specifications idx in the interval [imin, imax) 
+// note that idx can't have any value greater or EQUAL to imax
+function _valid_idx(idx,imin,imax) =
+    is_undef(idx) 
+    || ( is_finite(idx) && idx>=imin && idx< imax )
+    || ( is_list(idx) && min(idx)>=imin && max(idx)< imax )
+    || ( valid_range(idx) && idx[0]>=imin && idx[2]< imax );
     
-function _sort_general(arr, idx=undef) =
-    (len(arr)<=1) ? arr :
-    let(
-        pivot = arr[floor(len(arr)/2)],
-        pivotval = idx==undef? pivot : [for (i=idx) pivot[i]],
-        compare = [
-            for (entry = arr) let(
-                val = idx==undef? entry : [for (i=idx) entry[i]],
-                cmp = compare_vals(val, pivotval)
-            ) cmp
-        ],
-        lesser  = [ for (i = [0:1:len(arr)-1]) if (compare[i] < 0) arr[i] ],
-        equal   = [ for (i = [0:1:len(arr)-1]) if (compare[i] ==0) arr[i] ],
-        greater = [ for (i = [0:1:len(arr)-1]) if (compare[i] > 0) arr[i] ]
-    )
-    concat(_sort_general(lesser,idx), equal, _sort_general(greater,idx));
-
-
-
-
 
 // Function: sort()
 // Usage:
@@ -799,20 +757,21 @@ function _sort_general(arr, idx=undef) =
 //   sorted = sort(l);  // Returns [2,3,8,9,12,16,23,34,37,45,89]
 function sort(list, idx=undef) =
     !is_list(list) || len(list)<=1 ? list :
-    assert( is_undef(idx) || is_finite(idx) || is_vector(idx) || is_range(idx) , "Invalid indices.")
-    is_def(idx) ? _sort_general(list,idx) :
-    let(size = array_dim(list))
-    len(size)==1 ? _sort_scalars(list) :
-    len(size)==2 && size[1] <=4 
-    ? (
-        size[1]==0 ? list :
-        size[1]==1 ? _sort_vectors1(list) :
-        size[1]==2 ? _sort_vectors2(list) :
-        size[1]==3 ? _sort_vectors3(list)
-   /*size[1]==4*/  : _sort_vectors4(list)
-      ) 
-    : _sort_general(list);
-
+    is_def(idx) 
+    ?   assert( _valid_idx(idx,0,len(list)) , "Invalid indices.")
+        let( sarr = _sort_general(list,idx) )
+        [for(i=[0:len(sarr)-1]) list[sarr[i]] ] 
+    :   let(size = array_dim(list))
+        len(size)==1 ? _sort_scalars(list) :
+        len(size)==2 && size[1] <=4 
+        ? (
+            size[1]==0 ? list :
+            size[1]==1 ? _sort_vectors1(list) :
+            size[1]==2 ? _sort_vectors2(list) :
+            size[1]==3 ? _sort_vectors3(list)
+       /*size[1]==4*/  : _sort_vectors4(list)
+          ) 
+        : _sort_general(list);
 
 
 // Function: sortidx()
@@ -838,13 +797,13 @@ function sort(list, idx=undef) =
 //   idxs3 = sortidx(lst, idx=[1,3]); // Returns: [3,0,2,1]
 function sortidx(list, idx=undef) =
     assert( is_list(list) || is_string(list) , "Invalid input to sort." )
-    assert( is_undef(idx) || is_finite(idx) || is_vector(idx) , "Invalid indices.")
+    assert( _valid_idx(idx,0,len(list)) , "Invalid indices.")
     list==[] ? [] : 
     let(
         size = array_dim(list),
         aug = is_undef(idx) && (len(size) == 1 || (len(size) == 2 && size[1]<=4))
               ?  zip(list, list_range(len(list)))
-              :  enumerate(list,idx=idx)
+              :  0
     )
     is_undef(idx) && len(size) == 1? subindex(_sort_vectors1(aug),1) :
     is_undef(idx) && len(size) == 2 && size[1] <=4
@@ -856,25 +815,8 @@ function sortidx(list, idx=undef) =
     /*size[1]==4*/ : subindex(_sort_vectors4(aug),4)
       ) 
     :   // general case
-        subindex(_sort_general(aug, idx=list_range(s=1,n=len(aug)-1)), 0);
+        _sort_general(list,idx);
 
-function sortidx(list, idx=undef) =
-    list==[] ? [] : let(
-        size = array_dim(list),
-        aug = is_undef(idx) && (len(size) == 1 || (len(size) == 2 && size[1]<=4))?
-            zip(list, list_range(len(list))) :
-            enumerate(list,idx=idx)
-    )
-    is_undef(idx) && len(size) == 1? subindex(_sort_vectors1(aug),1) :
-    is_undef(idx) && len(size) == 2 && size[1] <=4? (
-        size[1]==0? list_range(len(arr)) :
-        size[1]==1? subindex(_sort_vectors1(aug),1) :
-        size[1]==2? subindex(_sort_vectors2(aug),2) :
-        size[1]==3? subindex(_sort_vectors3(aug),3) :
-        /*size[1]==4*/ subindex(_sort_vectors4(aug),4)
-    ) :
-    // general case
-    subindex(_sort_general(aug, idx=list_range(s=1,n=len(aug)-1)), 0);
 
 // sort() does not accept strings but sortidx does; isn't inconsistent ?
 
@@ -896,7 +838,6 @@ function unique(arr) =
     ];
 
 
-
 // Function: unique_count()
 // Usage:
 //   unique_count(arr);
@@ -913,8 +854,6 @@ function unique_count(arr) =
       [ select(arr,ind), deltas( concat(ind,[len(arr)]) ) ];
 
 
-
-
 // Section: List Iteration Helpers
 
 // Function: idx()
@@ -952,10 +891,10 @@ function idx(list, step=1, end=-1,start=0) =
 //   for (p=enumerate(colors)) right(20*p[0]) color(p[1]) circle(d=10);
 function enumerate(l,idx=undef) =
     assert(is_list(l)||is_string(list), "Invalid input." )
-    assert(is_undef(idx)||is_finite(idx)||is_vector(idx) ||is_range(idx), "Invalid index/indices." )
+    assert( _valid_idx(idx,0,len(l)), "Invalid index/indices." )
     (idx==undef)
     ?   [for (i=[0:1:len(l)-1]) [i,l[i]]]
-    :   [for (i=[0:1:len(l)-1]) concat([i], [for (j=idx) l[i][j]])];
+    :   [for (i=[0:1:len(l)-1]) [ i, for (j=idx) l[i][j]] ];
 
 
 // Function: force_list()
@@ -1109,8 +1048,6 @@ function set_union(a, b, get_indices=false) =
     ) [idxs, nset];
 
 
-
-
 // Function: set_difference()
 // Usage:
 //   s = set_difference(a, b);
@@ -1130,7 +1067,6 @@ function set_difference(a, b) =
     [ for (i=idx(a)) if(found[i]==[]) a[i] ];
 
 
-
 // Function: set_intersection()
 // Usage:
 //   s = set_intersection(a, b);
@@ -1151,7 +1087,6 @@ function set_intersection(a, b) =
 
 
 
-
 // Section: Array Manipulation
 
 // Function: add_scalar()
@@ -1170,26 +1105,60 @@ function add_scalar(v,s) =
     is_finite(s) ? [for (x=v) is_list(x)? add_scalar(x,s) : is_finite(x) ? x+s: x] : v;
 
 
-
 // Function: subindex()
+// Usage:
+//   subindex(M, idx)
 // Description:
-//   For each array item, return the indexed subitem.
-//   Returns a list of the values of each vector at the specfied
-//   index list or range.  If the index list or range has
-//   only one entry the output list is flattened.  
+//   Extracts the entries listed in idx from each entry in M.  For a matrix this means
+//   selecting a specified set of columns.  If idx is a number the return is a vector, 
+//   otherwise it is a list of lists (the submatrix).  
+//   This function will return `undef` at all entry positions indexed by idx not found in the input list M.
 // Arguments:
-//   v = The given list of lists.
+//   M = The given list of lists.
 //   idx = The index, list of indices, or range of indices to fetch.
 // Example:
-//   v = [[[1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]];
-//   subindex(v,2);      // Returns [3, 7, 11, 15]
-//   subindex(v,[2,1]);  // Returns [[3, 2], [7, 6], [11, 10], [15, 14]]
-//   subindex(v,[1:3]);  // Returns [[2, 3, 4], [6, 7, 8], [10, 11, 12], [14, 15, 16]]
-function subindex(v, idx) = 
-    [ for(val=v) 
-        let( value=[for(i=idx) val[i]] )
-        len(value)==1 ? value[0] : value
-    ];
+//   M = [[1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]];
+//   subindex(M,2);      // Returns [3, 7, 11, 15]
+//   subindex(M,[2]);    // Returns [[3], [7], [11], [15]]
+//   subindex(M,[2,1]);  // Returns [[3, 2], [7, 6], [11, 10], [15, 14]]
+//   subindex(M,[1:3]);  // Returns [[2, 3, 4], [6, 7, 8], [10, 11, 12], [14, 15, 16]]
+//   N = [ [1,2], [3], [4,5], [6,7,8] ];
+//   subindex(N,[0,1]);  // Returns [ [1,2], [3,undef], [4,5], [6,7] ]
+function subindex(M, idx) =
+    assert( is_list(M), "The input is not a list." )
+    assert( !is_undef(idx) && _valid_idx(idx,0,1/0), "Invalid index input." ) 
+    is_finite(idx)
+      ? [for(row=M) row[idx]]
+      : [for(row=M) [for(i=idx) row[i]]];
+
+
+// Function: submatrix()
+// Usage: submatrix(M, idx1, idx2)
+// Description:
+//   The input must be a list of lists (a matrix or 2d array).  Returns a submatrix by selecting the rows listed in idx1 and columsn listed in idx2.
+// Arguments:
+//   M = Given list of lists
+//   idx1 = rows index list or range
+//   idx2 = column index list or range
+// Example:
+//   M = [[ 1, 2, 3, 4, 5],
+//        [ 6, 7, 8, 9,10],
+//        [11,12,13,14,15],
+//        [16,17,18,19,20],
+//        [21,22,23,24,25]];
+//   submatrix(M,[1:2],[3:4]);  // Returns [[9, 10], [14, 15]]
+//   submatrix(M,[1], [3,4]));  // Returns [[9,10]]
+//   submatrix(M,1, [3,4]));  // Returns [[9,10]]
+//   submatrix(M,1,3));  // Returns [[9]]
+//   submatrix(M, [3,4],1); // Returns  [[17],[22]]);
+//   submatrix(M, [1,3],[2,4]); // Returns [[8,10],[18,20]]);
+//   A = [[true,    17, "test"],
+//        [[4,2],   91, false],
+//        [6,    [3,4], undef]];
+//   submatrix(A,[0,2],[1,2]);   // Returns [[17, "test"], [[3, 4], undef]]
+
+function submatrix(M,idx1,idx2) =
+    [for(i=idx1) [for(j=idx2) M[i][j] ] ];
 
 
 // Function: zip()
@@ -1318,6 +1287,10 @@ function array_dim(v, depth=undef) =
 
 // Function: transpose()
 // Description: Returns the transposition of the given array.
+//    When reverse=true, the transposition is done in respect to the secondary diagonal, that is:
+//    .
+//    reverse(transpose(reverse(arr))) == transpose(arr, reverse=true)
+//    By default, reverse=false.
 // Example:
 //   arr = [
 //       ["a", "b", "c"],
@@ -1344,16 +1317,32 @@ function array_dim(v, depth=undef) =
 //   //     ["c", "f"],
 //   // ]
 // Example:
+//   arr = [
+//       ["a", "b", "c"],
+//       ["d", "e", "f"],
+//       ["g", "h", "i"]
+//   ];
+//   t = transpose(arr, reverse=true);
+//   // Returns:
+//   // [
+//   //  ["i", "f", "c"],
+//   //  ["h", "e", "b"],
+//   //  ["g", "d", "a"]
+//   // ]
+// Example:
 //   transpose([3,4,5]);  // Returns: [3,4,5]
-function transpose(arr) =
-    let( a0 = arr[0] )
-    is_list(a0)
-    ?   assert([for(a=arr) if(len(a)!=len(a0)) 1]==[], "The array is not a matrix." )
-        [for (i=[0:1:len(a0)-1]) 
-          [ for (j=[0:1:len(arr)-1]) arr[j][i] ] ] 
-    :  arr;
-
-
+function transpose(arr, reverse=false) =
+    assert( is_list(arr) && len(arr)>0, "The array is not a vector neither a matrix." )
+    is_list(arr[0])
+    ?   let( l0 = len(arr[0]) )
+        assert([for(a=arr) if(!is_list(a) || len(a)!=l0) 1 ]==[], "The array is not a vector neither a matrix." )
+        reverse
+        ? [for (i=[0:1:l0-1]) 
+              [ for (j=[0:1:len(arr)-1]) arr[len(arr)-1-j][l0-1-i] ] ] 
+        : [for (i=[0:1:l0-1]) 
+              [ for (j=[0:1:len(arr)-1]) arr[j][i] ] ] 
+    :  assert( is_vector(arr), "The array is not a vector neither a matrix." )
+           arr;
 
 
 // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap
diff --git a/attachments.scad b/attachments.scad
index 8775832..48ed522 100644
--- a/attachments.scad
+++ b/attachments.scad
@@ -285,7 +285,7 @@ function attach_geom_size(geom) =
         ) [2*maxxr,2*maxyr,l]
     ) : type == "spheroid"? ( //r
         let( r=geom[1] )
-        is_num(r)? [2,2,2]*r : vmul([2,2,2],r)
+        is_num(r)? [2,2,2]*r : vmul([2,2,2],point3d(r))
     ) : type == "vnf_extent" || type=="vnf_isect"? ( //vnf
         let(
             mm = pointlist_bounds(geom[1][0]),
@@ -298,7 +298,7 @@ function attach_geom_size(geom) =
         ) [maxx, size.y]
     ) : type == "circle"? ( //r
         let( r=geom[1] )
-        is_num(r)? [2,2]*r : vmul([2,2],r)
+        is_num(r)? [2,2]*r : vmul([2,2],point2d(r))
     ) : type == "path_isect" || type == "path_extent"? ( //path
         let(
             mm = pointlist_bounds(geom[1]),
@@ -430,8 +430,8 @@ function find_anchor(anchor, geom) =
     ) : type == "cyl"? ( //r1, r2, l, shift
         let(
             rr1=geom[1], rr2=geom[2], l=geom[3], shift=point2d(geom[4]),
-            r1 = is_num(rr1)? [rr1,rr1] : rr1,
-            r2 = is_num(rr2)? [rr2,rr2] : rr2,
+            r1 = is_num(rr1)? [rr1,rr1] : point2d(rr1),
+            r2 = is_num(rr2)? [rr2,rr2] : point2d(rr2),
             u = (anchor.z+1)/2,
             axy = unit(point2d(anchor),[0,0]),
             bot = point3d(vmul(r1,axy), -l/2),
@@ -447,9 +447,9 @@ function find_anchor(anchor, geom) =
     ) : type == "spheroid"? ( //r
         let(
             rr = geom[1],
-            r = is_num(rr)? [rr,rr,rr] : rr,
+            r = is_num(rr)? [rr,rr,rr] : point3d(rr),
             anchor = unit(point3d(anchor),CENTER),
-            pos = point3d(cp) + vmul(r,anchor) + offset,
+            pos = point3d(cp) + vmul(r,anchor) + point3d(offset),
             vec = unit(vmul(r,anchor),UP)
         ) [anchor, pos, vec, oang]
     ) : type == "vnf_isect"? ( //vnf
@@ -458,10 +458,9 @@ function find_anchor(anchor, geom) =
             eps = 1/2048,
             points = vnf[0],
             faces = vnf[1],
-            rpts = rot(from=anchor, to=RIGHT, p=move(point3d(-cp), p=points)),
+            rpts = apply(rot(from=anchor, to=RIGHT) * move(point3d(-cp)), points),
             hits = [
-                for (i = idx(faces)) let(
-                    face = faces[i],
+                for (face = faces) let(
                     verts = select(rpts, face)
                 ) if (
                     max(subindex(verts,0)) >= -eps &&
@@ -470,35 +469,40 @@ function find_anchor(anchor, geom) =
                     min(subindex(verts,1)) <=  eps &&
                     min(subindex(verts,2)) <=  eps
                 ) let(
-                    pt = polygon_line_intersection(
-                        select(points, face),
-                        [CENTER,anchor], eps=eps
-                    )
-                ) if (!is_undef(pt)) [norm(pt), i, pt]
+                    poly = select(points, face),
+                    pt = polygon_line_intersection(poly, [cp,cp+anchor], bounded=[true,false], eps=eps)
+                ) if (!is_undef(pt)) let(
+                    plane = plane_from_polygon(poly),
+                    n = unit(plane_normal(plane))
+                )
+                [norm(pt-cp), n, pt]
             ]
         )
         assert(len(hits)>0, "Anchor vector does not intersect with the shape.  Attachment failed.")
         let(
             furthest = max_index(subindex(hits,0)),
-            pos = point3d(cp) + hits[furthest][2],
             dist = hits[furthest][0],
-            nfaces = [for (hit = hits) if(approx(hit[0],dist,eps=eps)) hit[1]],
-            n = unit(
-                sum([
-                    for (i = nfaces) let(
-                        faceverts = select(points, faces[i]),
-                        faceplane = plane_from_points(faceverts),
-                        nrm = plane_normal(faceplane)
-                    ) nrm
-                ]) / len(nfaces),
-                UP
-            )
+            pos = hits[furthest][2],
+            hitnorms = [for (hit = hits) if (approx(hit[0],dist,eps=eps)) hit[1]],
+            unorms = len(hitnorms) > 7
+              ? unique([for (nn = hitnorms) quant(nn,1e-9)])
+              : [
+                    for (i = idx(hitnorms)) let(
+                        nn = hitnorms[i],
+                        isdup = [
+                            for (j = [i+1:1:len(hitnorms)-1])
+                            if (approx(nn, hitnorms[j])) 1
+                        ] != []
+                    ) if (!isdup) nn
+                ],
+            n = unit(sum(unorms)),
+            oang = approx(point2d(n), [0,0])? 0 : atan2(n.y, n.x) + 90
         )
         [anchor, pos, n, oang]
     ) : type == "vnf_extent"? ( //vnf
         let(
             vnf=geom[1],
-            rpts = rot(from=anchor, to=RIGHT, p=move(point3d(-cp), p=vnf[0])),
+            rpts = apply(rot(from=anchor, to=RIGHT) * move(point3d(-cp)), vnf[0]),
             maxx = max(subindex(rpts,0)),
             idxs = [for (i = idx(rpts)) if (approx(rpts[i].x, maxx)) i],
             mm = pointlist_bounds(select(rpts,idxs)),
@@ -519,10 +523,10 @@ function find_anchor(anchor, geom) =
     ) : type == "circle"? ( //r
         let(
             rr = geom[1],
-            r = is_num(rr)? [rr,rr] : rr,
-            pos = point2d(cp) + vmul(r,anchor) + offset,
+            r = is_num(rr)? [rr,rr] : point2d(rr),
             anchor = unit(point2d(anchor),[0,0]),
-            vec = unit(vmul([r.y,r.x],anchor),[0,1])
+            pos = point2d(cp) + vmul(r,anchor) + point2d(offset),
+            vec = unit(vmul(r,anchor),[0,1])
         ) [anchor, pos, vec, 0]
     ) : type == "path_isect"? ( //path
         let(
@@ -849,7 +853,7 @@ module attachable(
 
 // Module: position()
 // Usage:
-//   position(from, [overlap]) ...
+//   position(from) ...
 // Description:
 //   Attaches children to a parent object at an anchor point.
 // Arguments:
@@ -970,7 +974,8 @@ module edge_profile(edges=EDGES_ALL, except=[], convexity=10) {
         $attach_anchor = anch;
         $attach_norot = true;
         $tags = "mask";
-        length = sum(vmul($parent_size, [for (x=vec) x?0:1]))+0.1;
+        psize = point3d($parent_size);
+        length = [for (i=[0:2]) if(!vec[i]) psize[i]][0]+0.1;
         rotang =
             vec.z<0? [90,0,180+vang(point2d(vec))] :
             vec.z==0 && sign(vec.x)==sign(vec.y)? 135+vang(point2d(vec)) :
@@ -1230,17 +1235,20 @@ module show(tags="")
 //   }
 module diff(neg, pos=undef, keep=undef)
 {
-    difference() {
-        if (pos != undef) {
-            show(pos) children();
-        } else {
-            if (keep == undef) {
-                hide(neg) children();
+    // Don't perform the operation if the current tags are hidden
+    if (attachment_is_shown($tags)) {
+        difference() {
+            if (pos != undef) {
+                show(pos) children();
             } else {
-                hide(str(neg," ",keep)) children();
+                if (keep == undef) {
+                    hide(neg) children();
+                } else {
+                    hide(str(neg," ",keep)) children();
+                }
             }
+            show(neg) children();
         }
-        show(neg) children();
     }
     if (keep!=undef) {
         show(keep) children();
@@ -1275,17 +1283,20 @@ module diff(neg, pos=undef, keep=undef)
 //   }
 module intersect(a, b=undef, keep=undef)
 {
-    intersection() {
-        if (b != undef) {
-            show(b) children();
-        } else {
-            if (keep == undef) {
-                hide(a) children();
+    // Don't perform the operation if the current tags are hidden
+    if (attachment_is_shown($tags)) {
+        intersection() {
+            if (b != undef) {
+                show(b) children();
             } else {
-                hide(str(a," ",keep)) children();
+                if (keep == undef) {
+                    hide(a) children();
+                } else {
+                    hide(str(a," ",keep)) children();
+                }
             }
+            show(a) children();
         }
-        show(a) children();
     }
     if (keep!=undef) {
         show(keep) children();
diff --git a/common.scad b/common.scad
index f0be6c1..db4f3bb 100644
--- a/common.scad
+++ b/common.scad
@@ -129,18 +129,13 @@ function is_list_of(list,pattern) =
     is_list(list) &&
     []==[for(entry=0*list) if (entry != pattern) entry];
 
-function _list_pattern(list) =
-  is_list(list) ? [for(entry=list) is_list(entry) ? _list_pattern(entry) : 0]
-                : 0;
-
-
 
 // Function: is_consistent()
 // Usage:
 //   is_consistent(list)
 // Description:
 //   Tests whether input is a list of entries which all have the same list structure
-//   and are filled with finite numerical data.
+//   and are filled with finite numerical data. It returns `true`for the empty list. 
 // Example:
 //   is_consistent([3,4,5]);              // Returns true
 //   is_consistent([[3,4],[4,5],[6,7]]);  // Returns true
@@ -148,7 +143,7 @@ function _list_pattern(list) =
 //   is_consistent([[3,[3,4,[5]]], [5,[2,9,[9]]]]); // Returns true
 //   is_consistent([[3,[3,4,[5]]], [5,[2,9,9]]]);   // Returns false
 function is_consistent(list) =
-  is_list(list) && is_list_of(list, _list_pattern(list[0]));
+  /*is_list(list) &&*/ is_list_of(list, _list_pattern(list[0]));
 
 
 //Internal function
@@ -198,11 +193,11 @@ function first_defined(v,recursive=false,_i=0) =
             is_undef(first_defined(v[_i],recursive=recursive))
         )
     )? first_defined(v,recursive=recursive,_i=_i+1) : v[_i];
-
+    
 
 // Function: one_defined()
 // Usage:
-//   one_defined(vars, names, [required])
+//   one_defined(vars, names, <required>)
 // Description:
 //   Examines the input list `vars` and returns the entry which is not `undef`.  If more
 //   than one entry is `undef` then issues an assertion specifying "Must define exactly one of" followed
@@ -221,8 +216,7 @@ function one_defined(vars, names, required=true) =
 
 // Function: num_defined()
 // Description: Counts how many items in list `v` are not `undef`.
-function num_defined(v,_i=0,_cnt=0) = _i>=len(v)? _cnt : num_defined(v,_i+1,_cnt+(is_undef(v[_i])? 0 : 1));
-
+function num_defined(v) = len([for(vi=v) if(!is_undef(vi)) 1]);
 
 // Function: any_defined()
 // Description:
@@ -239,8 +233,8 @@ function any_defined(v,recursive=false) = first_defined(v,recursive=recursive) !
 // Arguments:
 //   v = The list whose items are being checked.
 //   recursive = If true, any sublists are evaluated recursively.
-function all_defined(v,recursive=false) = max([for (x=v) is_undef(x)||(recursive&&is_list(x)&&!all_defined(x))? 1 : 0])==0;
-
+function all_defined(v,recursive=false) = 
+    []==[for (x=v) if(is_undef(x)||(recursive && is_list(x) && !all_defined(x,recursive))) 0 ];
 
 
 
@@ -249,7 +243,7 @@ function all_defined(v,recursive=false) = max([for (x=v) is_undef(x)||(recursive
 
 // Function: get_anchor()
 // Usage:
-//   get_anchor(anchor,center,[uncentered],[dflt]);
+//   get_anchor(anchor,center,<uncentered>,<dflt>);
 // Description:
 //   Calculated the correct anchor from `anchor` and `center`.  In order:
 //   - If `center` is not `undef` and `center` evaluates as true, then `CENTER` (`[0,0,0]`) is returned.
@@ -270,7 +264,7 @@ function get_anchor(anchor,center,uncentered=BOT,dflt=CENTER) =
 
 // Function: get_radius()
 // Usage:
-//   get_radius([r1], [r2], [r], [d1], [d2], [d], [dflt]);
+//   get_radius(<r1>, <r2>, <r>, <d1>, <d2>, <d>, <dflt>);
 // Description:
 //   Given various radii and diameters, returns the most specific radius.
 //   If a diameter is most specific, returns half its value, giving the radius.
@@ -288,19 +282,23 @@ function get_anchor(anchor,center,uncentered=BOT,dflt=CENTER) =
 //   r = Most general radius.
 //   d = Most general diameter.
 //   dflt = Value to return if all other values given are `undef`.
-function get_radius(r1=undef, r2=undef, r=undef, d1=undef, d2=undef, d=undef, dflt=undef) = (
-    !is_undef(r1)? assert(is_undef(r2)&&is_undef(d1)&&is_undef(d2), "Conflicting or redundant radius/diameter arguments given.") r1 :
-    !is_undef(r2)? assert(is_undef(d1)&&is_undef(d2), "Conflicting or redundant radius/diameter arguments given.") r2 :
-    !is_undef(d1)? d1/2 :
-    !is_undef(d2)? d2/2 :
-    !is_undef(r)? assert(is_undef(d), "Conflicting or redundant radius/diameter arguments given.") r :
-    !is_undef(d)? d/2 :
-    dflt
-);
+function get_radius(r1=undef, r2=undef, r=undef, d1=undef, d2=undef, d=undef, dflt=undef) = 
+    assert(num_defined([r1,d1,r2,d2])<2, "Conflicting or redundant radius/diameter arguments given.")
+    !is_undef(r1) ?   assert(is_finite(r1), "Invalid radius r1." ) r1 
+    : !is_undef(r2) ? assert(is_finite(r2), "Invalid radius r2." ) r2
+    : !is_undef(d1) ? assert(is_finite(d1), "Invalid diameter d1." ) d1/2
+    : !is_undef(d2) ? assert(is_finite(d2), "Invalid diameter d2." ) d2/2
+    : !is_undef(r)
+      ? assert(is_undef(d), "Conflicting or redundant radius/diameter arguments given.")
+        assert(is_finite(r) || is_vector(r,1) || is_vector(r,2), "Invalid radius r." )
+        r 
+    : !is_undef(d) ? assert(is_finite(d) || is_vector(d,1) || is_vector(d,2), "Invalid diameter d." ) d/2
+    : dflt;
+
 
 // Function: get_height()
 // Usage:
-//   get_height([h],[l],[height],[dflt])
+//   get_height(<h>,<l>,<height>,<dflt>)
 // Description:
 //   Given several different parameters for height check that height is not multiply defined
 //   and return a single value.  If the three values `l`, `h`, and `height` are all undefined
@@ -317,7 +315,7 @@ function get_height(h=undef,l=undef,height=undef,dflt=undef) =
 
 // Function: scalar_vec3()
 // Usage:
-//   scalar_vec3(v, [dflt]);
+//   scalar_vec3(v, <dflt>);
 // Description:
 //   If `v` is a scalar, and `dflt==undef`, returns `[v, v, v]`.
 //   If `v` is a scalar, and `dflt!=undef`, returns `[v, dflt, dflt]`.
@@ -346,6 +344,19 @@ function segs(r) =
 
 
 
+// Module: no_children()
+// Usage:
+//   no_children($children);
+// Description:
+//   Assert that the calling module does not support children.  Prints an error message to this effect and fails if children are present,
+//   as indicated by its argument.
+// Arguments:
+//   $children = number of children the module has.  
+module no_children(count) {
+  assert(count==0, str("Module ",parent_module(1),"() does not support child modules"));
+}
+
+
 // Section: Testing Helpers
 
 
@@ -356,7 +367,7 @@ function _valstr(x) =
 
 // Module: assert_approx()
 // Usage:
-//   assert_approx(got, expected, [info]);
+//   assert_approx(got, expected, <info>);
 // Description:
 //   Tests if the value gotten is what was expected.  If not, then
 //   the expected and received values are printed to the console and
@@ -383,7 +394,7 @@ module assert_approx(got, expected, info) {
 
 // Module: assert_equal()
 // Usage:
-//   assert_equal(got, expected, [info]);
+//   assert_equal(got, expected, <info>);
 // Description:
 //   Tests if the value gotten is what was expected.  If not, then
 //   the expected and received values are printed to the console and
@@ -410,7 +421,7 @@ module assert_equal(got, expected, info) {
 
 // Module: shape_compare()
 // Usage:
-//   shape_compare([eps]) {test_shape(); expected_shape();}
+//   shape_compare(<eps>) {test_shape(); expected_shape();}
 // Description:
 //   Compares two child shapes, returning empty geometry if they are very nearly the same shape and size.
 //   Returns the differential geometry if they are not nearly the same shape and size.
diff --git a/distributors.scad b/distributors.scad
index 5738033..787693e 100644
--- a/distributors.scad
+++ b/distributors.scad
@@ -949,7 +949,7 @@ module ovoid_spread(r=undef, d=undef, n=100, cone_ang=90, scale=[1,1,1], perp=tr
     for ($idx = idx(theta_phis)) {
         tp = theta_phis[$idx];
         xyz = spherical_to_xyz(r, tp[0], tp[1]);
-        $pos = vmul(xyz,scale);
+        $pos = vmul(xyz,point3d(scale,1));
         $theta = tp[0];
         $phi = tp[1];
         $rad = r;
diff --git a/errors.scad b/errors.scad
deleted file mode 100644
index 8c174d9..0000000
--- a/errors.scad
+++ /dev/null
@@ -1,126 +0,0 @@
-//////////////////////////////////////////////////////////////////////
-// LibFile: errors.scad
-//   Functions and modules to facilitate error reporting.
-//   To use, include this line at the top of your file:
-//   ```
-//   use <BOSL2/std.scad>
-//   ```
-//////////////////////////////////////////////////////////////////////
-
-
-
-// Section: Warnings and Errors
-
-
-// Module: no_children()
-// Usage:
-//   no_children($children);
-// Description:
-//   Assert that the calling module does not support children.  Prints an error message to this effect and fails if children are present,
-//   as indicated by its argument.
-// Arguments:
-//   $children = number of children the module has.  
-module no_children(count) {
-  assert(count==0, str("Module ",parent_module(1),"() does not support child modules"));
-}
-
-
-// Function&Module: echo_error()
-// Usage:
-//   echo_error(msg, [pfx]);
-// Description:
-//   Emulates printing of an error message.  The text will be shaded red.
-//   You can also use this as a function call from a function.
-// Arguments:
-//   msg = The message to print.
-//   pfx = The prefix to print before `msg`.  Default: `ERROR`
-module echo_error(msg, pfx="ERROR") {
-    echo(str("<p style=\"background-color: #ffb0b0\"><b>", pfx, ":</b> ", msg, "</p>"));
-}
-
-function echo_error(msg, pfx="ERROR") =
-    echo(str("<p style=\"background-color: #ffb0b0\"><b>", pfx, ":</b> ", msg, "</p>"));
-
-
-// Function&Module: echo_warning()
-// Usage:
-//   echo_warning(msg, [pfx]);
-// Description:
-//   Emulates printing of a warning message.  The text will be shaded yellow.
-//   You can also use this as a function call from a function.
-// Arguments:
-//   msg = The message to print.
-//   pfx = The prefix to print before `msg`.  Default: `WARNING`
-module echo_warning(msg, pfx="WARNING") {
-    echo(str("<p style=\"background-color: #ffffb0\"><b>", pfx, ":</b> ", msg, "</p>"));
-}
-
-function echo_warning(msg, pfx="WARNING") =
-    echo(str("<p style=\"background-color: #ffffb0\"><b>", pfx, ":</b> ", msg, "</p>"));
-
-
-// Function&Module: deprecate()
-// Usage:
-//   deprecate(name, [suggest]);
-// Description:
-//   Show module deprecation warnings.
-//   You can also use this as a function call from a function.
-// Arguments:
-//   name = The name of the module that is deprecated.
-//   suggest = If given, the module to recommend using instead.
-module deprecate(name, suggest=undef) {
-    echo_warning(pfx="DEPRECATED",
-        str(
-            "`<code>", name, "</code>` is deprecated and should not be used.",
-            is_undef(suggest)? "" : str(
-                "  You should use `<code>", suggest, "</code>` instead."
-            )
-        )
-    );
-}
-
-function deprecate(name, suggest=undef) =
-    echo_warning(pfx="DEPRECATED",
-        str(
-            "`<code>", name, "</code>` is deprecated and should not be used.",
-            is_undef(suggest)? "" : str(
-                "  You should use `<code>", suggest, "</code>` instead."
-            )
-        )
-    );
-
-
-// Function&Module: deprecate_argument()
-// Usage:
-//   deprecate(name, arg, [suggest]);
-// Description:
-//   Show argument deprecation warnings.
-//   You can also use this as a function call from a function.
-// Arguments:
-//   name = The name of the module/function the deprecated argument is used in.
-//   arg = The name of the deprecated argument.
-//   suggest = If given, the argument to recommend using instead.
-module deprecate_argument(name, arg, suggest=undef) {
-    echo_warning(pfx="DEPRECATED ARG", str(
-        "In `<code>", name, "</code>`, ",
-        "the argument `<code>", arg, "</code>` ",
-        "is deprecated and should not be used.",
-        is_undef(suggest)? "" : str(
-            "  You should use `<code>", suggest, "</code>` instead."
-        )
-    ));
-}
-
-function deprecate_argument(name, arg, suggest=undef) =
-    echo_warning(pfx="DEPRECATED ARG", str(
-        "In `<code>", name, "</code>`, ",
-        "the argument `<code>", arg, "</code>` ",
-        "is deprecated and should not be used.",
-        is_undef(suggest)? "" : str(
-            "  You should use `<code>", suggest, "</code>` instead."
-        )
-    ));
-
-
-
-// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap
diff --git a/geometry.scad b/geometry.scad
index 443c026..dc86187 100644
--- a/geometry.scad
+++ b/geometry.scad
@@ -21,85 +21,71 @@
 //   edge = Array of two points forming the line segment to test against.
 //   eps = Acceptable variance.  Default: `EPSILON` (1e-9)
 function point_on_segment2d(point, edge, eps=EPSILON) =
-    approx(point,edge[0],eps=eps) || approx(point,edge[1],eps=eps) ||  // The point is an endpoint
-    sign(edge[0].x-point.x)==sign(point.x-edge[1].x)  // point is in between the
-        && sign(edge[0].y-point.y)==sign(point.y-edge[1].y)  // edge endpoints
-        && approx(point_left_of_segment2d(point, edge),0,eps=eps);  // and on the line defined by edge
+    assert( is_vector(point,2), "Invalid point." )
+    assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." )
+    assert( _valid_line(edge,2,eps=eps), "Invalid segment." )
+    let( dp = point-edge[0],
+         de = edge[1]-edge[0],
+         ne = norm(de) ) 
+    ( dp*de >= -eps*ne ) 
+    && ( (dp-de)*de <= eps*ne )                  // point projects on the segment
+    && _dist2line(point-edge[0],unit(de))<eps;   // point is on the line 
+    
+    
+//Internal - distance from point `d` to the line passing through the origin with unit direction n
+//_dist2line works for any dimension
+function _dist2line(d,n) = norm(d-(d * n) * n);
+
+// Internal non-exposed function.
+function _point_above_below_segment(point, edge) =
+    let( edge = edge - [point, point] )
+    edge[0].y <= 0 
+    ?   (edge[1].y >  0 && cross(edge[0], edge[1]-edge[0]) > 0) ?  1 : 0
+    :   (edge[1].y <= 0 && cross(edge[0], edge[1]-edge[0]) < 0) ? -1 : 0 ;
+
+//Internal
+function _valid_line(line,dim,eps=EPSILON) = 
+    is_matrix(line,2,dim) 
+    && ! approx(norm(line[1]-line[0]), 0, eps); 
+    
+//Internal
+function _valid_plane(p, eps=EPSILON) = is_vector(p,4) && ! approx(norm(p),0,eps);
 
 
-// Function: point_left_of_segment2d()
+// Function: point_left_of_line2d()
 // Usage:
-//   point_left_of_segment2d(point, edge);
+//   point_left_of_line2d(point, line);
 // Description:
-//   Return >0 if point is left of the line defined by edge.
+//   Return >0 if point is left of the line defined by `line`.
 //   Return =0 if point is on the line.
 //   Return <0 if point is right of the line.
 // Arguments:
 //   point = The point to check position of.
-//   edge = Array of two points forming the line segment to test against.
-function point_left_of_segment2d(point, edge) =
-    (edge[1].x-edge[0].x) * (point.y-edge[0].y) - (point.x-edge[0].x) * (edge[1].y-edge[0].y);
-
-
-// Internal non-exposed function.
-function _point_above_below_segment(point, edge) =
-    edge[0].y <= point.y? (
-        (edge[1].y > point.y && point_left_of_segment2d(point, edge) > 0)? 1 : 0
-    ) : (
-        (edge[1].y <= point.y && point_left_of_segment2d(point, edge) < 0)? -1 : 0
-    );
+//   line  = Array of two points forming the line segment to test against.
+function point_left_of_line2d(point, line) =
+    assert( is_vector(point,2) && is_vector(line*point, 2), "Improper input." )
+    cross(line[0]-point, line[1]-line[0]);
 
 
 // Function: collinear()
 // Usage:
-//   collinear(a, b, c, [eps]);
+//   collinear(a, [b, c], [eps]);
 // Description:
-//   Returns true if three points are co-linear.
+//   Returns true if the points `a`, `b` and `c` are co-linear or if the list of points `a` is collinear.
 // Arguments:
-//   a = First point.
-//   b = Second point.
-//   c = Third point.
+//   a = First point or list of points.
+//   b = Second point or undef; it should be undef if `c` is undef
+//   c = Third point or undef.
 //   eps = Acceptable variance.  Default: `EPSILON` (1e-9)
 function collinear(a, b, c, eps=EPSILON) =
-    approx(a,b,eps=eps)? true :
-    distance_from_line([a,b], c) < eps;
-
-
-// Function: collinear_indexed()
-// Usage:
-//   collinear_indexed(points, a, b, c, [eps]);
-// Description:
-//   Returns true if three points are co-linear.
-// Arguments:
-//   points = A list of points.
-//   a = Index in `points` of first point.
-//   b = Index in `points` of second point.
-//   c = Index in `points` of third point.
-//   eps = Acceptable max angle variance.  Default: EPSILON (1e-9) degrees.
-function collinear_indexed(points, a, b, c, eps=EPSILON) =
-    let(
-        p1=points[a],
-        p2=points[b],
-        p3=points[c]
-    ) collinear(p1, p2, p3, eps);
-
-
-// Function: points_are_collinear()
-// Usage:
-//   points_are_collinear(points);
-// Description:
-//   Given a list of points, returns true if all points in the list are collinear.
-// Arguments:
-//   points = The list of points to test.
-//   eps = How much variance is allowed in testing that each point is on the same line.  Default: `EPSILON` (1e-9)
-function points_are_collinear(points, eps=EPSILON) =
-    let(
-        a = furthest_point(points[0], points),
-        b = furthest_point(points[a], points),
-        pa = points[a],
-        pb = points[b]
-    ) all([for (pt = points) collinear(pa, pb, pt, eps=eps)]);
-
+    assert( is_path([a,b,c],dim=undef)
+            || ( is_undef(b) && is_undef(c) && is_path(a,dim=undef) ), 
+            "Input should be 3 points or a list of points with same dimension.")
+    assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." )
+    let( points = is_def(c) ? [a,b,c]: a )
+    len(points)<3 ? true
+    : noncollinear_triple(points,error=false,eps=eps)==[];
+ 
 
 // Function: distance_from_line()
 // Usage:
@@ -112,10 +98,11 @@ function points_are_collinear(points, eps=EPSILON) =
 // Example:
 //   distance_from_line([[-10,0], [10,0]], [3,8]);  // Returns: 8
 function distance_from_line(line, pt) =
-    let(a=line[0], n=unit(line[1]-a), d=a-pt)
-    norm(d - ((d * n) * n));
-
-
+    assert( _valid_line(line) && is_vector(pt,len(line[0])), 
+            "Invalid line, invalid point or incompatible dimensions." )
+    _dist2line(pt-line[0],unit(line[1]-line[0]));
+    
+    
 // Function: line_normal()
 // Usage:
 //   line_normal([P1,P2])
@@ -133,9 +120,11 @@ function distance_from_line(line, pt) =
 //   color("green") stroke([p1,p1+10*n], endcap2="arrow2");
 //   color("blue") move_copies([p1,p2]) circle(d=2, $fn=12);
 function line_normal(p1,p2) =
-    is_undef(p2)?
-        assert(is_path(p1,2)) line_normal(p1[0],p1[1]) :
-        assert(is_vector(p1,2)&&is_vector(p2,2)) unit([p1.y-p2.y,p2.x-p1.x]);
+    is_undef(p2)
+    ?   assert( len(p1)==2 && !is_undef(p1[1]) , "Invalid input." ) 
+        line_normal(p1[0],p1[1]) 
+    :   assert( _valid_line([p1,p2],dim=2), "Invalid line." ) 
+        unit([p1.y-p2.y,p2.x-p1.x]);
 
 
 // 2D Line intersection from two segments.
@@ -166,7 +155,10 @@ function _general_line_intersection(s1,s2,eps=EPSILON) =
 //   l2 = Second 2D line, given as a list of two 2D points on the line.
 //   eps = Acceptable variance.  Default: `EPSILON` (1e-9)
 function line_intersection(l1,l2,eps=EPSILON) =
-    let(isect = _general_line_intersection(l1,l2,eps=eps)) isect[0];
+    assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." )
+    assert( _valid_line(l1,dim=2,eps=eps) &&_valid_line(l2,dim=2,eps=eps), "Invalid line(s)." )
+    let(isect = _general_line_intersection(l1,l2,eps=eps)) 
+    isect[0];
 
 
 // Function: line_ray_intersection()
@@ -180,9 +172,12 @@ function line_intersection(l1,l2,eps=EPSILON) =
 //   ray = The 2D ray, given as a list `[START,POINT]` of the 2D start-point START, and a 2D point POINT on the ray.
 //   eps = Acceptable variance.  Default: `EPSILON` (1e-9)
 function line_ray_intersection(line,ray,eps=EPSILON) =
+    assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." )
+    assert( _valid_line(line,dim=2,eps=eps) && _valid_line(ray,dim=2,eps=eps), "Invalid line or ray." )
     let(
         isect = _general_line_intersection(line,ray,eps=eps)
-    ) isect[2]<0-eps? undef : isect[0];
+    ) 
+    (isect[2]<0-eps) ? undef : isect[0];
 
 
 // Function: line_segment_intersection()
@@ -196,6 +191,8 @@ function line_ray_intersection(line,ray,eps=EPSILON) =
 //   segment = The bounded 2D line segment, given as a list of the two 2D endpoints of the segment.
 //   eps = Acceptable variance.  Default: `EPSILON` (1e-9)
 function line_segment_intersection(line,segment,eps=EPSILON) =
+    assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." )
+    assert( _valid_line(line,  dim=2,eps=eps) &&_valid_line(segment,dim=2,eps=eps), "Invalid line or segment." )
     let(
         isect = _general_line_intersection(line,segment,eps=eps)
     ) isect[2]<0-eps || isect[2]>1+eps ? undef : isect[0];
@@ -212,9 +209,12 @@ function line_segment_intersection(line,segment,eps=EPSILON) =
 //   r2 = Second 2D ray, given as a list `[START,POINT]` of the 2D start-point START, and a 2D point POINT on the ray.
 //   eps = Acceptable variance.  Default: `EPSILON` (1e-9)
 function ray_intersection(r1,r2,eps=EPSILON) =
+    assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." )
+    assert( _valid_line(r1,dim=2,eps=eps) && _valid_line(r2,dim=2,eps=eps), "Invalid ray(s)." )
     let(
         isect = _general_line_intersection(r1,r2,eps=eps)
-    ) isect[1]<0-eps || isect[2]<0-eps? undef : isect[0];
+    ) 
+    isect[1]<0-eps || isect[2]<0-eps ? undef : isect[0];
 
 
 // Function: ray_segment_intersection()
@@ -228,9 +228,16 @@ function ray_intersection(r1,r2,eps=EPSILON) =
 //   segment = The bounded 2D line segment, given as a list of the two 2D endpoints of the segment.
 //   eps = Acceptable variance.  Default: `EPSILON` (1e-9)
 function ray_segment_intersection(ray,segment,eps=EPSILON) =
+    assert( _valid_line(ray,dim=2,eps=eps) && _valid_line(segment,dim=2,eps=eps), "Invalid ray or segment." )
+    assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." )
     let(
         isect = _general_line_intersection(ray,segment,eps=eps)
-    ) isect[1]<0-eps || isect[2]<0-eps || isect[2]>1+eps ? undef : isect[0];
+    ) 
+    isect[1]<0-eps 
+    || isect[2]<0-eps 
+    || isect[2]>1+eps
+    ? undef 
+    : isect[0];
 
 
 // Function: segment_intersection()
@@ -244,9 +251,17 @@ function ray_segment_intersection(ray,segment,eps=EPSILON) =
 //   s2 = Second 2D segment, given as a list of the two 2D endpoints of the line segment.
 //   eps = Acceptable variance.  Default: `EPSILON` (1e-9)
 function segment_intersection(s1,s2,eps=EPSILON) =
+    assert( _valid_line(s1,dim=2,eps=eps) && _valid_line(s2,dim=2,eps=eps), "Invalid segment(s)." )
+    assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." )
     let(
         isect = _general_line_intersection(s1,s2,eps=eps)
-    ) isect[1]<0-eps || isect[1]>1+eps || isect[2]<0-eps || isect[2]>1+eps ? undef : isect[0];
+    ) 
+    isect[1]<0-eps 
+    || isect[1]>1+eps 
+    || isect[2]<0-eps 
+    || isect[2]>1+eps 
+    ? undef 
+    : isect[0];
 
 
 // Function: line_closest_point()
@@ -301,15 +316,10 @@ function segment_intersection(s1,s2,eps=EPSILON) =
 //   color("blue") translate(pt) sphere(r=1,$fn=12);
 //   color("red") translate(p2) sphere(r=1,$fn=12);
 function line_closest_point(line,pt) =
-    assert(is_path(line)&&len(line)==2)
-    assert(same_shape(pt,line[0]))
-    assert(!approx(line[0],line[1]))
-    let(
-        seglen = norm(line[1]-line[0]),
-        segvec = (line[1]-line[0])/seglen,
-        projection = (pt-line[0]) * segvec
-    )
-    line[0] + projection*segvec;
+    assert(_valid_line(line), "Invalid line." )
+    assert( is_vector(pt,len(line[0])), "Invalid point or incompatible dimensions." )
+    let( n = unit( line[0]- line[1]) )
+    line[1]+((pt- line[1]) * n) * n;
 
 
 // Function: ray_closest_point()
@@ -364,9 +374,8 @@ function line_closest_point(line,pt) =
 //   color("blue") translate(pt) sphere(r=1,$fn=12);
 //   color("red") translate(p2) sphere(r=1,$fn=12);
 function ray_closest_point(ray,pt) =
-    assert(is_path(ray)&&len(ray)==2)
-    assert(same_shape(pt,ray[0]))
-    assert(!approx(ray[0],ray[1]))
+    assert( _valid_line(ray), "Invalid ray." )
+    assert(is_vector(pt,len(ray[0])), "Invalid point or incompatible dimensions." )
     let(
         seglen = norm(ray[1]-ray[0]),
         segvec = (ray[1]-ray[0])/seglen,
@@ -428,8 +437,8 @@ function ray_closest_point(ray,pt) =
 //   color("blue") translate(pt) sphere(r=1,$fn=12);
 //   color("red") translate(p2) sphere(r=1,$fn=12);
 function segment_closest_point(seg,pt) =
-    assert(is_path(seg)&&len(seg)==2)
-    assert(same_shape(pt,seg[0]))
+    assert(_valid_line(seg), "Invalid segment." )
+    assert(len(pt)==len(seg[0]), "Incompatible dimensions." )
     approx(seg[0],seg[1])? seg[0] :
     let(
         seglen = norm(seg[1]-seg[0]),
@@ -440,6 +449,26 @@ function segment_closest_point(seg,pt) =
     projection>=seglen ? seg[1] :
     seg[0] + projection*segvec;
 
+    
+// Function: line_from_points()
+// Usage:
+//   line_from_points(points, [fast], [eps]);
+// Description:
+//   Given a list of 2 or more colinear points, returns a line containing them.
+//   If `fast` is false and the points are coincident, then `undef` is returned.
+//   if `fast` is true, then the collinearity test is skipped and a line passing through 2 distinct arbitrary points is returned.
+// Arguments:
+//   points = The list of points to find the line through.
+//   fast = If true, don't verify that all points are collinear.  Default: false
+//   eps = How much variance is allowed in testing each point against the line.  Default: `EPSILON` (1e-9)
+function line_from_points(points, fast=false, eps=EPSILON) =
+    assert( is_path(points,dim=undef), "Improper point list." )
+    assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." )
+    let( pb = furthest_point(points[0],points) )
+    approx(norm(points[pb]-points[0]),0) ? undef :
+    fast || collinear(points) ? [points[pb], points[0]] : undef;
+
+
 
 // Section: 2D Triangles
 
@@ -486,21 +515,30 @@ function segment_closest_point(seg,pt) =
 //   ang = tri_calc(adj=20,hyp=30)[3];
 //   ang2 = tri_calc(adj=20,hyp=40)[4];
 function tri_calc(ang,ang2,adj,opp,hyp) =
-    assert(ang==undef || ang2==undef,"You cannot specify both ang and ang2.")
-    assert(num_defined([ang,ang2,adj,opp,hyp])==2, "You must specify exactly two arguments.")
+    assert(ang==undef || ang2==undef,"At most one angle is allowed.")
+    assert(num_defined([ang,ang2,adj,opp,hyp])==2, "Exactly two arguments must be given.")
     let(
-        ang = ang!=undef? assert(ang>0&&ang<90) ang :
-            ang2!=undef? (90-ang2) :
-            adj==undef? asin(constrain(opp/hyp,-1,1)) :
-            opp==undef? acos(constrain(adj/hyp,-1,1)) :
-            atan2(opp,adj),
-        ang2 = ang2!=undef? assert(ang2>0&&ang2<90) ang2 : (90-ang),
-        adj = adj!=undef? assert(adj>0) adj :
-            (opp!=undef? (opp/tan(ang)) : (hyp*cos(ang))),
-        opp = opp!=undef? assert(opp>0) opp :
-            (adj!=undef? (adj*tan(ang)) : (hyp*sin(ang))),
-        hyp = hyp!=undef? assert(hyp>0) assert(adj<hyp) assert(opp<hyp) hyp :
-            (adj!=undef? (adj/cos(ang)) : (opp/sin(ang)))
+        ang   = ang!=undef
+                ? assert(ang>0&&ang<90, "The input angles should be acute angles." ) ang 
+                : ang2!=undef ? (90-ang2) 
+                : adj==undef ? asin(constrain(opp/hyp,-1,1)) 
+                : opp==undef ? acos(constrain(adj/hyp,-1,1)) 
+                : atan2(opp,adj),
+        ang2 =  ang2!=undef
+                ? assert(ang2>0&&ang2<90, "The input angles should be acute angles." ) ang2 
+                : (90-ang),
+        adj  =  adj!=undef
+                ? assert(adj>0, "Triangle side lengths should be positive." ) adj 
+                : (opp!=undef? (opp/tan(ang)) : (hyp*cos(ang))),
+        opp  =  opp!=undef
+                ? assert(opp>0, "Triangle side lengths should be positive." )  opp 
+                : (adj!=undef? (adj*tan(ang)) : (hyp*sin(ang))),
+        hyp  =  hyp!=undef
+                ? assert(hyp>0, "Triangle side lengths should be positive." ) 
+                  assert(adj<hyp && opp<hyp, "Hyphotenuse length should be greater than the other sides." )
+                  hyp 
+                : (adj!=undef? (adj/cos(ang)) 
+                : (opp/sin(ang)))
     )
     [adj, opp, hyp, ang, ang2];
 
@@ -517,8 +555,8 @@ function tri_calc(ang,ang2,adj,opp,hyp) =
 // Example:
 //   hyp = hyp_opp_to_adj(5,3);  // Returns: 4
 function hyp_opp_to_adj(hyp,opp) =
-    assert(is_num(hyp)&&hyp>=0)
-    assert(is_num(opp)&&opp>=0)
+    assert(is_finite(hyp+opp) && hyp>=0 && opp>=0, 
+           "Triangle side lengths should be a positive numbers." )
     sqrt(hyp*hyp-opp*opp);
 
 
@@ -534,8 +572,8 @@ function hyp_opp_to_adj(hyp,opp) =
 // Example:
 //   adj = hyp_ang_to_adj(8,60);  // Returns: 4
 function hyp_ang_to_adj(hyp,ang) =
-    assert(is_num(hyp)&&hyp>=0)
-    assert(is_num(ang)&&ang>0&&ang<90)
+    assert(is_finite(hyp) && hyp>=0, "Triangle side length should be a positive number." )
+    assert(is_finite(ang) && ang>0 && ang<90, "The angle should be an acute angle." )
     hyp*cos(ang);
 
 
@@ -551,8 +589,8 @@ function hyp_ang_to_adj(hyp,ang) =
 // Example:
 //   adj = opp_ang_to_adj(8,30);  // Returns: 4
 function opp_ang_to_adj(opp,ang) =
-    assert(is_num(opp)&&opp>=0)
-    assert(is_num(ang)&&ang>0&&ang<90)
+    assert(is_finite(opp) && opp>=0, "Triangle side length should be a positive number." )
+    assert(is_finite(ang) && ang>0 && ang<90, "The angle should be an acute angle." )
     opp/tan(ang);
 
 
@@ -567,8 +605,8 @@ function opp_ang_to_adj(opp,ang) =
 // Example:
 //   opp = hyp_adj_to_opp(5,4);  // Returns: 3
 function hyp_adj_to_opp(hyp,adj) =
-    assert(is_num(hyp)&&hyp>=0)
-    assert(is_num(adj)&&adj>=0)
+    assert(is_finite(hyp) && hyp>=0 && is_finite(adj) && adj>=0, 
+           "Triangle side lengths should be a positive numbers." )
     sqrt(hyp*hyp-adj*adj);
 
 
@@ -583,8 +621,8 @@ function hyp_adj_to_opp(hyp,adj) =
 // Example:
 //   opp = hyp_ang_to_opp(8,30);  // Returns: 4
 function hyp_ang_to_opp(hyp,ang) =
-    assert(is_num(hyp)&&hyp>=0)
-    assert(is_num(ang)&&ang>0&&ang<90)
+    assert(is_finite(hyp)&&hyp>=0, "Triangle side length should be a positive number." )
+    assert(is_finite(ang) && ang>0 && ang<90, "The angle should be an acute angle." )
     hyp*sin(ang);
 
 
@@ -599,8 +637,8 @@ function hyp_ang_to_opp(hyp,ang) =
 // Example:
 //   opp = adj_ang_to_opp(8,45);  // Returns: 8
 function adj_ang_to_opp(adj,ang) =
-    assert(is_num(adj)&&adj>=0)
-    assert(is_num(ang)&&ang>0&&ang<90)
+    assert(is_finite(adj)&&adj>=0, "Triangle side length should be a positive number." )
+    assert(is_finite(ang) && ang>0 && ang<90, "The angle should be an acute angle." )
     adj*tan(ang);
 
 
@@ -615,8 +653,8 @@ function adj_ang_to_opp(adj,ang) =
 // Example:
 //   hyp = adj_opp_to_hyp(3,4);  // Returns: 5
 function adj_opp_to_hyp(adj,opp) =
-    assert(is_num(adj)&&adj>=0)
-    assert(is_num(opp)&&opp>=0)
+    assert(is_finite(opp) && opp>=0 && is_finite(adj) && adj>=0, 
+           "Triangle side lengths should be a positive numbers." )
     norm([opp,adj]);
 
 
@@ -631,8 +669,8 @@ function adj_opp_to_hyp(adj,opp) =
 // Example:
 //   hyp = adj_ang_to_hyp(4,60);  // Returns: 8
 function adj_ang_to_hyp(adj,ang) =
-    assert(is_num(adj)&&adj>=0)
-    assert(is_num(ang)&&ang>=0&&ang<90)
+    assert(is_finite(adj) && adj>=0, "Triangle side length should be a positive number." )
+    assert(is_finite(ang) && ang>0 && ang<90, "The angle should be an acute angle." )
     adj/cos(ang);
 
 
@@ -647,8 +685,8 @@ function adj_ang_to_hyp(adj,ang) =
 // Example:
 //   hyp = opp_ang_to_hyp(4,30);  // Returns: 8
 function opp_ang_to_hyp(opp,ang) =
-    assert(is_num(opp)&&opp>=0)
-    assert(is_num(ang)&&ang>0&&ang<=90)
+    assert(is_finite(opp) && opp>=0, "Triangle side length should be a positive number." )
+    assert(is_finite(ang) && ang>0 && ang<90, "The angle should be an acute angle." )
     opp/sin(ang);
 
 
@@ -663,8 +701,8 @@ function opp_ang_to_hyp(opp,ang) =
 // Example:
 //   ang = hyp_adj_to_ang(8,4);  // Returns: 60 degrees
 function hyp_adj_to_ang(hyp,adj) =
-    assert(is_num(hyp)&&hyp>0)
-    assert(is_num(adj)&&adj>=0)
+    assert(is_finite(hyp) && hyp>0 && is_finite(adj) && adj>=0, 
+            "Triangle side lengths should be positive numbers." )
     acos(adj/hyp);
 
 
@@ -679,8 +717,8 @@ function hyp_adj_to_ang(hyp,adj) =
 // Example:
 //   ang = hyp_opp_to_ang(8,4);  // Returns: 30 degrees
 function hyp_opp_to_ang(hyp,opp) =
-    assert(is_num(hyp)&&hyp>0)
-    assert(is_num(opp)&&opp>=0)
+    assert(is_finite(hyp+opp) && hyp>0 && opp>=0, 
+            "Triangle side lengths should be positive numbers." )
     asin(opp/hyp);
 
 
@@ -695,8 +733,8 @@ function hyp_opp_to_ang(hyp,opp) =
 // Example:
 //   ang = adj_opp_to_ang(sqrt(3)/2,0.5);  // Returns: 30 degrees
 function adj_opp_to_ang(adj,opp) =
-    assert(is_num(adj)&&adj>=0)
-    assert(is_num(opp)&&opp>=0)
+    assert(is_finite(adj+opp) && adj>0 && opp>=0, 
+            "Triangle side lengths should be positive numbers." )
     atan2(opp,adj);
 
 
@@ -709,55 +747,63 @@ function adj_opp_to_ang(adj,opp) =
 // Examples:
 //   triangle_area([0,0], [5,10], [10,0]);  // Returns -50
 //   triangle_area([10,0], [5,10], [0,0]);  // Returns 50
-function triangle_area(a,b,c) =
-    len(a)==3? 0.5*norm(cross(c-a,c-b)) : (
-        a.x * (b.y - c.y) +
-        b.x * (c.y - a.y) +
-        c.x * (a.y - b.y)
-    ) / 2;
+function triangle_area(a,b,c) = 
+    assert( is_path([a,b,c]), "Invalid points or incompatible dimensions." )    
+    len(a)==3 
+    ? 0.5*norm(cross(c-a,c-b)) 
+    : 0.5*cross(c-a,c-b);
 
 
 
 // Section: Planes
 
+
 // Function: plane3pt()
 // Usage:
 //   plane3pt(p1, p2, p3);
 // Description:
-//   Generates the cartesian equation of a plane from three non-collinear points on the plane.
+//   Generates the cartesian equation of a plane from three 3d points.
 //   Returns [A,B,C,D] where Ax + By + Cz = D is the equation of a plane.
+//   Returns [], if the points are collinear.
 // Arguments:
 //   p1 = The first point on the plane.
 //   p2 = The second point on the plane.
 //   p3 = The third point on the plane.
 function plane3pt(p1, p2, p3) =
+    assert( is_path([p1,p2,p3],dim=3) && len(p1)==3,
+            "Invalid points or incompatible dimensions." )    
     let(
-        p1=point3d(p1),
-        p2=point3d(p2),
-        p3=point3d(p3),
-        normal = unit(cross(p3-p1, p2-p1))
-    ) concat(normal, [normal*p1]);
+        crx = cross(p3-p1, p2-p1),
+        nrm = norm(crx)
+    ) 
+    approx(nrm,0) ? [] :
+    concat(crx/nrm, [crx*p1]/nrm);
 
 
 // Function: plane3pt_indexed()
 // Usage:
 //   plane3pt_indexed(points, i1, i2, i3);
 // Description:
-//   Given a list of points, and the indices of three of those points,
+//   Given a list of 3d points, and the indices of three of those points,
 //   generates the cartesian equation of a plane that those points all
-//   lie on.  Requires that the three indexed points be non-collinear.
-//   Returns [A,B,C,D] where Ax+By+Cz=D is the equation of a plane.
+//   lie on. If the points are not collinear, returns [A,B,C,D] where Ax+By+Cz=D is the equation of a plane.
+//   If they are collinear, returns [].
 // Arguments:
 //   points = A list of points.
 //   i1 = The index into `points` of the first point on the plane.
 //   i2 = The index into `points` of the second point on the plane.
 //   i3 = The index into `points` of the third point on the plane.
 function plane3pt_indexed(points, i1, i2, i3) =
+    assert( is_vector([i1,i2,i3]) && min(i1,i2,i3)>=0 && is_list(points) && max(i1,i2,i3)<len(points),
+            "Invalid or out of range indices." )
+    assert( is_path([points[i1], points[i2], points[i3]],dim=3),
+            "Improper points or improper dimensions." )
     let(
         p1 = points[i1],
         p2 = points[i2],
         p3 = points[i3]
-    ) plane3pt(p1,p2,p3);
+    ) 
+    plane3pt(p1,p2,p3);
 
 
 // Function: plane_from_normal()
@@ -768,18 +814,19 @@ function plane3pt_indexed(points, i1, i2, i3) =
 // Example:
 //   plane_from_normal([0,0,1], [2,2,2]);  // Returns the xy plane passing through the point (2,2,2)
 function plane_from_normal(normal, pt=[0,0,0]) =
+  assert( is_matrix([normal,pt],2,3) && !approx(norm(normal),0), 
+          "Inputs `normal` and `pt` should 3d vectors/points and `normal` cannot be zero." )
   concat(normal, [normal*pt]);
 
 
 // Function: plane_from_points()
 // Usage:
-//   plane_from_points(points, [fast], [eps]);
+//   plane_from_points(points, <fast>, <eps>);
 // Description:
-//   Given a list of 3 or more coplanar points, returns the cartesian equation of a plane.
-//   Returns [A,B,C,D] where Ax+By+Cz=D is the equation of the plane.
-//   If not all the points in the points list are coplanar, then `undef` is returned.
-//   If `fast` is true, then a list where not all points are coplanar will result
-//   in an invalid plane value, as all coplanar checks are skipped.
+//   Given a list of 3 or more coplanar 3D points, returns the coefficients of the cartesian equation of a plane,
+//   that is [A,B,C,D] where Ax+By+Cz=D is the equation of the plane.
+//   If `fast` is false and the points in the list are collinear or not coplanar, then `undef` is returned.
+//   if `fast` is true, then the coplanarity test is skipped and a plane passing through 3 non-collinear arbitrary points is returned.
 // Arguments:
 //   points = The list of points to find the plane of.
 //   fast = If true, don't verify that all points in the list are coplanar.  Default: false
@@ -791,31 +838,33 @@ function plane_from_normal(normal, pt=[0,0,0]) =
 //   cp = centroid(xyzpath);
 //   move(cp) rot(from=UP,to=plane_normal(plane)) anchor_arrow();
 function plane_from_points(points, fast=false, eps=EPSILON) =
+    assert( is_path(points,dim=3), "Improper 3d point list." )
+    assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." )
     let(
         points = deduplicate(points),
-        indices = sort(find_noncollinear_points(points)),
+        indices = noncollinear_triple(points,error=false)
+    )
+    indices==[] ? undef :
+    let(
         p1 = points[indices[0]],
         p2 = points[indices[1]],
         p3 = points[indices[2]],
-        plane = plane3pt(p1,p2,p3),
-        all_coplanar = fast || all([
-            for (pt = points) coplanar(plane,pt,eps=eps)
-        ])
-    ) all_coplanar? plane : undef;
+        plane = plane3pt(p1,p2,p3)
+    )
+    fast || points_on_plane(points,plane,eps=eps) ? plane : undef;
 
 
 // Function: plane_from_polygon()
 // Usage:
 //   plane_from_polygon(points, [fast], [eps]);
 // Description:
-//   Given a 3D planar polygon, returns the cartesian equation of a plane.
+//   Given a 3D planar polygon, returns the cartesian equation of its plane.
 //   Returns [A,B,C,D] where Ax+By+Cz=D is the equation of the plane.
-//   If not all the points in the polygon are coplanar, then `undef` is returned.
-//   If `fast` is true, then a polygon where not all points are coplanar will
-//   result in an invalid plane value, as all coplanar checks are skipped.
+//   If not all the points in the polygon are coplanar, then [] is returned.
+//   If `fast` is true, the polygon coplanarity check is skipped and the plane may not contain all polygon points.
 // Arguments:
 //   poly = The planar 3D polygon to find the plane of.
-//   fast = If true, don't verify that all points in the polygon are coplanar.  Default: false
+//   fast = If true, doesn't verify that all points in the polygon are coplanar.  Default: false
 //   eps = How much variance is allowed in testing that each point is on the same plane.  Default: `EPSILON` (1e-9)
 // Example(3D):
 //   xyzpath = rot(45, v=[0,1,0], p=path3d(star(n=5,step=2,d=100), 70));
@@ -824,17 +873,14 @@ function plane_from_points(points, fast=false, eps=EPSILON) =
 //   cp = centroid(xyzpath);
 //   move(cp) rot(from=UP,to=plane_normal(plane)) anchor_arrow();
 function plane_from_polygon(poly, fast=false, eps=EPSILON) =
+    assert( is_path(poly,dim=3), "Invalid polygon." )
+    assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." )
     let(
         poly = deduplicate(poly),
         n = polygon_normal(poly),
         plane = [n.x, n.y, n.z, n*poly[0]]
-    ) fast? plane : let(
-        all_coplanar = [
-            for (pt = poly)
-            if (!coplanar(plane,pt,eps=eps)) 1
-        ] == []
-    ) all_coplanar? plane :
-    undef;
+    ) 
+    fast? plane: coplanar(poly,eps=eps)? plane: [];
 
 
 // Function: plane_normal()
@@ -842,7 +888,9 @@ function plane_from_polygon(poly, fast=false, eps=EPSILON) =
 //   plane_normal(plane);
 // Description:
 //   Returns the unit length normal vector for the given plane.
-function plane_normal(plane) = unit([for (i=[0:2]) plane[i]]);
+function plane_normal(plane) = 
+    assert( _valid_plane(plane), "Invalid input plane." )
+    unit([plane.x, plane.y, plane.z]);
 
 
 // Function: plane_offset()
@@ -851,7 +899,9 @@ function plane_normal(plane) = unit([for (i=[0:2]) plane[i]]);
 // Description:
 //   Returns D, or the scalar offset of the plane from the origin. This can be a negative value.
 //   The absolute value of this is the distance of the plane from the origin at its closest approach.
-function plane_offset(plane) = plane[3];
+function plane_offset(plane) =  
+    assert( _valid_plane(plane), "Invalid input plane." )
+    plane[3]/norm([plane.x, plane.y, plane.z]);
 
 
 // Function: plane_transform()
@@ -859,8 +909,8 @@ function plane_offset(plane) = plane[3];
 //   mat = plane_transform(plane);
 // Description:
 //   Given a plane definition `[A,B,C,D]`, where `Ax+By+Cz=D`, returns a 3D affine
-//   transformation matrix that will rotate and translate from points on that plane
-//   to points on the XY plane.  You can generally then use `path2d()` to drop the
+//   transformation matrix that will linear transform points on that plane
+//   into points on the XY plane.  You can generally then use `path2d()` to drop the
 //   Z coordinates, so you can work with the points in 2D.
 // Arguments:
 //   plane = The `[A,B,C,D]` plane definition where `Ax+By+Cz=D` is the formula of the plane.
@@ -875,7 +925,34 @@ function plane_transform(plane) =
     let(
         n = plane_normal(plane),
         cp = n * plane[3]
-    ) rot(from=n, to=UP) * move(-cp);
+        ) 
+    rot(from=n, to=UP) * move(-cp);
+
+
+// Function: projection_on_plane()
+// Usage:
+//   projection_on_plane(points);
+// Description:
+//   Given a plane definition `[A,B,C,D]`, where `Ax+By+Cz=D`, and a list of 2d or 3d points, return the 3D orthogonal 
+//   projection of the points on the plane.
+// Arguments:
+//   plane = The `[A,B,C,D]` plane definition where `Ax+By+Cz=D` is the formula of the plane.
+//   points = List of points to project
+// Example(3D):
+//   points = move([10,20,30], p=yrot(25, p=path3d(circle(d=100))));
+//   plane = plane3pt([1,0,0],[0,1,0],[0,0,1]);
+//   proj = projection_on_plane(plane,points);
+function projection_on_plane(plane, points) =
+    assert( _valid_plane(plane), "Invalid plane." )
+    assert( is_path(points), "Invalid list of points or dimension." )
+    let( 
+        p  = len(points[0])==2
+             ? [for(pi=points) point3d(pi) ]
+             : points, 
+        plane = plane/norm([plane.x,plane.y,plane.z]),
+        n = [plane.x,plane.y,plane.z]
+        ) 
+    [for(pi=p) pi - (pi*n - plane[3])*n];
 
 
 // Function: plane_point_nearest_origin()
@@ -899,9 +976,12 @@ function plane_point_nearest_origin(plane) =
 //   will be negative.  The normal of the plane is the same as [A,B,C].
 // Arguments:
 //   plane = The [A,B,C,D] values for the equation of the plane.
-//   point = The point to test.
+//   point = The distance evaluation point.
 function distance_from_plane(plane, point) =
-    [plane.x, plane.y, plane.z] * point3d(point) - plane[3];
+    assert( _valid_plane(plane), "Invalid input plane." )
+    assert( is_vector(point,3), "The point should be a 3D point." )
+    let( nrml = [plane.x, plane.y, plane.z] )
+    ( nrml* point - plane[3])/norm(nrml);
 
 
 // Function: closest_point_on_plane()
@@ -911,34 +991,27 @@ function distance_from_plane(plane, point) =
 //   Takes a point, and a plane [A,B,C,D] where the equation of that plane is `Ax+By+Cz=D`.
 //   Returns the coordinates of the closest point on that plane to the given `point`.
 // Arguments:
-//   plane = The [A,B,C,D] values for the equation of the plane.
+//   plane = The [A,B,C,D] coefficients for the equation of the plane.
 //   point = The 3D point to find the closest point to.
 function closest_point_on_plane(plane, point) =
+    assert( _valid_plane(plane), "Invalid input plane." )
+    assert( is_vector(point,3), "Invalid point." )
     let(
-        n = unit(plane_normal(plane)),
+        n = unit([plane.x, plane.y, plane.z]),
         d = distance_from_plane(plane, point)
-    ) point - n*d;
+        ) 
+    point - n*d;
 
 
 // Returns [POINT, U] if line intersects plane at one point.
 // Returns [LINE, undef] if the line is on the plane.
 // Returns undef if line is parallel to, but not on the given plane.
 function _general_plane_line_intersection(plane, line, eps=EPSILON) =
-    let(
-        p0 = line[0],
-        p1 = line[1],
-        n = plane_normal(plane),
-        u = p1 - p0,
-        d = n * u
-    ) abs(d)<eps? (
-        coplanar(plane, p0)? [line,undef] :  // Line on plane
-        undef  // Line parallel to plane
-    ) : let(
-        v0 = closest_point_on_plane(plane, [0,0,0]),
-        w = p0 - v0,
-        s1 = (-n * w) / d,
-        pt = s1 * u + p0
-    ) [pt, s1];
+    let( a = plane*[each line[0],-1],
+         b = plane*[each(line[1]-line[0]),-1] )
+    approx(b,0,eps) 
+    ? points_on_plane(line[0],plane,eps)? [line,undef]: undef
+    : [ line[0]+a/b*(line[1]-line[0]), a/b ];
 
 
 // Function: plane_line_angle()
@@ -948,38 +1021,41 @@ function _general_plane_line_intersection(plane, line, eps=EPSILON) =
 //   The resulting angle is signed, with the sign positive if the vector p2-p1 lies on 
 //   the same side of the plane as the plane's normal vector.  
 function plane_line_angle(plane, line) =
+    assert( _valid_plane(plane), "Invalid plane." )
+    assert( _valid_line(line), "Invalid line." )
     let(
         vect = line[1]-line[0],
         zplane = plane_normal(plane),
         sin_angle = vect*zplane/norm(zplane)/norm(vect)
-    ) asin(constrain(sin_angle,-1,1));
+        ) 
+    asin(constrain(sin_angle,-1,1));
 
 
 // Function: plane_line_intersection()
 // Usage:
-//   pt = plane_line_intersection(plane, line, [eps]);
+//   pt = plane_line_intersection(plane, line, [bounded], [eps]);
 // Description:
 //   Takes a line, and a plane [A,B,C,D] where the equation of that plane is `Ax+By+Cz=D`.
 //   If `line` intersects `plane` at one point, then that intersection point is returned.
 //   If `line` lies on `plane`, then the original given `line` is returned.
-//   If `line` is parallel to, but not on `plane`, then `undef` is returned.
+//   If `line` is parallel to, but not on `plane`, then undef is returned.
 // Arguments:
 //   plane = The [A,B,C,D] values for the equation of the plane.
-//   line = A list of two 3D points that are on the line.
+//   line = A list of two distinct 3D points that are on the line.
 //   bounded = If false, the line is considered unbounded.  If true, it is treated as a bounded line segment.  If given as `[true, false]` or `[false, true]`, the boundedness of the points are specified individually, allowing the line to be treated as a half-bounded ray.  Default: false (unbounded)
-//   eps = The epsilon error value to determine whether the line is too close to parallel to the plane.  Default: `EPSILON` (1e-9)
+//   eps = The tolerance value in determining whether the line is parallel to the plane.  Default: `EPSILON` (1e-9)
 function plane_line_intersection(plane, line, bounded=false, eps=EPSILON) =
-    assert(is_vector(plane)&&len(plane)==4, "Invalid plane value.")
-    assert(is_path(line)&&len(line)==2, "Invalid line value.")
-    assert(!approx(line[0],line[1]), "The two points defining the line must not be the same point.")
+    assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." )
+    assert(_valid_plane(plane,eps=eps) && _valid_line(line,dim=3,eps=eps), "Invalid plane and/or line.")
+    assert(is_bool(bounded) || (is_list(bounded) && len(bounded)==2), "Invalid bound condition(s).")
     let(
         bounded = is_list(bounded)? bounded : [bounded, bounded],
         res = _general_plane_line_intersection(plane, line, eps=eps)
     )
-    is_undef(res)? undef :
-    is_undef(res[1])? res[0] :
-    bounded[0]&&res[1]<0? undef :
-    bounded[1]&&res[1]>1? undef :
+    is_undef(res) ? undef :
+    is_undef(res[1]) ? res[0] :
+    bounded[0] && res[1]<0 ? undef :
+    bounded[1] && res[1]>1 ? undef :
     res[0];
 
 
@@ -988,22 +1064,28 @@ function plane_line_intersection(plane, line, bounded=false, eps=EPSILON) =
 //   pt = polygon_line_intersection(poly, line, [bounded], [eps]);
 // Description:
 //   Takes a possibly bounded line, and a 3D planar polygon, and finds their intersection point.
-//   If the line is on the plane as the polygon, and intersects, then a list of 3D line
-//   segments is returned, one for each section of the line that is inside the polygon.
-//   If the line is not on the plane of the polygon, but intersects, then the 3D intersection
-//   point is returned.  If the line does not intersect the polygon, then `undef` is returned.
+//   If the line and the polygon are on the same plane then returns a list, possibly empty, of 3D line
+//   segments, one for each section of the line that is inside the polygon.
+//   If the line is not on the plane of the polygon, but intersects it, then returns the 3D intersection
+//   point.  If the line does not intersect the polygon, then `undef` is returned.
 // Arguments:
 //   poly = The 3D planar polygon to find the intersection with.
-//   line = A list of two 3D points that are on the line.
+//   line = A list of two distinct 3D points on the line.
 //   bounded = If false, the line is considered unbounded.  If true, it is treated as a bounded line segment.  If given as `[true, false]` or `[false, true]`, the boundedness of the points are specified individually, allowing the line to be treated as a half-bounded ray.  Default: false (unbounded)
-//   eps = The epsilon error value to determine whether the line is too close to parallel to the plane.  Default: `EPSILON` (1e-9)
+//   eps = The tolerance value in determining whether the line is parallel to the plane.  Default: `EPSILON` (1e-9)
 function polygon_line_intersection(poly, line, bounded=false, eps=EPSILON) =
-    assert(is_path(poly))
-    assert(is_path(line)&&len(line)==2)
+    assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." )
+    assert(is_path(poly,dim=3), "Invalid polygon." )
+    assert(is_bool(bounded) || (is_list(bounded) && len(bounded)==2), "Invalid bound condition(s).")
+    assert(_valid_line(line,dim=3,eps=eps), "Invalid line." )
     let(
         bounded = is_list(bounded)? bounded : [bounded, bounded],
         poly = deduplicate(poly),
-        indices = sort(find_noncollinear_points(poly)),
+        indices = noncollinear_triple(poly)
+    )
+    indices==[] ? undef :
+    let(
+        indices = sort(indices),
         p1 = poly[indices[0]],
         p2 = poly[indices[1]],
         p3 = poly[indices[2]],
@@ -1011,34 +1093,31 @@ function polygon_line_intersection(poly, line, bounded=false, eps=EPSILON) =
         res = _general_plane_line_intersection(plane, line, eps=eps)
     )
     is_undef(res)? undef :
-    is_undef(res[1])? (
-        let(
-            // Line is on polygon plane.
+    is_undef(res[1])
+    ? ( let(// Line is on polygon plane.
             linevec = unit(line[1] - line[0]),
             lp1 = line[0] + (bounded[0]? 0 : -1000000) * linevec,
             lp2 = line[1] + (bounded[1]? 0 :  1000000) * linevec,
             poly2d = clockwise_polygon(project_plane(poly, p1, p2, p3)),
             line2d = project_plane([lp1,lp2], p1, p2, p3),
             parts = split_path_at_region_crossings(line2d, [poly2d], closed=false),
-            inside = [
-                for (part = parts)
-                if (point_in_polygon(mean(part), poly2d)>0) part
-            ]
-        ) !inside? undef :
+            inside = [for (part = parts)
+                          if (point_in_polygon(mean(part), poly2d)>0) part
+                     ]
+        ) 
+        !inside? undef :
         let(
-            isegs = [
-                for (seg = inside)
-                lift_plane(seg, p1, p2, p3)
-            ]
-        ) isegs
-    ) :
-    bounded[0]&&res[1]<0? undef :
-    bounded[1]&&res[1]>1? undef :
-    let(
-        proj = clockwise_polygon(project_plane(poly, p1, p2, p3)),
-        pt = project_plane(res[0], p1, p2, p3)
-    ) point_in_polygon(pt, proj) < 0? undef :
-    res[0];
+            isegs = [for (seg = inside) lift_plane(seg, p1, p2, p3) ]
+        ) 
+        isegs
+    ) 
+    :   bounded[0]&&res[1]<0? [] :
+        bounded[1]&&res[1]>1? [] :
+        let(
+            proj = clockwise_polygon(project_plane(poly, p1, p2, p3)),
+            pt = project_plane(res[0], p1, p2, p3)
+        ) 
+        point_in_polygon(pt, proj) < 0 ? undef : res[0];
 
 
 // Function: plane_intersection()
@@ -1047,53 +1126,62 @@ function polygon_line_intersection(poly, line, bounded=false, eps=EPSILON) =
 // Description:
 //   Compute the point which is the intersection of the three planes, or the line intersection of two planes.
 //   If you give three planes the intersection is returned as a point.  If you give two planes the intersection
-//   is returned as a list of two points on the line of intersection.  If any of the input planes are parallel
-//   then returns undef.  
+//   is returned as a list of two points on the line of intersection.  If any two input planes are parallel
+//   or coincident then returns undef.  
 function plane_intersection(plane1,plane2,plane3) =
-    is_def(plane3)? let(
-        matrix = [for(p=[plane1,plane2,plane3]) select(p,0,2)],
-        rhs = [for(p=[plane1,plane2,plane3]) p[3]]
-    ) linear_solve(matrix,rhs) :
-    let(
-        normal = cross(plane_normal(plane1), plane_normal(plane2))
-    ) approx(norm(normal),0) ? undef :
-    let(
-        matrix = [for(p=[plane1,plane2]) select(p,0,2)],
-        rhs = [for(p=[plane1,plane2]) p[3]],
-        point = linear_solve(matrix,rhs)
-    ) is_undef(point)? undef :
-    [point, point+normal];
+    assert( _valid_plane(plane1) && _valid_plane(plane2) && (is_undef(plane3) ||_valid_plane(plane3)),
+                "The input must be 2 or 3 planes." )
+    is_def(plane3)
+    ?   let(
+          matrix = [for(p=[plane1,plane2,plane3]) select(p,0,2)],
+          rhs = [for(p=[plane1,plane2,plane3]) p[3]]
+        ) 
+        linear_solve(matrix,rhs)
+    :   let( normal = cross(plane_normal(plane1), plane_normal(plane2)) ) 
+        approx(norm(normal),0) ? undef :
+        let(
+            matrix = [for(p=[plane1,plane2]) select(p,0,2)],
+            rhs = [for(p=[plane1,plane2]) p[3]],
+            point = linear_solve(matrix,rhs)
+        ) 
+        point==[]? undef: [point, point+normal];
 
 
 // Function: coplanar()
 // Usage:
-//   coplanar(plane, point);
+//   coplanar(points,<eps>);
 // Description:
-//   Given a plane as [A,B,C,D] where the cartesian equation for that plane
-//   is Ax+By+Cz=D, determines if the given point is on that plane.
-//   Returns true if the point is on that plane.
+//   Returns true if the given 3D points are non-collinear and are on a plane.
 // Arguments:
-//   plane = The [A,B,C,D] values for the equation of the plane.
-//   point = The point to test.
-//   eps = How much variance is allowed in testing that each point is on the same plane.  Default: `EPSILON` (1e-9)
-function coplanar(plane, point, eps=EPSILON) =
-    abs(distance_from_plane(plane, point)) <= eps;
+//   points = The points to test.
+//   eps = How much variance is allowed in the planarity test.  Default: `EPSILON` (1e-9)
+function coplanar(points, eps=EPSILON) =
+    assert( is_path(points,dim=3) , "Input should be a list of 3D points." )
+    assert( is_finite(eps) && eps>=0, "The tolerance should be a non-negative number." )
+    len(points)<=2 ? false
+    :   let( ip = noncollinear_triple(points,error=false,eps=eps) ) 
+        ip == [] ? false : 
+        let( plane  = plane3pt(points[ip[0]],points[ip[1]],points[ip[2]]), 
+             normal = point3d(plane) )
+        max( points*normal ) - plane[3]< eps*norm(normal);
 
-
-// Function: points_are_coplanar()
+    
+// Function: points_on_plane()
 // Usage:
-//   points_are_coplanar(points, [eps]);
+//   points_on_plane(points, plane, <eps>);
 // Description:
-//   Given a list of points, returns true if all points in the list are coplanar.
+//   Returns true if the given 3D points are on the given plane.
 // Arguments:
-//   points = The list of points to test.
-//   eps = How much variance is allowed in testing that each point is on the same plane.  Default: `EPSILON` (1e-9)
-function points_are_coplanar(points, eps=EPSILON) =
-    points_are_collinear(points, eps=eps)? true :
-    let(
-        plane = plane_from_points(points, fast=true, eps=eps)
-    ) all([for (pt = points) coplanar(plane, pt, eps=eps)]);
-
+//   plane = The plane to test the points on.
+//   points = The list of 3D points to test.
+//   eps = How much variance is allowed in the planarity testing.  Default: `EPSILON` (1e-9)
+function points_on_plane(points, plane, eps=EPSILON) =
+    assert( _valid_plane(plane), "Invalid plane." )
+    assert( is_matrix(points,undef,3) && len(points)>0, "Invalid pointlist." ) // using is_matrix it accepts len(points)==1
+    assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." )
+    let( normal = point3d(plane),
+         pt_nrm = points*normal )
+    abs(max( max(pt_nrm) - plane[3], -min(pt_nrm)+plane[3]))< eps*norm(normal);
 
 
 // Function: in_front_of_plane()
@@ -1101,12 +1189,12 @@ function points_are_coplanar(points, eps=EPSILON) =
 //   in_front_of_plane(plane, point);
 // Description:
 //   Given a plane as [A,B,C,D] where the cartesian equation for that plane
-//   is Ax+By+Cz=D, determines if the given point is on the side of that
+//   is Ax+By+Cz=D, determines if the given 3D point is on the side of that
 //   plane that the normal points towards.  The normal of the plane is the
 //   same as [A,B,C].
 // Arguments:
-//   plane = The [A,B,C,D] values for the equation of the plane.
-//   point = The point to test.
+//   plane = The [A,B,C,D] coefficients for the equation of the plane.
+//   point = The 3D point to test.
 function in_front_of_plane(plane, point) =
     distance_from_plane(plane, point) > EPSILON;
 
@@ -1116,7 +1204,7 @@ function in_front_of_plane(plane, point) =
 
 // Function: find_circle_2tangents()
 // Usage:
-//   find_circle_2tangents(pt1, pt2, pt3, r|d, [tangents]);
+//   find_circle_2tangents(pt1, pt2, pt3, r|d, <tangents>);
 // Description:
 //   Given a pair of rays with a common origin, and a known circle radius/diameter, finds
 //   the centerpoint for the circle of that size that touches both rays tangentally.
@@ -1156,37 +1244,46 @@ function in_front_of_plane(plane, point) =
 function find_circle_2tangents(pt1, pt2, pt3, r, d, tangents=false) =
     let(r = get_radius(r=r, d=d, dflt=undef))
     assert(r!=undef, "Must specify either r or d.")
-    (is_undef(pt2) && is_undef(pt3) && is_list(pt1))? find_circle_2tangents(pt1[0], pt1[1], pt1[2], r=r) :
-    collinear(pt1, pt2, pt3)? undef :
-    let(
-        v1 = unit(pt1 - pt2),
-        v2 = unit(pt3 - pt2),
-        vmid = unit(mean([v1, v2])),
-        n = vector_axis(v1, v2),
-        a = vector_angle(v1, v2),
-        hyp = r / sin(a/2),
-        cp = pt2 + hyp * vmid
-    ) !tangents? [cp, n] :
-    let(
-        x = hyp * cos(a/2),
-        tp1 = pt2 + x * v1,
-        tp2 = pt2 + x * v2,
-        fff=echo(tp1=tp1,cp=cp,pt2=pt2),
-        dang1 = vector_angle(tp1-cp,pt2-cp),
-        dang2 = vector_angle(tp2-cp,pt2-cp)
-    ) [cp, n, tp1, tp2, dang1, dang2];
+    assert( ( is_path(pt1) && len(pt1)==3 && is_undef(pt2) && is_undef(pt3)) 
+            || (is_matrix([pt1,pt2,pt3]) && (len(pt1)==2 || len(pt1)==3) ),
+            "Invalid input points." )
+    is_undef(pt2) 
+    ? find_circle_2tangents(pt1[0], pt1[1], pt1[2], r=r, tangents=tangents) 
+    : collinear(pt1, pt2, pt3)? undef :
+        let(
+            v1 = unit(pt1 - pt2),
+            v2 = unit(pt3 - pt2),
+            vmid = unit(mean([v1, v2])),
+            n = vector_axis(v1, v2),
+            a = vector_angle(v1, v2),
+            hyp = r / sin(a/2),
+            cp = pt2 + hyp * vmid
+            ) 
+        !tangents ? [cp, n] :
+        let(
+            x = hyp * cos(a/2),
+            tp1 = pt2 + x * v1,
+            tp2 = pt2 + x * v2,
+//            fff=echo(tp1=tp1,cp=cp,pt2=pt2),
+            dang1 = vector_angle(tp1-cp,pt2-cp),
+            dang2 = vector_angle(tp2-cp,pt2-cp)
+            ) 
+        [cp, n, tp1, tp2, dang1, dang2];
 
 
 // Function: find_circle_3points()
 // Usage:
 //   find_circle_3points(pt1, pt2, pt3);
+//   find_circle_3points([pt1, pt2, pt3]);
 // Description:
 //   Returns the [CENTERPOINT, RADIUS, NORMAL] of the circle that passes through three non-collinear
-//   points.  The centerpoint will be a 2D or 3D vector, depending on the points input.  If all three
+//   points where NORMAL is the normal vector of the plane that the circle is on (UP or DOWN if the points are 2D).
+//   The centerpoint will be a 2D or 3D vector, depending on the points input.  If all three
 //   points are 2D, then the resulting centerpoint will be 2D, and the normal will be UP ([0,0,1]).
 //   If any of the points are 3D, then the resulting centerpoint will be 3D.  If the three points are
 //   collinear, then `[undef,undef,undef]` will be returned.  The normal will be a normalized 3D
 //   vector with a non-negative Z axis.
+//   Instead of 3 arguments, it is acceptable to input the 3 points in a list `pt1`, leaving `pt2`and `pt3` as undef.
 // Arguments:
 //   pt1 = The first point.
 //   pt2 = The second point.
@@ -1198,49 +1295,43 @@ function find_circle_2tangents(pt1, pt2, pt3, r, d, tangents=false) =
 //   translate(circ[0]) color("red") circle(d=3, $fn=12);
 //   move_copies(pts) color("blue") circle(d=3, $fn=12);
 function find_circle_3points(pt1, pt2, pt3) =
-    (is_undef(pt2) && is_undef(pt3) && is_list(pt1))? find_circle_3points(pt1[0], pt1[1], pt1[2]) :
-    collinear(pt1,pt2,pt3)? [undef,undef,undef] :
-    let(
-        v1 = pt1-pt2,
-        v2 = pt3-pt2,
-        n = vector_axis(v1,v2),
-        n2 = n.z<0? -n : n
-    ) len(pt1)+len(pt2)+len(pt3)>6? (
+    (is_undef(pt2) && is_undef(pt3) && is_list(pt1))
+    ? find_circle_3points(pt1[0], pt1[1], pt1[2]) 
+    :   assert( is_vector(pt1) && is_vector(pt2) && is_vector(pt3) 
+                && max(len(pt1),len(pt2),len(pt3))<=3 && min(len(pt1),len(pt2),len(pt3))>=2,
+                "Invalid point(s)." )
+        collinear(pt1,pt2,pt3)? [undef,undef,undef] :
         let(
-            a = project_plane(pt1, pt1, pt2, pt3),
-            b = project_plane(pt2, pt1, pt2, pt3),
-            c = project_plane(pt3, pt1, pt2, pt3),
-            res = find_circle_3points(a, b, c)
-        ) res[0]==undef? [undef,undef,undef] : let(
-            cp = lift_plane(res[0], pt1, pt2, pt3),
-            r = norm(pt2-cp)
-        ) [cp, r, n2]
-    ) : let(
-        mp1 = pt2 + v1/2,
-        mp2 = pt2 + v2/2,
-        mpv1 = rot(90, v=n, p=v1),
-        mpv2 = rot(90, v=n, p=v2),
-        l1 = [mp1, mp1+mpv1],
-        l2 = [mp2, mp2+mpv2],
-        isect = line_intersection(l1,l2)
-    ) is_undef(isect)? [undef,undef,undef] : let(
-        r = norm(pt2-isect)
-    ) [isect, r, n2];
-
-
+            v  = [ point3d(pt1), point3d(pt2), point3d(pt3) ], // triangle vertices
+            ed = [for(i=[0:2]) v[(i+1)%3]-v[i] ],    // triangle edge vectors
+            pm = [for(i=[0:2]) v[(i+1)%3]+v[i] ]/2,  // edge mean points
+            es = sortidx( [for(di=ed) norm(di) ] ),   
+            e1 = ed[es[1]],                          // take the 2 longest edges
+            e2 = ed[es[2]],
+            n0 = vector_axis(e1,e2),                 // normal standardization 
+            n  = n0.z<0? -n0 : n0,
+            sc = plane_intersection(                 
+                    [ each e1, e1*pm[es[1]] ],       // planes orthogonal to 2 edges
+                    [ each e2, e2*pm[es[2]] ],
+                    [ each n,  n*v[0] ] ) ,          // triangle plane
+            cp = len(pt1)+len(pt2)+len(pt3)>6 ? sc: [sc.x, sc.y], 
+            r  = norm(sc-v[0])
+            )
+        [ cp, r, n ];
+     
 
 // Function: circle_point_tangents()
 // Usage:
 //   tangents = circle_point_tangents(r|d, cp, pt);
 // Description:
-//   Given a circle and a point outside that circle, finds the tangent point(s) on the circle for a
+//   Given a 2d circle and a 2d point outside that circle, finds the 2d tangent point(s) on the circle for a
 //   line passing through the point.  Returns list of zero or more sublists of [ANG, TANGPT]
 // Arguments:
 //   r = Radius of the circle.
 //   d = Diameter of the circle.
-//   cp = The coordinates of the circle centerpoint.
-//   pt = The coordinates of the external point.
-// Example(2D):
+//   cp = The coordinates of the 2d circle centerpoint.
+//   pt = The coordinates of the 2d external point.
+// Example:
 //   cp = [-10,-10];  r = 30;  pt = [30,10];
 //   tanpts = subindex(circle_point_tangents(r=r, cp=cp, pt=pt),1);
 //   color("yellow") translate(cp) circle(r=r);
@@ -1248,9 +1339,8 @@ function find_circle_3points(pt1, pt2, pt3) =
 //   color("red") move_copies(tanpts) circle(d=3,$fn=12);
 //   color("blue") move_copies([cp,pt]) circle(d=3,$fn=12);
 function circle_point_tangents(r, d, cp, pt) =
-    assert(is_num(r) || is_num(d))
-    assert(is_vector(cp))
-    assert(is_vector(pt))
+    assert(is_finite(r) || is_finite(d), "Invalid radius or diameter." )
+    assert(is_path([cp, pt],dim=2), "Invalid center point or external point.")
     let(
         r = get_radius(r=r, d=d, dflt=1),
         delta = pt - cp,
@@ -1264,11 +1354,10 @@ function circle_point_tangents(r, d, cp, pt) =
     ) [for (ang=angs) [ang, cp + r*[cos(ang),sin(ang)]]];
 
 
-
 // Function: circle_circle_tangents()
 // Usage: circle_circle_tangents(c1, r1|d1, c2, r2|d2)
 // Description:
-//   Computes lines tangents to a pair of circles.  Returns a list of line endpoints [p1,p2] where
+//   Computes 2d lines tangents to a pair of circles in 2d.  Returns a list of line endpoints [p1,p2] where
 //   p2 is the tangent point on circle 1 and p2 is the tangent point on circle 2.
 //   If four tangents exist then the first one the left hand exterior tangent as regarded looking from
 //   circle 1 toward circle 2.  The second value is the right hand exterior tangent.  The third entry
@@ -1314,6 +1403,7 @@ function circle_point_tangents(r, d, cp, pt) =
 //   move(c2) stroke(circle(r=r2), width=.1, closed=true);
 //   echo(pts);   // Returns []
 function circle_circle_tangents(c1,r1,c2,r2,d1,d2) =
+    assert( is_path([c1,c2],dim=2), "Invalid center point(s)." )
     let(
         r1 = get_radius(r1=r1,d1=d1),
         r2 = get_radius(r1=r2,d1=d2),
@@ -1342,24 +1432,32 @@ function circle_circle_tangents(c1,r1,c2,r2,d1,d2) =
 // Section: Pointlists
 
 
-// Function: find_noncollinear_points()
+// Function: noncollinear_triple()
 // Usage:
-//   find_noncollinear_points(points);
+//   noncollinear_triple(points);
 // Description:
 //   Finds the indices of three good non-collinear points from the points list `points`.
-function find_noncollinear_points(points,error=true,eps=EPSILON) =
+//   If all points are collinear, returns [].
+function noncollinear_triple(points,error=true,eps=EPSILON) =
+    assert( is_path(points), "Invalid input points." )
+    assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative number." )
     let(
         pa = points[0],
-        b = furthest_point(pa, points),
-        n = unit(points[b]-pa),
-        relpoints = [for(pt=points) pt-pa],
-        proj = relpoints * n,
-        distlist = [for(i=[0:len(points)-1]) norm(relpoints[i]-proj[i]*n)]
-    )
-    max(distlist)<eps
+        b  = furthest_point(pa, points),
+        pb = points[b],
+        nrm = norm(pa-pb)
+        )
+    approx(nrm, 0)
     ? assert(!error, "Cannot find three noncollinear points in pointlist.")
-      []
-    : [0,b,max_index(distlist)];
+        []
+    :   let(
+            n = (pb-pa)/nrm,
+            distlist = [for(i=[0:len(points)-1]) _dist2line(points[i]-pa, n)]   
+           )
+        max(distlist)<eps
+        ?  assert(!error, "Cannot find three noncollinear points in pointlist.")
+           []
+        :  [0,b,max_index(distlist)];    
 
 
 // Function: pointlist_bounds()
@@ -1372,13 +1470,14 @@ function find_noncollinear_points(points,error=true,eps=EPSILON) =
 // Arguments:
 //   pts = List of points.
 function pointlist_bounds(pts) =
-    assert(is_matrix(pts))
+    assert(is_matrix(pts) && len(pts)>0 && len(pts[0])>0 , "Invalid pointlist." ) 
     let(ptsT = transpose(pts))
     [
       [for(row=ptsT) min(row)],
       [for(row=ptsT) max(row)]
     ];
 
+
 // Function: closest_point()
 // Usage:
 //   closest_point(pt, points);
@@ -1388,6 +1487,8 @@ function pointlist_bounds(pts) =
 //   pt = The point to find the closest point to.
 //   points = The list of points to search.
 function closest_point(pt, points) =
+    assert( is_vector(pt), "Invalid point." )
+    assert(is_path(points,dim=len(pt)), "Invalid pointlist or incompatible dimensions." )
     min_index([for (p=points) norm(p-pt)]);
 
 
@@ -1400,6 +1501,8 @@ function closest_point(pt, points) =
 //   pt = The point to find the farthest point from.
 //   points = The list of points to search.
 function furthest_point(pt, points) =
+    assert( is_vector(pt), "Invalid point." )
+    assert(is_path(points,dim=len(pt)), "Invalid pointlist or incompatible dimensions." )
     max_index([for (p=points) norm(p-pt)]);
 
 
@@ -1410,37 +1513,48 @@ function furthest_point(pt, points) =
 // Usage:
 //   area = polygon_area(poly);
 // Description:
-//   Given a 2D or 3D planar polygon, returns the area of that polygon.  If the polygon is self-crossing, the results are undefined.
-function polygon_area(poly) =
-    len(poly)<3? 0 :
-    len(poly[0])==2? 0.5*sum([for(i=[0:1:len(poly)-1]) det2(select(poly,i,i+1))]) :
-    let(
-        plane = plane_from_points(poly)
-    ) plane==undef? undef :
-    let(
-        n = unit(plane_normal(plane)),
-        total = sum([for (i=[0:1:len(poly)-1]) cross(poly[i], select(poly,i+1))]),
-        res = abs(total * n) / 2
-    ) res;
+//   Given a 2D or 3D planar polygon, returns the area of that polygon.  
+//   If the polygon is self-crossing, the results are undefined. For non-planar points the result is undef.
+//   When `signed` is true, a signed area is returned; a positive area indicates a counterclockwise polygon.
+// Arguments:
+//   poly = polygon to compute the area of.
+//   signed = if true, a signed area is returned (default: false)
+function polygon_area(poly, signed=false) =
+    assert(is_path(poly), "Invalid polygon." )
+    len(poly)<3 ? 0 :
+    len(poly[0])==2
+    ? sum([for(i=[1:1:len(poly)-2]) cross(poly[i]-poly[0],poly[i+1]-poly[0]) ])/2
+    : let( plane = plane_from_points(poly) )
+      plane==undef? undef :
+      let( n = unit(plane_normal(plane)), 
+           total = sum([for(i=[1:1:len(poly)-1]) cross(poly[i]-poly[0],poly[i+1]-poly[0])*n ])/2
+          )
+      signed ? total : abs(total);
 
 
-// Function: polygon_is_convex()
+// Function: is_convex_polygon()
 // Usage:
-//   polygon_is_convex(poly);
+//   is_convex_polygon(poly);
 // Description:
-//   Returns true if the given polygon is convex.  Result is undefined if the polygon is self-intersecting.
+//   Returns true if the given 2D polygon is convex.  The result is meaningless if the polygon is not simple (self-intersecting).
+//   If the points are collinear the result is true. 
 // Example:
-//   polygon_is_convex(circle(d=50));  // Returns: true
+//   is_convex_polygon(circle(d=50));  // Returns: true
 // Example:
 //   spiral = [for (i=[0:36]) let(a=-i*10) (10+i)*[cos(a),sin(a)]];
-//   polygon_is_convex(spiral);  // Returns: false
-function polygon_is_convex(poly) =
-    let(
-        l = len(poly),
-        c = [for (i=idx(poly)) cross(poly[(i+1)%l]-poly[i],poly[(i+2)%l]-poly[(i+1)%l])]
-    )
-    len([for (x=c) if(x>0) 1])==0 ||
-    len([for (x=c) if(x<0) 1])==0;
+//   is_convex_polygon(spiral);  // Returns: false
+function is_convex_polygon(poly) =
+    assert(is_path(poly,dim=2), "The input should be a 2D polygon." )
+    let( l = len(poly) )
+    len([for( i = l-1, 
+              c = cross(poly[(i+1)%l]-poly[i], poly[(i+2)%l]-poly[(i+1)%l]), 
+              s = sign(c);
+            i>=0 && sign(c)==s;
+              i = i-1, 
+              c = i<0? 0: cross(poly[(i+1)%l]-poly[i],poly[(i+2)%l]-poly[(i+1)%l]), 
+              s = s==0 ? sign(c) : s
+            ) i 
+        ])== l;
 
 
 // Function: polygon_shift()
@@ -1454,6 +1568,7 @@ function polygon_is_convex(poly) =
 // Example:
 //   polygon_shift([[3,4], [8,2], [0,2], [-4,0]], 2);   // Returns [[0,2], [-4,0], [3,4], [8,2]]
 function polygon_shift(poly, i) =
+    assert(is_path(poly), "Invalid polygon." )
     list_rotate(cleanup_path(poly), i);
 
 
@@ -1461,13 +1576,15 @@ function polygon_shift(poly, i) =
 // Usage:
 //   polygon_shift_to_closest_point(path, pt);
 // Description:
-//   Given a polygon `path`, rotates the point ordering so that the first point in the path is the one closest to the given point `pt`.
-function polygon_shift_to_closest_point(path, pt) =
+//   Given a polygon `poly`, rotates the point ordering so that the first point in the path is the one closest to the given point `pt`.
+function polygon_shift_to_closest_point(poly, pt) =
+    assert(is_vector(pt), "Invalid point." )
+    assert(is_path(poly,dim=len(pt)), "Invalid polygon or incompatible dimension with the point." )
     let(
-        path = cleanup_path(path),
-        dists = [for (p=path) norm(p-pt)],
+        poly = cleanup_path(poly),
+        dists = [for (p=poly) norm(p-pt)],
         closest = min_index(dists)
-    ) select(path,closest,closest+len(path)-1);
+    ) select(poly,closest,closest+len(poly)-1);
 
 
 // Function: reindex_polygon()
@@ -1500,38 +1617,32 @@ function polygon_shift_to_closest_point(path, pt) =
 //   color("red") move_copies([pent[0],circ[0]]) circle(r=.1,$fn=32);
 //   color("blue") translate(reindexed[0])circle(r=.1,$fn=32);
 function reindex_polygon(reference, poly, return_error=false) = 
-    assert(is_path(reference) && is_path(poly))
-    assert(len(reference)==len(poly), "Polygons must be the same length in reindex_polygon")
+    assert(is_path(reference) && is_path(poly,dim=len(reference[0])),
+           "Invalid polygon(s) or incompatible dimensions. " )
+    assert(len(reference)==len(poly), "The polygons must have the same length.")
     let(
         dim = len(reference[0]),
         N = len(reference),
         fixpoly = dim != 2? poly :
-            polygon_is_clockwise(reference)? clockwise_polygon(poly) :
-            ccw_polygon(poly),
-        dist = [
-            // Matrix of all pairwise distances
-            for (p1=reference) [
-                for (p2=fixpoly) norm(p1-p2)
-            ]
-        ],
-        // Compute the sum of all distance pairs for a each shift
-        sums = [
-            for(shift=[0:1:N-1]) sum([
-                for(i=[0:1:N-1]) dist[i][(i+shift)%N]
-            ])
-        ],
-        optimal_poly = polygon_shift(fixpoly,min_index(sums))
-    )
-    return_error? [optimal_poly, min(sums)] :
+                  polygon_is_clockwise(reference)
+                  ? clockwise_polygon(poly)
+                  : ccw_polygon(poly),
+        I   = [for(i=[0:N-1]) 1],
+        val = [ for(k=[0:N-1]) 
+                  [for(i=[0:N-1]) 
+                     (reference[i]*poly[(i+k)%N]) ] ]*I,
+        optimal_poly = polygon_shift(fixpoly, max_index(val))
+      )
+    return_error? [optimal_poly, min(poly*(I*poly)-2*val)] :
     optimal_poly;
-
+    
 
 // Function: align_polygon()
 // Usage:
-//   newpoly = align_polygon(reference, poly, angles, [cp]);
+//   newpoly = align_polygon(reference, poly, angles, <cp>);
 // Description:
-//   Tries the list or range of angles to find a rotation of the specified polygon that best aligns
-//   with the reference polygon.  For each angle, the polygon is reindexed, which is a costly operation
+//   Tries the list or range of angles to find a rotation of the specified 2D polygon that best aligns
+//   with the reference 2D polygon.  For each angle, the polygon is reindexed, which is a costly operation
 //   so if run time is a problem, use a smaller sampling of angles.  Returns the rotated and reindexed
 //   polygon.
 // Arguments:
@@ -1546,9 +1657,11 @@ function reindex_polygon(reference, poly, return_error=false) =
 //   color("red") move_copies(scale(1.4,p=align_polygon(pentagon,hexagon,[0:10:359]))) circle(r=.1);
 //   move_copies(concat(pentagon,hexagon))circle(r=.1);
 function align_polygon(reference, poly, angles, cp) =
-    assert(is_path(reference) && is_path(poly))
-    assert(len(reference)==len(poly), "Polygons must be the same length to be aligned in align_polygon")
-    assert(is_num(angles[0]), "The `angle` parameter to align_polygon must be a range or vector")
+    assert(is_path(reference,dim=2) && is_path(poly,dim=2),
+           "Invalid polygon(s). " )
+    assert(len(reference)==len(poly), "The polygons must have the same length.")
+    assert( (is_vector(angles) && len(angles)>0) || valid_range(angles),
+            "The `angle` parameter must be a range or a non void list of numbers.")
     let(     // alignments is a vector of entries of the form: [polygon, error]
         alignments = [
             for(angle=angles) reindex_polygon(
@@ -1569,30 +1682,37 @@ function align_polygon(reference, poly, angles, cp) =
 //   Given a simple 3D planar polygon, returns the 3D coordinates of the polygon's centroid.
 //   If the polygon is self-intersecting, the results are undefined.
 function centroid(poly) =
-    len(poly[0])==2? (
-        sum([
+    assert( is_path(poly,dim=[2,3]), "The input must be a 2D or 3D polygon." )
+    len(poly[0])==2
+    ? sum([
             for(i=[0:len(poly)-1])
             let(segment=select(poly,i,i+1))
             det2(segment)*sum(segment)
-        ]) / 6 / polygon_area(poly)
-    ) : (
-        let(
-            n = plane_normal(plane_from_points(poly)),
-            p1 = vector_angle(n,UP)>15? vector_axis(n,UP) : vector_axis(n,RIGHT),
-            p2 = vector_axis(n,p1),
-            cp = mean(poly),
-            proj = project_plane(poly,cp,cp+p1,cp+p2),
-            cxy = centroid(proj)
-        ) lift_plane(cxy,cp,cp+p1,cp+p2)
-    );
-
+          ]) / 6 / polygon_area(poly)
+    : let( plane = plane_from_points(poly, fast=true) )
+      assert( !is_undef(plane), "The polygon must be planar." )
+      let(
+        n = plane_normal(plane),
+        val = sum([for(i=[1:len(poly)-2])
+                      let(
+                         v0 = poly[0],
+                         v1 = poly[i],
+                         v2 = poly[i+1],
+                         area = cross(v2-v0,v1-v0)*n
+                         )
+                      [ area, (v0+v1+v2)*area ]
+                  ] )
+          )
+      val[1]/val[0]/3;
+            
 
 // Function: point_in_polygon()
 // Usage:
-//   point_in_polygon(point, path, [eps])
+//   point_in_polygon(point, poly, <eps>)
 // Description:
 //   This function tests whether the given 2D point is inside, outside or on the boundary of
-//   the specified 2D polygon using the Winding Number method.
+//   the specified 2D polygon using either the Nonzero Winding rule or the Even-Odd rule.
+//   See https://en.wikipedia.org/wiki/Nonzero-rule and https://en.wikipedia.org/wiki/Even–odd_rule.
 //   The polygon is given as a list of 2D points, not including the repeated end point.
 //   Returns -1 if the point is outside the polyon.
 //   Returns 0 if the point is on the boundary.
@@ -1602,52 +1722,81 @@ function centroid(poly) =
 //   Rounding error may give mixed results for points on or near the boundary.
 // Arguments:
 //   point = The 2D point to check position of.
-//   path = The list of 2D path points forming the perimeter of the polygon.
+//   poly = The list of 2D path points forming the perimeter of the polygon.
+//   nonzero = The rule to use: true for "Nonzero" rule and false for "Even-Odd" (Default: true )
 //   eps = Acceptable variance.  Default: `EPSILON` (1e-9)
-function point_in_polygon(point, path, eps=EPSILON) =
-    // Original algorithm from http://geomalgorithms.com/a03-_inclusion.html
+function point_in_polygon(point, poly, eps=EPSILON, nonzero=true) =
+    // Original algorithms from http://geomalgorithms.com/a03-_inclusion.html
+    assert( is_vector(point,2) && is_path(poly,dim=2) && len(poly)>2,
+            "The point and polygon should be in 2D. The polygon should have more that 2 points." )
+    assert( is_finite(eps) && eps>=0, "Invalid tolerance." )
     // Does the point lie on any edges?  If so return 0.
-    sum([for(i=[0:1:len(path)-1]) let(seg=select(path,i,i+1)) if(!approx(seg[0],seg[1],eps=eps)) point_on_segment2d(point, seg, eps=eps)?1:0]) > 0? 0 :
-    // Otherwise compute winding number and return 1 for interior, -1 for exterior
-    sum([for(i=[0:1:len(path)-1]) let(seg=select(path,i,i+1)) if(!approx(seg[0],seg[1],eps=eps)) _point_above_below_segment(point, seg)]) != 0? 1 : -1;
+    let(
+        on_brd = [for(i=[0:1:len(poly)-1]) 
+                    let( seg = select(poly,i,i+1) ) 
+                    if( !approx(seg[0],seg[1],eps=EPSILON) ) 
+                        point_on_segment2d(point, seg, eps=eps)? 1:0 ]
+        )
+    sum(on_brd) > 0
+    ? 0 
+    :   nonzero
+        ?    // Compute winding number and return 1 for interior, -1 for exterior
+            let(
+                windchk = [for(i=[0:1:len(poly)-1]) 
+                            let(seg=select(poly,i,i+1)) 
+                            if(!approx(seg[0],seg[1],eps=eps)) 
+                                _point_above_below_segment(point, seg)
+                          ]
+                )
+            sum(windchk) != 0 ? 1 : -1
+        :   // or compute the crossings with the ray [point, point+[1,0]]
+            let( 
+              n  = len(poly),
+              cross = 
+                [for(i=[0:n-1])
+                    let( 
+                      p0 = poly[i]-point, 
+                      p1 = poly[(i+1)%n]-point
+                      )
+                    if( ( (p1.y>eps && p0.y<=0) || (p1.y<=0 && p0.y>eps) )
+                        &&  0 < p0.x - p0.y *(p1.x - p0.x)/(p1.y - p0.y) )
+                    1
+                ]
+            )
+            2*(len(cross)%2)-1;;
 
 
 // Function: polygon_is_clockwise()
 // Usage:
-//   polygon_is_clockwise(path);
+//   polygon_is_clockwise(poly);
 // Description:
 //   Return true if the given 2D simple polygon is in clockwise order, false otherwise.
 //   Results for complex (self-intersecting) polygon are indeterminate.
 // Arguments:
-//   path = The list of 2D path points for the perimeter of the polygon.
-function polygon_is_clockwise(path) =
-    assert(is_path(path) && len(path[0])==2, "Input must be a 2d path")
-    let(
-        minx = min(subindex(path,0)),
-        lowind = search(minx, path, 0, 0),
-        lowpts = select(path, lowind),
-        miny = min(subindex(lowpts, 1)),
-        extreme_sub = search(miny, lowpts, 1, 1)[0],
-        extreme = select(lowind,extreme_sub)
-    ) det2([select(path,extreme+1)-path[extreme], select(path, extreme-1)-path[extreme]])<0;
+//   poly = The list of 2D path points for the perimeter of the polygon.
+function polygon_is_clockwise(poly) =
+    assert(is_path(poly,dim=2), "Input should be a 2d path")
+    polygon_area(poly, signed=true)<0;
 
 
 // Function: clockwise_polygon()
 // Usage:
-//   clockwise_polygon(path);
+//   clockwise_polygon(poly);
 // Description:
 //   Given a 2D polygon path, returns the clockwise winding version of that path.
-function clockwise_polygon(path) =
-    polygon_is_clockwise(path)? path : reverse_polygon(path);
+function clockwise_polygon(poly) =
+    assert(is_path(poly,dim=2), "Input should be a 2d polygon")
+    polygon_area(poly, signed=true)<0 ? poly : reverse_polygon(poly);
 
 
 // Function: ccw_polygon()
 // Usage:
-//   ccw_polygon(path);
+//   ccw_polygon(poly);
 // Description:
-//   Given a 2D polygon path, returns the counter-clockwise winding version of that path.
-function ccw_polygon(path) =
-    polygon_is_clockwise(path)? reverse_polygon(path) : path;
+//   Given a 2D polygon poly, returns the counter-clockwise winding version of that poly.
+function ccw_polygon(poly) =
+    assert(is_path(poly,dim=2), "Input should be a 2d polygon")
+    polygon_area(poly, signed=true)<0 ? reverse_polygon(poly) : poly;
 
 
 // Function: reverse_polygon()
@@ -1656,6 +1805,7 @@ function ccw_polygon(path) =
 // Description:
 //   Reverses a polygon's winding direction, while still using the same start point.
 function reverse_polygon(poly) =
+    assert(is_path(poly), "Input should be a polygon")
     let(lp=len(poly)) [for (i=idx(poly)) poly[(lp-i)%lp]];
 
 
@@ -1666,8 +1816,9 @@ function reverse_polygon(poly) =
 //   Given a 3D planar polygon, returns a unit-length normal vector for the
 //   clockwise orientation of the polygon. 
 function polygon_normal(poly) =
+    assert(is_path(poly,dim=3), "Invalid 3D polygon." )
     let(
-        poly = path3d(cleanup_path(poly)),
+        poly = cleanup_path(poly),
         p0 = poly[0],
         n = sum([
             for (i=[1:1:len(poly)-2])
@@ -1772,6 +1923,9 @@ function _split_polygon_at_z(poly, z) =
 //   polys = A list of 3D polygons to split.
 //   xs = A list of scalar X values to split at.
 function split_polygons_at_each_x(polys, xs, _i=0) =
+    assert( is_consistent(polys) && is_path(poly[0],dim=3) ,
+            "The input list should contains only 3D polygons." )
+    assert( is_finite(xs), "The split value list should contain only numbers." )  
     _i>=len(xs)? polys :
     split_polygons_at_each_x(
         [
@@ -1779,7 +1933,7 @@ function split_polygons_at_each_x(polys, xs, _i=0) =
             each _split_polygon_at_x(poly, xs[_i])
         ], xs, _i=_i+1
     );
-
+    
 
 // Function: split_polygons_at_each_y()
 // Usage:
@@ -1790,6 +1944,9 @@ function split_polygons_at_each_x(polys, xs, _i=0) =
 //   polys = A list of 3D polygons to split.
 //   ys = A list of scalar Y values to split at.
 function split_polygons_at_each_y(polys, ys, _i=0) =
+//    assert( is_consistent(polys) && is_path(polys[0],dim=3) , // not all polygons should have the same length!!!
+  //          "The input list should contains only 3D polygons." )
+    assert( is_finite(ys) || is_vector(ys), "The split value list should contain only numbers." )  //***
     _i>=len(ys)? polys :
     split_polygons_at_each_y(
         [
@@ -1808,6 +1965,9 @@ function split_polygons_at_each_y(polys, ys, _i=0) =
 //   polys = A list of 3D polygons to split.
 //   zs = A list of scalar Z values to split at.
 function split_polygons_at_each_z(polys, zs, _i=0) =
+    assert( is_consistent(polys) && is_path(poly[0],dim=3) ,
+            "The input list should contains only 3D polygons." )
+    assert( is_finite(zs), "The split value list should contain only numbers." )  
     _i>=len(zs)? polys :
     split_polygons_at_each_z(
         [
@@ -1817,5 +1977,4 @@ function split_polygons_at_each_z(polys, zs, _i=0) =
     );
 
 
-
 // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap
diff --git a/hull.scad b/hull.scad
index da7d2c9..874e2c8 100644
--- a/hull.scad
+++ b/hull.scad
@@ -92,7 +92,7 @@ function hull2d_path(points) =
     assert(is_path(points,2),"Invalid input to hull2d_path")
     len(points) < 2 ? []
   : len(points) == 2 ? [0,1]
-  : let(tri=find_noncollinear_points(points, error=false))
+  : let(tri=noncollinear_triple(points, error=false))
     tri == [] ? _hull_collinear(points)
   : let(
         remaining = [ for (i = [0:1:len(points)-1]) if (i != tri[0] && i!=tri[1] && i!=tri[2]) i ],
@@ -170,7 +170,7 @@ function hull3d_faces(points) =
     assert(is_path(points,3),"Invalid input to hull3d_faces")
     len(points) < 3 ? list_range(len(points))
   : let ( // start with a single non-collinear triangle
-          tri = find_noncollinear_points(points, error=false)
+          tri = noncollinear_triple(points, error=false)
         )
     tri==[] ? _hull_collinear(points)
   : let(
@@ -250,7 +250,7 @@ function _find_conflicts(point, planes) = [
 
 
 function _find_first_noncoplanar(plane, points, i) = 
-    (i >= len(points) || !coplanar(plane, points[i]))? i :
+    (i >= len(points) || !points_on_plane([points[i]],plane))? i :
     _find_first_noncoplanar(plane, points, i+1);
 
 
diff --git a/joiners.scad b/joiners.scad
index 1ab34e6..1178d4e 100644
--- a/joiners.scad
+++ b/joiners.scad
@@ -564,6 +564,7 @@ module dovetail(gender, length, l, width, w, height, h, angle, slope, taper, bac
 }
 
 
+// Section: Tension Clips
 
 // h is total height above 0 of the nub
 // nub extends below xy plane by distance nub/2
@@ -780,5 +781,244 @@ module snap_pin_socket(size, r, radius, l,length, d,diameter,nub_depth, snap, fi
 
 
 
+// Module: rabbit_clip()
+// Usage:
+//   rabbit_clip(type, length, width, snap, thickness, depth, [compression], [clearance], [lock],
+//               [lock_clearance], [splineteps], [anchor], [orient], [spin])
+// Description:
+//   Creates a clip with two flexible ears to lock into a mating socket, or create a mask to produce the appropriate
+//   mating socket.  The clip can be made to insert and release easily, or to hold much better, or it can be
+//   created with locking flanges that will make it very hard or impossible to remove.  Unlike the snap pin, this clip
+//   is rectangular and can be made at any height, so a suitable clip could be very thin.  It's also possible to get a
+//   solid connection with a short pin.
+//   .
+//   The type parameters specifies whether to make a clip, a socket mask, or a double clip.  The length is the
+//   total nominal length of the clip.  (The actual length will be very close, but not equal to this.)  The width
+//   gives the nominal width of the clip, which is the actual width of the clip at its base.  The snap parameter
+//   gives the depth of the clip sides, which controls how easy the clip is to insert and remove.  The clip "ears" are
+//   made over-wide by the compression value.  A nonzero compression helps make the clip secure in its socket.
+//   The socket's width and length are increased by the clearance value which creates some space and can compensate
+//   for printing inaccuracy.  The socket will be slightly longer than the nominal width.  The thickness is the thickness
+//   curved line that forms the clip.  The clip depth is the amount the basic clip shape is extruded.  Be sure that you
+//   make the socket with a larger depth than the clip (try 0.4 mm) to allow ease of insertion of the clip.  The clearance
+//   value does not apply to the depth.  The splinesteps parameter increases the sampling of the clip curves.
+//   .
+//   By default clips appear with orient=UP and sockets with orient=DOWN.  
+//   .
+//   The first figure shows the dimensions of the rabbit clip.  The second figure shows the clip in red overlayed on
+//   its socket in yellow.  The left clip has a nonzero clearance, so its socket is bigger than the clip all around.
+//   The right hand locking clip has no clearance, but it has a lock clearance, which provides some space behind
+//   the lock to allow the clip to fit.  (Note that depending on your printer, this can be set to zero.)  
+//
+// Figure(2DMed):
+//   snap=1.5;
+//   comp=0.75;
+//   mid = 8.053;  // computed in rabbit_clip
+//   tip = [-4.58,18.03];
+//   translate([9,3]){
+//   back_half()
+//      rabbit_clip("pin",width=12, length=18, depth=1, thickness = 1, compression=comp, snap=snap, orient=BACK);
+//   color("blue"){
+//      stroke([[6,0],[6,18]],width=0.1);
+//      stroke([[6+comp, 12], [6+comp, 18]], width=.1);
+//   }
+//   color("red"){
+//      stroke([[6-snap,mid], [6,mid]], endcaps="arrow2",width=0.15);
+//      translate([6+.4,mid-.15])text("snap",size=1,valign="center");
+//      translate([6+comp/2,19.5])text("compression", size=1, halign="center");
+//      stroke([[6+comp/2,19.3], [6+comp/2,17.7]], endcap2="arrow2", width=.15);
+//      fwd(1.1)text("width",size=1,halign="center");
+//      xflip_copy()stroke([[2,-.7], [6,-.7]], endcap2="arrow2", width=.15);     
+//      move([-6.7,mid])rot(90)text("length", size=1, halign="center");
+//      stroke([[-7,10.3], [-7,18]], width=.15, endcap2="arrow2");
+//      stroke([[-7,0], [-7,5.8]], width=.15,endcap1="arrow2");
+//      stroke([tip, tip-[0,1]], width=.15);
+//      move([tip.x+2,19.5])text("thickness", halign="center",size=1);
+//      stroke([[tip.x+2, 19.3], tip+[.1,.1]], width=.15, endcap2="arrow2");
+//   }
+//   }
+//
+// Figure(2DMed):
+//   snap=1.5;
+//   comp=0;
+//   translate([29,3]){
+//   back_half()
+//      rabbit_clip("socket", width=12, length=18, depth=1, thickness = 1, compression=comp, snap=snap, orient=BACK,lock=true);
+//   color("red")back_half()
+//      rabbit_clip("pin",width=12, length=18, depth=1, thickness = 1, compression=comp, snap=snap,
+//               orient=BACK,lock=true,lock_clearance=1);
+//   }
+//   translate([9,3]){
+//   back_half()
+//      rabbit_clip("socket", clearance=.5,width=12, length=18, depth=1, thickness = 1,
+//                  compression=comp, snap=snap, orient=BACK,lock=false);
+//   color("red")back_half()
+//      rabbit_clip("pin",width=12, length=18, depth=1, thickness = 1, compression=comp, snap=snap,
+//               orient=BACK,lock=false,lock_clearance=1);
+//   }
+// Arguments:
+//   type = One of "pin",  "socket", "male", "female" or "double" to specify what to make.
+//   length = nominal clip length
+//   width = nominal clip width
+//   snap = depth of hollow on the side of the clip
+//   thickness = thickness of the clip "line"
+//   depth = amount to extrude clip (give extra room for the socket, about 0.4mm)
+//   compression = excess width at the "ears" to lock more tightly.  Default: 0.1
+//   clearance = extra space in the socket for easier insertion.  Default: 0.1
+//   lock = set to true to make a locking clip that may be irreversible.  Default: false
+//   lock_clearance = give clearance for the lock.  Default: 0
+//   splinesteps = number of samples in the curves of the clip.  Default: 8
+//   anchor = anchor point for clip
+//   orient = clip orientation.  Default: UP for pins, DOWN for sockets
+//   spin = spin the clip.  Default: 0
+//
+// Example:  Here are several sizes that work printed in PLA on a Prusa MK3, with default clearance of 0.1 and a depth of 5
+//   module test_pair(length, width, snap, thickness, compression, lock=false)
+//   {
+//     depth = 5;
+//     extra_depth = 10;// Change this to 0.4 for closed sockets
+//     cuboid([max(width+5,12),12, depth], chamfer=.5, edges=[FRONT,"Y"], anchor=BOTTOM)
+//         attach(BACK)
+//           rabbit_clip(type="pin",length=length, width=width,snap=snap,thickness=thickness,depth=depth,
+//                       compression=compression,lock=lock);
+//     right(width+13)
+//     diff("remove")
+//         cuboid([width+8,max(12,length+2),depth+3], chamfer=.5, edges=[FRONT,"Y"], anchor=BOTTOM)
+//             attach(BACK)
+//               rabbit_clip(type="socket",length=length, width=width,snap=snap,thickness=thickness,depth=depth+extra_depth,
+//                           lock=lock,compression=0,$tags="remove");
+//   }
+//   left(37)ydistribute(spacing=28){
+//     test_pair(length=6, width=7, snap=0.25, thickness=0.8, compression=0.1);
+//     test_pair(length=3.5, width=7, snap=0.1, thickness=0.8, compression=0.1);  // snap = 0.2 gives a firmer connection
+//     test_pair(length=3.5, width=5, snap=0.1, thickness=0.8, compression=0.1);  // hard to take apart
+//   }
+//   right(17)ydistribute(spacing=28){
+//     test_pair(length=12, width=10, snap=1, thickness=1.2, compression=0.2);
+//     test_pair(length=8, width=7, snap=0.75, thickness=0.8, compression=0.2, lock=true); // With lock, very firm and irreversible
+//     test_pair(length=8, width=7, snap=0.75, thickness=0.8, compression=0.2, lock=true); // With lock, very firm and irreversible
+//   }
+// Example: Double clip to connect two sockets
+//   rabbit_clip("double",length=8, width=7, snap=0.75, thickness=0.8, compression=0.2,depth=5);
+// Example:  A modified version of the clip that acts like a backpack strap clip, where it locks tightly but you can squeeze to release. 
+//   cuboid([25,15,5],anchor=BOTTOM)
+//       attach(BACK)rabbit_clip("pin", length=25, width=25, thickness=1.5, snap=2, compression=0, lock=true, depth=5, lock_clearance=3);
+//   left(32)
+//   diff("remove")
+//   cuboid([30,30,11],orient=BACK,anchor=BACK){
+//       attach(BACK)rabbit_clip("socket", length=25, width=25, thickness=1.5, snap=2, compression=0, lock=true, depth=5.5, lock_clearance=3,$tags="remove");
+//       xflip_copy()
+//         position(FRONT+LEFT)
+//         xscale(0.8)
+//         zcyl(l=20,r=13.5, $tags="remove",$fn=64);
+//   }
+module rabbit_clip(type, length, width,  snap, thickness, depth, compression=0.1,  clearance=.1, lock=false, lock_clearance=0,
+                   splinesteps=8, anchor, orient, spin=0)
+{
+  assert(is_num(width) && width>0,"Width must be a positive value");
+  assert(is_num(length) && length>0, "Length must be a positive value");
+  assert(is_num(thickness) && thickness>0, "Thickness must be a positive value");  
+  assert(is_num(snap) && snap>=0, "Snap must be a non-negative value");
+  assert(is_num(depth) && depth>0, "Depth must be a positive value");
+  assert(is_num(compression) && compression >= 0, "Compression must be a nonnegative value");
+  assert(is_bool(lock));
+  assert(is_num(lock_clearance));
+  legal_types = ["pin","socket","male","female","double"];
+  assert(in_list(type,legal_types),str("type must be one of ",legal_types));
+  
+  if (type=="double") {
+    attachable(size=[width+2*compression, depth, 2*length], anchor=default(anchor,BACK), spin=spin, orient=default(orient,BACK)){
+      union(){
+        rabbit_clip("pin", length=length, width=width, snap=snap, thickness=thickness, depth=depth, compression=compression,
+                    lock=lock, anchor=BOTTOM, orient=UP);
+        rabbit_clip("pin", length=length, width=width, snap=snap, thickness=thickness, depth=depth, compression=compression,
+                    lock=lock, anchor=BOTTOM, orient=DOWN);
+        cuboid([width-thickness, depth, thickness]);
+      }
+      children();
+    }
+  } else {
+    anchor = default(anchor,BOTTOM);
+    is_pin = in_list(type,["pin","male"]);
+    default_overlap = 0.01 * (is_pin?1:-1);    // Shift by this much to undo default overlap
+    extra = 0.02;  // Amount of extension below nominal based position for the socket, must exceed default overlap of 0.01
+    clearance = is_pin ? 0 : clearance;
+    compression = is_pin ? compression : 0;
+    orient =  is_def(orient) ? orient
+            : is_pin ? UP
+            : DOWN;
+    earwidth = 2*thickness+snap;
+    point_length = earwidth/2.15;
+    // The adjustment is using cos(theta)*earwidth/2 and sin(theta)*point_length, but the computation
+    // is obscured because theta is atan(length/2/snap)
+    scaled_len = length - 0.5 * (earwidth * snap + point_length * length) / sqrt(sqr(snap)+sqr(length/2));
+    bottom_pt = [0,max(scaled_len*0.15+thickness, 2*thickness)];
+    ctr = [width/2,scaled_len] + line_normal([width/2-snap, scaled_len/2], [width/2, scaled_len]) * earwidth/2;
+    inside_pt = circle_circle_tangents(bottom_pt, 0, ctr, earwidth/2)[0][1];
+    sidepath =[
+               [width/2,0],
+               [width/2-snap,scaled_len/2],
+               [width/2+(is_pin?compression:0), scaled_len],
+               ctr - point_length * line_normal([width/2,scaled_len], inside_pt),
+               inside_pt
+              ];
+    fullpath = concat(
+                      sidepath,
+                      [bottom_pt], 
+                      reverse(apply(xflip(),sidepath))
+                      );
+    assert(fullpath[4].y < fullpath[3].y, "Pin is too wide for its length");
+    
+    snapmargin = -snap + select(sidepath,-1).x;// - compression;
+    if (is_pin){
+      if (snapmargin<0) echo("WARNING: The snap is too large for the clip to squeeze to fit its socket")
+      echo(snapmargin=snapmargin);
+    }
+    // Force tangent to be vertical at the outer edge of the clip to avoid overshoot
+    fulltangent = list_set(path_tangents(fullpath, uniform=false),[2,8], [[0,1],[0,-1]]);
+
+    subset = is_pin ? [0:10] : [0,1,2,3, 7,8,9,10];  // Remove internal points from the socket
+    tangent = select(fulltangent, subset);
+    path = select(fullpath, subset);
+    
+    socket_smooth = .04;
+    pin_smooth = [.075, .075, .15, .12, .06];
+    smoothing = is_pin
+                  ? concat(pin_smooth, reverse(pin_smooth))
+                  : let(side_smooth=select(pin_smooth, 0, 2))
+                    concat(side_smooth, [socket_smooth], reverse(side_smooth));
+    bez = path_to_bezier(path,relsize=smoothing,tangents=tangent);
+    rounded = bezier_polyline(bez,splinesteps=splinesteps);
+    bounds = pointlist_bounds(rounded);
+    kk = search([bounds[1].y], subindex(rounded,1));
+    echo(rounded[kk[0]]);
+    extrapt = is_pin ? [] : [rounded[0] - [0,extra]];
+    finalpath = is_pin ? rounded
+                       : let(withclearance=offset(rounded, r=-clearance))
+                         concat( [[withclearance[0].x,-extra]],
+                                 withclearance,
+                                 [[-withclearance[0].x,-extra]]);
+    attachable(size=[bounds[1].x-bounds[0].x, depth, bounds[1].y-bounds[0].y], anchor=anchor, spin=spin, orient=orient){
+      xrot(90)
+        translate([0,-(bounds[1].y-bounds[0].y)/2+default_overlap,-depth/2])
+        linear_extrude(height=depth, convexity=10) {
+            if (lock)
+              xflip_copy()
+              right(clearance)
+              polygon([sidepath[1]+[-thickness/10,lock_clearance],
+                       sidepath[2],
+                       [sidepath[2].x,sidepath[1].y+lock_clearance]]);
+            if (is_pin)
+              offset_stroke(finalpath, width=[thickness,0]);
+            else
+              polygon(finalpath);
+        }
+      children();
+    }
+  }
+}
+
+
+
 
 // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap
diff --git a/masks.scad b/masks.scad
index 4a75886..894e50d 100644
--- a/masks.scad
+++ b/masks.scad
@@ -12,8 +12,8 @@
 
 // Module: angle_pie_mask()
 // Usage:
-//   angle_pie_mask(r|d, l, ang);
-//   angle_pie_mask(r1|d1, r2|d2, l, ang);
+//   angle_pie_mask(r|d, l, ang, [excess]);
+//   angle_pie_mask(r1|d1, r2|d2, l, ang, [excess]);
 // Description:
 //   Creates a pie wedge shape that can be used to mask other shapes.
 // Arguments:
@@ -25,6 +25,7 @@
 //   d = Diameter of circle wedge is created from. (optional)
 //   d1 = Bottom diameter of cone that wedge is created from.  (optional)
 //   d2 = Upper diameter of cone that wedge is created from. (optional)
+//   excess = The extra thickness of the mask.  Default: `0.1`.
 //   anchor = Translate so anchor point is at origin (0,0,0).  See [anchor](attachments.scad#anchor).  Default: `CENTER`
 //   spin = Rotate this many degrees around the Z axis after anchor.  See [spin](attachments.scad#spin).  Default: `0`
 //   orient = Vector to rotate top towards, after spin.  See [orient](attachments.scad#orient).  Default: `UP`
@@ -34,14 +35,14 @@ module angle_pie_mask(
     ang=45, l=undef,
     r=undef, r1=undef, r2=undef,
     d=undef, d1=undef, d2=undef,
-    h=undef,
+    h=undef, excess=0.1,
     anchor=CENTER, spin=0, orient=UP
 ) {
     l = first_defined([l, h, 1]);
     r1 = get_radius(r1=r1, r=r, d1=d1, d=d, dflt=10);
     r2 = get_radius(r1=r2, r=r, d1=d2, d=d, dflt=10);
     attachable(anchor,spin,orient, r1=r1, r2=r2, l=l) {
-        pie_slice(ang=ang, l=l+0.1, r1=r1, r2=r2, anchor=CENTER);
+        pie_slice(ang=ang, l=l+excess, r1=r1, r2=r2, anchor=CENTER);
         children();
     }
 }
@@ -49,13 +50,13 @@ module angle_pie_mask(
 
 // Module: cylinder_mask()
 // Usage: Mask objects
-//   cylinder_mask(l, r|d, chamfer, [chamfang], [from_end], [circum], [overage], [ends_only]);
-//   cylinder_mask(l, r|d, rounding, [circum], [overage], [ends_only]);
-//   cylinder_mask(l, r|d, [chamfer1|rounding1], [chamfer2|rounding2], [chamfang1], [chamfang2], [from_end], [circum], [overage], [ends_only]);
+//   cylinder_mask(l, r|d, chamfer, [chamfang], [from_end], [circum], [excess], [ends_only]);
+//   cylinder_mask(l, r|d, rounding, [circum], [excess], [ends_only]);
+//   cylinder_mask(l, r|d, [chamfer1|rounding1], [chamfer2|rounding2], [chamfang1], [chamfang2], [from_end], [circum], [excess], [ends_only]);
 // Usage: Masking operators
-//   cylinder_mask(l, r|d, chamfer, [chamfang], [from_end], [circum], [overage], [ends_only]) ...
-//   cylinder_mask(l, r|d, rounding, [circum], [overage], [ends_only]) ...
-//   cylinder_mask(l, r|d, [chamfer1|rounding1], [chamfer2|rounding2], [chamfang1], [chamfang2], [from_end], [circum], [overage], [ends_only]) ...
+//   cylinder_mask(l, r|d, chamfer, [chamfang], [from_end], [circum], [excess], [ends_only]) ...
+//   cylinder_mask(l, r|d, rounding, [circum], [excess], [ends_only]) ...
+//   cylinder_mask(l, r|d, [chamfer1|rounding1], [chamfer2|rounding2], [chamfang1], [chamfang2], [from_end], [circum], [excess], [ends_only]) ...
 // Description:
 //   If passed children, bevels/chamfers and/or rounds one or both
 //   ends of the origin-centered cylindrical region specified.  If
@@ -83,7 +84,7 @@ module angle_pie_mask(
 //   rounding2 = The radius of the rounding on the axis-positive end of the region.
 //   circum = If true, region will circumscribe the circle of the given radius/diameter.
 //   from_end = If true, chamfer/bevel size is measured from end of region.  If false, chamfer/bevel is measured outset from the radius of the region.  (Default: false)
-//   overage = The extra thickness of the mask.  Default: `10`.
+//   excess = The extra thickness of the mask.  Default: `10`.
 //   ends_only = If true, only mask the ends and not around the middle of the cylinder.
 //   anchor = Translate so anchor point is at origin (0,0,0).  See [anchor](attachments.scad#anchor).  Default: `CENTER`
 //   spin = Rotate this many degrees around the Z axis after anchor.  See [spin](attachments.scad#spin).  Default: `0`
@@ -105,7 +106,7 @@ module cylinder_mask(
     chamfang=undef, chamfang1=undef, chamfang2=undef,
     rounding=undef, rounding1=undef, rounding2=undef,
     circum=false, from_end=false,
-    overage=10, ends_only=false,
+    excess=10, ends_only=false,
     anchor=CENTER, spin=0, orient=UP
 ) {
     r1 = get_radius(r=r, d=d, r1=r1, d1=d1, dflt=1);
@@ -132,12 +133,12 @@ module cylinder_mask(
                     chlen1 = cham1 / (from_end? 1 : tan(ang1));
                     chlen2 = cham2 / (from_end? 1 : tan(ang2));
                     if (!ends_only) {
-                        cylinder(r=maxd+overage, h=l+2*overage, center=true);
+                        cylinder(r=maxd+excess, h=l+2*excess, center=true);
                     } else {
-                        if (cham2>0) up(l/2-chlen2) cylinder(r=maxd+overage, h=chlen2+overage, center=false);
-                        if (cham1>0) down(l/2+overage) cylinder(r=maxd+overage, h=chlen1+overage, center=false);
-                        if (fil2>0) up(l/2-fil2) cylinder(r=maxd+overage, h=fil2+overage, center=false);
-                        if (fil1>0) down(l/2+overage) cylinder(r=maxd+overage, h=fil1+overage, center=false);
+                        if (cham2>0) up(l/2-chlen2) cylinder(r=maxd+excess, h=chlen2+excess, center=false);
+                        if (cham1>0) down(l/2+excess) cylinder(r=maxd+excess, h=chlen1+excess, center=false);
+                        if (fil2>0) up(l/2-fil2) cylinder(r=maxd+excess, h=fil2+excess, center=false);
+                        if (fil1>0) down(l/2+excess) cylinder(r=maxd+excess, h=fil1+excess, center=false);
                     }
                 }
                 cyl(r1=sc*r1, r2=sc*r2, l=l, chamfer1=cham1, chamfer2=cham2, chamfang1=ang1, chamfang2=ang2, from_end=from_end, rounding1=fil1, rounding2=fil2);
@@ -154,14 +155,15 @@ module cylinder_mask(
 
 // Module: chamfer_mask()
 // Usage:
-//   chamfer_mask(l, chamfer);
+//   chamfer_mask(l, chamfer, [excess]);
 // Description:
 //   Creates a shape that can be used to chamfer a 90 degree edge.
 //   Difference it from the object to be chamfered.  The center of
 //   the mask object should align exactly with the edge to be chamfered.
 // Arguments:
 //   l = Length of mask.
-//   chamfer = Size of chamfer
+//   chamfer = Size of chamfer.
+//   excess = The extra amount to add to the length of the mask so that it differences away from other shapes cleanly.  Default: `0.1`
 //   anchor = Translate so anchor point is at origin (0,0,0).  See [anchor](attachments.scad#anchor).  Default: `CENTER`
 //   spin = Rotate this many degrees around the Z axis after anchor.  See [spin](attachments.scad#spin).  Default: `0`
 //   orient = Vector to rotate top towards, after spin.  See [orient](attachments.scad#orient).  Default: `UP`
@@ -170,9 +172,9 @@ module cylinder_mask(
 //       cube(50, anchor=BOTTOM+FRONT);
 //       #chamfer_mask(l=50, chamfer=10, orient=RIGHT);
 //   }
-module chamfer_mask(l=1, chamfer=1, anchor=CENTER, spin=0, orient=UP) {
+module chamfer_mask(l=1, chamfer=1, excess=0.1, anchor=CENTER, spin=0, orient=UP) {
     attachable(anchor,spin,orient, size=[chamfer*2, chamfer*2, l]) {
-        cylinder(r=chamfer, h=l+0.1, center=true, $fn=4);
+        cylinder(r=chamfer, h=l+excess, center=true, $fn=4);
         children();
     }
 }
@@ -180,14 +182,15 @@ module chamfer_mask(l=1, chamfer=1, anchor=CENTER, spin=0, orient=UP) {
 
 // Module: chamfer_mask_x()
 // Usage:
-//   chamfer_mask_x(l, chamfer, [anchor]);
+//   chamfer_mask_x(l, chamfer, [excess]);
 // Description:
 //   Creates a shape that can be used to chamfer a 90 degree edge along the X axis.
 //   Difference it from the object to be chamfered.  The center of the mask
 //   object should align exactly with the edge to be chamfered.
 // Arguments:
-//   l = Height of mask
-//   chamfer = size of chamfer
+//   l = Length of mask.
+//   chamfer = Size of chamfer.
+//   excess = The extra amount to add to the length of the mask so that it differences away from other shapes cleanly.  Default: `0.1`
 //   anchor = Translate so anchor point is at origin (0,0,0).  See [anchor](attachments.scad#anchor).  Default: `CENTER`
 //   spin = Rotate this many degrees around the X axis after anchor.  See [spin](attachments.scad#spin).  Default: `0`
 // Example:
@@ -195,21 +198,22 @@ module chamfer_mask(l=1, chamfer=1, anchor=CENTER, spin=0, orient=UP) {
 //       cube(50, anchor=BOTTOM+FRONT);
 //       #chamfer_mask_x(l=50, chamfer=10);
 //   }
-module chamfer_mask_x(l=1.0, chamfer=1.0, anchor=CENTER, spin=0) {
-    chamfer_mask(l=l, chamfer=chamfer, anchor=anchor, spin=spin, orient=RIGHT) children();
+module chamfer_mask_x(l=1.0, chamfer=1.0, excess=0.1, anchor=CENTER, spin=0) {
+    chamfer_mask(l=l, chamfer=chamfer, excess=excess, anchor=anchor, spin=spin, orient=RIGHT) children();
 }
 
 
 // Module: chamfer_mask_y()
 // Usage:
-//   chamfer_mask_y(l, chamfer, [anchor]);
+//   chamfer_mask_y(l, chamfer, [excess]);
 // Description:
 //   Creates a shape that can be used to chamfer a 90 degree edge along the Y axis.
 //   Difference it from the object to be chamfered.  The center of the mask
 //   object should align exactly with the edge to be chamfered.
 // Arguments:
-//   l = Height of mask
-//   chamfer = size of chamfer
+//   l = Length of mask.
+//   chamfer = Size of chamfer.
+//   excess = The extra amount to add to the length of the mask so that it differences away from other shapes cleanly.  Default: `0.1`
 //   anchor = Translate so anchor point is at origin (0,0,0).  See [anchor](attachments.scad#anchor).  Default: `CENTER`
 //   spin = Rotate this many degrees around the Y axis after anchor.  See [spin](attachments.scad#spin).  Default: `0`
 // Example:
@@ -217,21 +221,22 @@ module chamfer_mask_x(l=1.0, chamfer=1.0, anchor=CENTER, spin=0) {
 //       cube(50, anchor=BOTTOM+RIGHT);
 //       #chamfer_mask_y(l=50, chamfer=10);
 //   }
-module chamfer_mask_y(l=1.0, chamfer=1.0, anchor=CENTER, spin=0) {
-    chamfer_mask(l=l, chamfer=chamfer, anchor=anchor, spin=spin, orient=BACK) children();
+module chamfer_mask_y(l=1.0, chamfer=1.0, excess=0.1, anchor=CENTER, spin=0) {
+    chamfer_mask(l=l, chamfer=chamfer, excess=excess, anchor=anchor, spin=spin, orient=BACK) children();
 }
 
 
 // Module: chamfer_mask_z()
 // Usage:
-//   chamfer_mask_z(l, chamfer, [anchor]);
+//   chamfer_mask_z(l, chamfer, [excess]);
 // Description:
 //   Creates a shape that can be used to chamfer a 90 degree edge along the Z axis.
 //   Difference it from the object to be chamfered.  The center of the mask
 //   object should align exactly with the edge to be chamfered.
 // Arguments:
-//   l = Height of mask
-//   chamfer = size of chamfer
+//   l = Length of mask.
+//   chamfer = Size of chamfer.
+//   excess = The extra amount to add to the length of the mask so that it differences away from other shapes cleanly.  Default: `0.1`
 //   anchor = Translate so anchor point is at origin (0,0,0).  See [anchor](attachments.scad#anchor).  Default: `CENTER`
 //   spin = Rotate this many degrees around the Z axis after anchor.  See [spin](attachments.scad#spin).  Default: `0`
 // Example:
@@ -239,8 +244,8 @@ module chamfer_mask_y(l=1.0, chamfer=1.0, anchor=CENTER, spin=0) {
 //       cube(50, anchor=FRONT+RIGHT);
 //       #chamfer_mask_z(l=50, chamfer=10);
 //   }
-module chamfer_mask_z(l=1.0, chamfer=1.0, anchor=CENTER, spin=0) {
-    chamfer_mask(l=l, chamfer=chamfer, anchor=anchor, spin=spin, orient=UP) children();
+module chamfer_mask_z(l=1.0, chamfer=1.0, excess=0.1, anchor=CENTER, spin=0) {
+    chamfer_mask(l=l, chamfer=chamfer, excess=excess, anchor=anchor, spin=spin, orient=UP) children();
 }
 
 
@@ -313,7 +318,7 @@ module chamfer_cylinder_mask(r=undef, d=undef, chamfer=0.25, ang=45, from_end=fa
 
 // Module: chamfer_hole_mask()
 // Usage:
-//   chamfer_hole_mask(r|d, chamfer, [ang], [from_end]);
+//   chamfer_hole_mask(r|d, chamfer, [ang], [from_end], [excess]);
 // Description:
 //   Create a mask that can be used to bevel/chamfer the end of a cylindrical hole.
 //   Difference it from the hole to be chamfered.  The center of the mask object
@@ -324,7 +329,7 @@ module chamfer_cylinder_mask(r=undef, d=undef, chamfer=0.25, ang=45, from_end=fa
 //   chamfer = Size of the chamfer. (Default: 0.25)
 //   ang = Angle of chamfer in degrees from vertical.  (Default: 45)
 //   from_end = If true, chamfer size is measured from end of hole.  If false, chamfer is measured outset from the radius of the hole.  (Default: false)
-//   overage = The extra thickness of the mask.  Default: `0.1`.
+//   excess = The extra thickness of the mask.  Default: `0.1`.
 //   anchor = Translate so anchor point is at origin (0,0,0).  See [anchor](attachments.scad#anchor).  Default: `CENTER`
 //   spin = Rotate this many degrees around the Z axis after anchor.  See [spin](attachments.scad#spin).  Default: `0`
 //   orient = Vector to rotate top towards, after spin.  See [orient](attachments.scad#orient).  Default: `UP`
@@ -341,8 +346,8 @@ module chamfer_cylinder_mask(r=undef, d=undef, chamfer=0.25, ang=45, from_end=fa
 //       up(50) chamfer_hole_mask(d=50, chamfer=10);
 //   }
 // Example:
-//   chamfer_hole_mask(d=100, chamfer=25, ang=30, overage=10);
-module chamfer_hole_mask(r=undef, d=undef, chamfer=0.25, ang=45, from_end=false, overage=0.1, anchor=CENTER, spin=0, orient=UP)
+//   chamfer_hole_mask(d=100, chamfer=25, ang=30, excess=10);
+module chamfer_hole_mask(r=undef, d=undef, chamfer=0.25, ang=45, from_end=false, excess=0.1, anchor=CENTER, spin=0, orient=UP)
 {
     r = get_radius(r=r, d=d, dflt=1);
     h = chamfer * (from_end? 1 : tan(90-ang));
@@ -350,7 +355,7 @@ module chamfer_hole_mask(r=undef, d=undef, chamfer=0.25, ang=45, from_end=false,
     $fn = segs(r);
     attachable(anchor,spin,orient, r1=r, r2=r2, l=h*2) {
         union() {
-            cylinder(r=r2, h=overage, center=false);
+            cylinder(r=r2, h=excess, center=false);
             down(h) cylinder(r1=r, r2=r2, h=h, center=false);
         }
         children();
@@ -735,14 +740,14 @@ module rounding_corner_mask(r=1.0, anchor=CENTER, spin=0, orient=UP)
 //   }
 module rounding_cylinder_mask(r=1.0, rounding=0.25)
 {
-    cylinder_mask(l=rounding*3, r=r, rounding2=rounding, overage=rounding, ends_only=true, anchor=TOP);
+    cylinder_mask(l=rounding*3, r=r, rounding2=rounding, excess=rounding, ends_only=true, anchor=TOP);
 }
 
 
 
 // Module: rounding_hole_mask()
 // Usage:
-//   rounding_hole_mask(r|d, rounding);
+//   rounding_hole_mask(r|d, rounding, [excess]);
 // Description:
 //   Create a mask that can be used to round the edge of a circular hole.
 //   Difference it from the hole to be rounded.  The center of the
@@ -752,7 +757,7 @@ module rounding_cylinder_mask(r=1.0, rounding=0.25)
 //   r = Radius of hole.
 //   d = Diameter of hole to rounding.
 //   rounding = Radius of the rounding. (Default: 0.25)
-//   overage = The extra thickness of the mask.  Default: `0.1`.
+//   excess = The extra thickness of the mask.  Default: `0.1`.
 //   anchor = Translate so anchor point is at origin (0,0,0).  See [anchor](attachments.scad#anchor).  Default: `CENTER`
 //   spin = Rotate this many degrees around the Z axis after anchor.  See [spin](attachments.scad#spin).  Default: `0`
 //   orient = Vector to rotate top towards, after spin.  See [orient](attachments.scad#orient).  Default: `UP`
@@ -770,13 +775,13 @@ module rounding_cylinder_mask(r=1.0, rounding=0.25)
 //   }
 // Example:
 //   rounding_hole_mask(r=40, rounding=20, $fa=2, $fs=2);
-module rounding_hole_mask(r=undef, d=undef, rounding=0.25, overage=0.1, anchor=CENTER, spin=0, orient=UP)
+module rounding_hole_mask(r=undef, d=undef, rounding=0.25, excess=0.1, anchor=CENTER, spin=0, orient=UP)
 {
     r = get_radius(r=r, d=d, dflt=1);
     attachable(anchor,spin,orient, r=r+rounding, l=2*rounding) {
         rotate_extrude(convexity=4) {
             difference() {
-                right(r-overage) fwd(rounding) square(rounding+overage, center=false);
+                right(r-excess) fwd(rounding) square(rounding+excess, center=false);
                 right(r+rounding) fwd(rounding) circle(r=rounding);
             }
         }
diff --git a/math.scad b/math.scad
index 55e662b..90182f1 100644
--- a/math.scad
+++ b/math.scad
@@ -36,7 +36,7 @@ NAN = acos(2);  // The value `nan`, useful for comparisons.
 function sqr(x) = 
     is_list(x) ? [for(val=x) sqr(val)] : 
     is_finite(x) ? x*x :
-    assert(is_finite(x) || is_vector(x), "Input is not neither a number nor a list of numbers.");
+    assert(is_finite(x) || is_vector(x), "Input is not a number nor a list of numbers.");
 
 
 // Function: log2()
@@ -84,7 +84,7 @@ function hypot(x,y,z=0) =
 //   y = factorial(6);  // Returns: 720
 //   z = factorial(9);  // Returns: 362880
 function factorial(n,d=0) =
-    assert(is_int(n) && is_int(d) && n>=0 && d>=0, "Factorial is not defined for negative numbers")
+    assert(is_int(n) && is_int(d) && n>=0 && d>=0, "Factorial is defined only for non negative integers")
     assert(d<=n, "d cannot be larger than n")
     product([1,for (i=[n:-1:d+1]) i]);
 
@@ -164,7 +164,7 @@ function binomial_coefficient(n,k) =
 function lerp(a,b,u) =
     assert(same_shape(a,b), "Bad or inconsistent inputs to lerp")
     is_finite(u)? (1-u)*a + u*b :
-    assert(is_finite(u) || is_vector(u) || valid_range(u), "Input u to lerp must be a number, vector, or range.")
+    assert(is_finite(u) || is_vector(u) || valid_range(u), "Input u to lerp must be a number, vector, or valid range.")
     [for (v = u) (1-v)*a + v*b ];
 
 
@@ -387,12 +387,13 @@ function modang(x) =
 //   modrange(90,270,360, step=-45);  // Returns: [90,45,0,315,270]
 //   modrange(270,90,360, step=-45);  // Returns: [270,225,180,135,90]
 function modrange(x, y, m, step=1) =
-    assert( is_finite(x+y+step+m) && !approx(m,0), "Input must be finite numbers. The module value cannot be zero.")
+    assert( is_finite(x+y+step+m) && !approx(m,0), "Input must be finite numbers and the module value cannot be zero." )
     let(
         a = posmod(x, m),
         b = posmod(y, m),
-        c = step>0? (a>b? b+m : b) : (a<b? b-m : b)
-    ) [for (i=[a:step:c]) (i%m+m)%m];
+        c = step>0? (a>b? b+m : b) 
+            : (a<b? b-m : b)
+    ) [for (i=[a:step:c]) (i%m+m)%m ];
 
 
 
@@ -536,9 +537,13 @@ function _sum(v,_total,_i=0) = _i>=len(v) ? _total : _sum(v,_total+v[_i], _i+1);
 //   cumsum([2,2,2]);  // returns [2,4,6]
 //   cumsum([1,2,3]);  // returns [1,3,6]
 //   cumsum([[1,2,3], [3,4,5], [5,6,7]]);  // returns [[1,2,3], [4,6,8], [9,12,15]]
-function cumsum(v,_i=0,_acc=[]) =
+function cumsum(v) =
+    assert(is_consistent(v), "The input is not consistent." )
+    _cumsum(v,_i=0,_acc=[]);
+
+function _cumsum(v,_i=0,_acc=[]) =
     _i==len(v) ? _acc :
-    cumsum(
+    _cumsum(
         v, _i+1,
         concat(
             _acc,
@@ -598,7 +603,7 @@ function deltas(v) =
 // Description:
 //   Returns the product of all entries in the given list.
 //   If passed a list of vectors of same dimension, returns a vector of products of each part.
-//   If passed a list of square matrices, returns a the resulting product matrix.
+//   If passed a list of square matrices, returns the resulting product matrix.
 // Arguments:
 //   v = The list to get the product of.
 // Example:
@@ -606,7 +611,7 @@ function deltas(v) =
 //   product([[1,2,3], [3,4,5], [5,6,7]]);  // returns [15, 48, 105]
 function product(v) = 
     assert( is_vector(v) || is_matrix(v) || ( is_matrix(v[0],square=true) && is_consistent(v)), 
-            "Invalid input.")
+    "Invalid input.")
     _product(v, 1, v[0]);
 
 function _product(v, i=0, _tot) = 
@@ -641,17 +646,6 @@ function mean(v) =
     sum(v)/len(v);
 
 
-// Function: median()
-// Usage:
-//   x = median(v);
-// Description:
-//   Given a list of numbers or vectors, finds the median value or midpoint.
-//   If passed a list of vectors, returns the vector of the median of each component.
-function median(v) =
-    is_vector(v) ? (min(v)+max(v))/2 :
-    is_matrix(v) ? [for(ti=transpose(v))  (min(ti)+max(ti))/2 ]
-    :   assert(false , "Invalid input.");
-
 // Function: convolve()
 // Usage:
 //   x = convolve(p,q);
@@ -681,7 +675,7 @@ function convolve(p,q) =
 // Usage: linear_solve(A,b)
 // Description:
 //   Solves the linear system Ax=b.  If A is square and non-singular the unique solution is returned.  If A is overdetermined
-//   the least squares solution is returned.  If A is underdetermined, the minimal norm solution is returned.
+//   the least squares solution is returned. If A is underdetermined, the minimal norm solution is returned.
 //   If A is rank deficient or singular then linear_solve returns [].  If b is a matrix that is compatible with A
 //   then the problem is solved for the matrix valued right hand side and a matrix is returned.  Note that if you 
 //   want to solve Ax=b1 and Ax=b2 that you need to form the matrix transpose([b1,b2]) for the right hand side and then
@@ -692,7 +686,7 @@ function linear_solve(A,b) =
         m = len(A),
         n = len(A[0])
     )
-    assert(is_vector(b,m) || is_matrix(b,m),"Incompatible matrix and right hand side")
+    assert(is_vector(b,m) || is_matrix(b,m),"Invalid right hand side or incompatible with the matrix")
     let (
         qr = m<n? qr_factor(transpose(A)) : qr_factor(A),
         maxdim = max(n,m),
@@ -702,9 +696,11 @@ function linear_solve(A,b) =
         zeros = [for(i=[0:mindim-1]) if (approx(R[i][i],0)) i]
     )
     zeros != [] ? [] :
-    m<n ? Q*back_substitute(R,b,transpose=true) :
-    back_substitute(R, transpose(Q)*b);
-
+    m<n 
+    // avoiding input validation in back_substitute
+    ? let( n  = len(R) )
+      Q*reverse(_back_substitute(transpose(R, reverse=true), reverse(b))) 
+    : _back_substitute(R, transpose(Q)*b);
 
 // Function: matrix_inverse()
 // Usage:
@@ -719,18 +715,6 @@ function matrix_inverse(A) =
     linear_solve(A,ident(len(A)));
 
 
-// Function: submatrix()
-// Usage: submatrix(M, ind1, ind2)
-// Description:
-//   Returns a submatrix with the specified index ranges or index sets.  
-function submatrix(M,ind1,ind2) =
-    assert( is_matrix(M), "Input must be a matrix." )
-    [for(i=ind1) 
-        [for(j=ind2) 
-            assert( ! is_undef(M[i][j]), "Invalid indexing." )
-            M[i][j] ] ];
-    
-
 // Function: qr_factor()
 // Usage: qr = qr_factor(A)
 // Description:
@@ -743,13 +727,14 @@ function qr_factor(A) =
       n = len(A[0])
     )
     let(
-        qr =_qr_factor(A, column=0, m = m, n=n, Q=ident(m)),
-        Rzero = [
-            for(i=[0:m-1]) [
-                for(j=[0:n-1])
-                i>j ? 0 : qr[1][i][j]
+        qr = _qr_factor(A, Q=ident(m), column=0, m = m, n=n),
+        Rzero = 
+          let( R = qr[1] )
+          [ for(i=[0:m-1]) [
+              let( ri = R[i] )
+              for(j=[0:n-1]) i>j ? 0 : ri[j]
             ]
-        ]
+          ]
     ) [qr[0],Rzero];
 
 function _qr_factor(A,Q, column, m, n) =
@@ -760,7 +745,13 @@ function _qr_factor(A,Q, column, m, n) =
         u = x - concat([alpha],repeat(0,m-1)),
         v = alpha==0 ? u : u / norm(u),
         Qc = ident(len(x)) - 2*outer_product(v,v),
-        Qf = [for(i=[0:m-1]) [for(j=[0:m-1]) i<column || j<column ? (i==j ? 1 : 0) : Qc[i-column][j-column]]]
+        Qf = [for(i=[0:m-1]) 
+                [for(j=[0:m-1]) 
+                    i<column || j<column 
+                    ? (i==j ? 1 : 0) 
+                    : Qc[i-column][j-column]
+                ]
+             ]
     )
     _qr_factor(Qf*A, Q*Qf, column+1, m, n);
 
@@ -773,26 +764,25 @@ function _qr_factor(A,Q, column, m, n) =
 //   You can supply a compatible matrix b and it will produce the solution for every column of b.  Note that if you want to
 //   solve Rx=b1 and Rx=b2 you must set b to transpose([b1,b2]) and then take the transpose of the result.  If the matrix
 //   is singular (e.g. has a zero on the diagonal) then it returns [].  
-function back_substitute(R, b, x=[],transpose = false) =
+function back_substitute(R, b, transpose = false) =
     assert(is_matrix(R, square=true))
     let(n=len(R))
     assert(is_vector(b,n) || is_matrix(b,n),str("R and b are not compatible in back_substitute ",n, len(b)))
-    !is_vector(b) ? transpose([for(i=[0:len(b[0])-1]) back_substitute(R,subindex(b,i),transpose=transpose)]) :
-    transpose?
-        reverse(back_substitute(
-            [for(i=[0:n-1]) [for(j=[0:n-1]) R[n-1-j][n-1-i]]],
-            reverse(b), x, false
-        )) :
-    len(x) == n ? x :
-    let(
-        ind = n - len(x) - 1
-    )
-    R[ind][ind] == 0 ? [] :
-    let(
-        newvalue =
-            len(x)==0? b[ind]/R[ind][ind] : 
-            (b[ind]-select(R[ind],ind+1,-1) * x)/R[ind][ind]
-    ) back_substitute(R, b, concat([newvalue],x));
+    transpose
+      ? reverse(_back_substitute(transpose(R, reverse=true), reverse(b)))  
+      : _back_substitute(R,b);
+
+function _back_substitute(R, b, x=[]) =
+    let(n=len(R))
+    len(x) == n ? x
+    : let(ind = n - len(x) - 1)
+      R[ind][ind] == 0 ? []
+    : let(
+          newvalue = len(x)==0
+            ? b[ind]/R[ind][ind]
+            : (b[ind]-select(R[ind],ind+1,-1) * x)/R[ind][ind]
+      )
+      _back_substitute(R, b, concat([newvalue],x));
 
 
 // Function: det2()
@@ -804,7 +794,7 @@ function back_substitute(R, b, x=[],transpose = false) =
 //   M = [ [6,-2], [1,8] ];
 //   det = det2(M);  // Returns: 50
 function det2(M) = 
-    assert( is_matrix(M,2,2), "Matrix should be 2x2." )
+    assert( 0*M==[[0,0],[0,0]], "Matrix should be 2x2." )
     M[0][0] * M[1][1] - M[0][1]*M[1][0];
 
 
@@ -817,7 +807,7 @@ function det2(M) =
 //   M = [ [6,4,-2], [1,-2,8], [1,5,7] ];
 //   det = det3(M);  // Returns: -334
 function det3(M) =
-    assert( is_matrix(M,3,3), "Matrix should be 3x3." )
+    assert( 0*M==[[0,0,0],[0,0,0],[0,0,0]], "Matrix should be 3x3." )
     M[0][0] * (M[1][1]*M[2][2]-M[2][1]*M[1][2]) -
     M[1][0] * (M[0][1]*M[2][2]-M[2][1]*M[0][2]) +
     M[2][0] * (M[0][1]*M[1][2]-M[1][1]*M[0][2]);
@@ -857,7 +847,7 @@ function determinant(M) =
 // Description:
 //   Returns true if A is a numeric matrix of height m and width n.  If m or n
 //   are omitted or set to undef then true is returned for any positive dimension.
-//   If `square` is true then the matrix is required to be square.  Note if you
+//   If `square` is true then the matrix is required to be square. 
 //   specify m != n and require a square matrix then the result will always be false.
 // Arguments:
 //   A = matrix to test
@@ -866,14 +856,131 @@ function determinant(M) =
 //   square = set to true to require a square matrix.  Default: false        
 function is_matrix(A,m,n,square=false) =
     is_list(A[0]) 
-    && ( let(v = A*A[0]) is_num(0*(v*v)) ) // a matrix of finite numbers 
-    && (is_undef(n) || len(A[0])==n )
-    && (is_undef(m) || len(A)==m )
-    && ( !square || len(A)==len(A[0]));
+    && ( let(v = A*A[0]) is_num(0*(v*v)) ) // a matrix of finite numbers 
+    && (is_undef(n) || len(A[0])==n )
+    && (is_undef(m) || len(A)==m )
+    && ( !square || len(A)==len(A[0]));
 
 
 // Section: Comparisons and Logic
 
+// Function: is_zero()
+// Usage:
+//   is_zero(x);
+// Description:
+//   Returns true if the number passed to it is approximately zero, to within `eps`.
+//   If passed a list, recursively checks if all items in the list are approximately zero.
+//   Otherwise, returns false.
+// Arguments:
+//   x = The value to check.
+//   eps = The maximum allowed variance.  Default: `EPSILON` (1e-9)
+// Example:
+//   is_zero(0);  // Returns: true.
+//   is_zero(1e-3);  // Returns: false.
+//   is_zero([0,0,0]);  // Returns: true.
+//   is_zero([0,0,1e-3]);  // Returns: false.
+function is_zero(x, eps=EPSILON) =
+    is_list(x)? (x != [] && [for (xx=x) if(!is_zero(xx,eps=eps)) 1] == []) :
+    is_num(x)? approx(x,eps) :
+    false;
+
+
+// Function: is_positive()
+// Usage:
+//   is_positive(x);
+// Description:
+//   Returns true if the number passed to it is greater than zero.
+//   If passed a list, recursively checks if all items in the list are positive.
+//   Otherwise, returns false.
+// Arguments:
+//   x = The value to check.
+// Example:
+//   is_positive(-2);  // Returns: false.
+//   is_positive(0);  // Returns: false.
+//   is_positive(2);  // Returns: true.
+//   is_positive([0,0,0]);  // Returns: false.
+//   is_positive([0,1,2]);  // Returns: false.
+//   is_positive([3,1,2]);  // Returns: true.
+//   is_positive([3,-1,2]);  // Returns: false.
+function is_positive(x) =
+    is_list(x)? (x != [] && [for (xx=x) if(!is_positive(xx)) 1] == []) :
+    is_num(x)? x>0 :
+    false;
+
+
+// Function: is_negative()
+// Usage:
+//   is_negative(x);
+// Description:
+//   Returns true if the number passed to it is less than zero.
+//   If passed a list, recursively checks if all items in the list are negative.
+//   Otherwise, returns false.
+// Arguments:
+//   x = The value to check.
+// Example:
+//   is_negative(-2);  // Returns: true.
+//   is_negative(0);  // Returns: false.
+//   is_negative(2);  // Returns: false.
+//   is_negative([0,0,0]);  // Returns: false.
+//   is_negative([0,1,2]);  // Returns: false.
+//   is_negative([3,1,2]);  // Returns: false.
+//   is_negative([3,-1,2]);  // Returns: false.
+//   is_negative([-3,-1,-2]);  // Returns: true.
+function is_negative(x) =
+    is_list(x)? (x != [] && [for (xx=x) if(!is_negative(xx)) 1] == []) :
+    is_num(x)? x<0 :
+    false;
+
+
+// Function: is_nonpositive()
+// Usage:
+//   is_nonpositive(x);
+// Description:
+//   Returns true if the number passed to it is less than or equal to zero.
+//   If passed a list, recursively checks if all items in the list are nonpositive.
+//   Otherwise, returns false.
+// Arguments:
+//   x = The value to check.
+// Example:
+//   is_nonpositive(-2);  // Returns: true.
+//   is_nonpositive(0);  // Returns: true.
+//   is_nonpositive(2);  // Returns: false.
+//   is_nonpositive([0,0,0]);  // Returns: true.
+//   is_nonpositive([0,1,2]);  // Returns: false.
+//   is_nonpositive([3,1,2]);  // Returns: false.
+//   is_nonpositive([3,-1,2]);  // Returns: false.
+//   is_nonpositive([-3,-1,-2]);  // Returns: true.
+function is_nonpositive(x) =
+    is_list(x)? (x != [] && [for (xx=x) if(!is_nonpositive(xx)) 1] == []) :
+    is_num(x)? x<=0 :
+    false;
+
+
+// Function: is_nonnegative()
+// Usage:
+//   is_nonnegative(x);
+// Description:
+//   Returns true if the number passed to it is greater than or equal to zero.
+//   If passed a list, recursively checks if all items in the list are nonnegative.
+//   Otherwise, returns false.
+// Arguments:
+//   x = The value to check.
+// Example:
+//   is_nonnegative(-2);  // Returns: false.
+//   is_nonnegative(0);  // Returns: true.
+//   is_nonnegative(2);  // Returns: true.
+//   is_nonnegative([0,0,0]);  // Returns: true.
+//   is_nonnegative([0,1,2]);  // Returns: true.
+//   is_nonnegative([0,-1,-2]);  // Returns: false.
+//   is_nonnegative([3,1,2]);  // Returns: true.
+//   is_nonnegative([3,-1,2]);  // Returns: false.
+//   is_nonnegative([-3,-1,-2]);  // Returns: false.
+function is_nonnegative(x) =
+    is_list(x)? (x != [] && [for (xx=x) if(!is_nonnegative(xx)) 1] == []) :
+    is_num(x)? x>=0 :
+    false;
+
+
 // Function: approx()
 // Usage:
 //   approx(a,b,[eps])
@@ -959,13 +1066,16 @@ function compare_lists(a, b) =
 //   any([1,5,true]);       // Returns true.
 //   any([[0,0], [0,0]]);   // Returns false.
 //   any([[0,0], [1,0]]);   // Returns true.
-function any(l, i=0, succ=false) =
-    (i>=len(l) || succ)? succ :
-    any( l, 
-         i+1, 
-         succ = is_list(l[i]) ? any(l[i]) : !(!l[i])
-        );
+function any(l) =
+    assert(is_list(l), "The input is not a list." )
+    _any(l, i=0, succ=false);
 
+function _any(l, i=0, succ=false) =
+    (i>=len(l) || succ)? succ :
+    _any( l, 
+         i+1, 
+         succ = is_list(l[i]) ? _any(l[i]) : !(!l[i])
+        );
 
 
 // Function: all()
@@ -982,12 +1092,15 @@ function any(l, i=0, succ=false) =
 //   all([[0,0], [1,0]]);   // Returns false.
 //   all([[1,1], [1,1]]);   // Returns true.
 function all(l, i=0, fail=false) =
-    (i>=len(l) || fail)? !fail :
-    all( l, 
-         i+1,
-         fail = is_list(l[i]) ? !all(l[i]) : !l[i]
-        ) ;
+    assert( is_list(l), "The input is not a list." )
+    _all(l, i=0, fail=false);
 
+function _all(l, i=0, fail=false) =
+    (i>=len(l) || fail)? !fail :
+    _all( l, 
+         i+1,
+         fail = is_list(l[i]) ? !_all(l[i]) : !l[i]
+        ) ;
 
 
 // Function: count_true()
@@ -1010,16 +1123,6 @@ function all(l, i=0, fail=false) =
 //   count_true([[0,0], [1,0]]);   // Returns 1.
 //   count_true([[1,1], [1,1]]);   // Returns 4.
 //   count_true([[1,1], [1,1]], nmax=3);  // Returns 3.
-function count_true(l, nmax=undef, i=0, cnt=0) =
-    (i>=len(l) || (nmax!=undef && cnt>=nmax))? cnt :
-    count_true(
-        l=l, nmax=nmax, i=i+1, cnt=cnt+(
-            is_list(l[i])? count_true(l[i], nmax=nmax-cnt) :
-            (l[i]? 1 : 0)
-        )
-    );
-
-
 function count_true(l, nmax) = 
     !is_list(l) ? !(!l) ? 1: 0 :
     let( c = [for( i = 0,
@@ -1120,19 +1223,21 @@ function _deriv_nonuniform(data, h, closed) =
 //   closed = boolean to indicate if the data set should be wrapped around from the end to the start.
 function deriv2(data, h=1, closed=false) =
     assert( is_consistent(data) , "Input list is not consistent or not numerical.") 
-    assert( len(data)>=3, "Input list has less than 3 elements.") 
     assert( is_finite(h), "The sampling `h` must be a number." )
     let( L = len(data) )
-    closed? [
+    assert( L>=3, "Input list has less than 3 elements.") 
+    closed
+    ? [
         for(i=[0:1:L-1])
         (data[(i+1)%L]-2*data[i]+data[(L+i-1)%L])/h/h
-    ] :
+      ]
+    :
     let(
-        first = L<3? undef : 
+        first = 
             L==3? data[0] - 2*data[1] + data[2] :
             L==4? 2*data[0] - 5*data[1] + 4*data[2] - data[3] :
             (35*data[0] - 104*data[1] + 114*data[2] - 56*data[3] + 11*data[4])/12, 
-        last = L<3? undef :
+        last = 
             L==3? data[L-1] - 2*data[L-2] + data[L-3] :
             L==4? -2*data[L-1] + 5*data[L-2] - 4*data[L-3] + data[L-4] :
             (35*data[L-1] - 104*data[L-2] + 114*data[L-3] - 56*data[L-4] + 11*data[L-5])/12
@@ -1212,34 +1317,13 @@ function C_div(z1,z2) =
 //   The polynomial is specified as p=[a_n, a_{n-1},...,a_1,a_0]
 //   where a_n is the z^n coefficient.  Polynomial coefficients are real.
 //   The result is a number if `z` is a number and a complex number otherwise.
-
-// Note: this should probably be recoded to use division by [1,-z], which is more accurate
-// and avoids overflow with large coefficients, but requires poly_div to support complex coefficients.  
-function polynomial(p, z, _k, _zk, _total) =
-    is_undef(_k)  
-    ?   assert( is_vector(p), "Input polynomial coefficients must be a vector." )
-        let(p = _poly_trim(p))
-        assert( is_finite(z) || is_vector(z,2), "The value of `z` must be a real or a complex number." )
-        polynomial( p, 
-                    z, 
-                    len(p)-1, 
-                    is_num(z)? 1 : [1,0], 
-                    is_num(z) ? 0 : [0,0]) 
-    :   _k==0 
-        ? _total + +_zk*p[0]
-        : polynomial( p, 
-                      z, 
-                      _k-1, 
-                      is_num(z) ? _zk*z : C_times(_zk,z), 
-                      _total+_zk*p[_k]);
-
 function polynomial(p,z,k,total) =
-     is_undef(k)
-   ?    assert( is_vector(p) , "Input polynomial coefficients must be a vector." )
+    is_undef(k)
+    ?   assert( is_vector(p) , "Input polynomial coefficients must be a vector." )
         assert( is_finite(z) || is_vector(z,2), "The value of `z` must be a real or a complex number." )
         polynomial( _poly_trim(p), z, 0, is_num(z) ? 0 : [0,0])
-   : k==len(p) ? total
-   : polynomial(p,z,k+1, is_num(z) ? total*z+p[k] : C_times(total,z)+[p[k],0]);
+    : k==len(p) ? total
+    : polynomial(p,z,k+1, is_num(z) ? total*z+p[k] : C_times(total,z)+[p[k],0]);
 
 // Function: poly_mult()
 // Usage:
@@ -1248,36 +1332,16 @@ function polynomial(p,z,k,total) =
 // Description:
 //   Given a list of polynomials represented as real coefficient lists, with the highest degree coefficient first, 
 //   computes the coefficient list of the product polynomial.  
-function poly_mult(p,q) = 
-    is_undef(q) ?
-       assert( is_list(p) 
-               && []==[for(pi=p) if( !is_vector(pi) && pi!=[]) 0], 
-               "Invalid arguments to poly_mult")
-       len(p)==2 ? poly_mult(p[0],p[1]) 
-                 : poly_mult(p[0], poly_mult(select(p,1,-1)))
-    :
-    _poly_trim(
-    [
-    for(n = [len(p)+len(q)-2:-1:0])
-        sum( [for(i=[0:1:len(p)-1])
-             let(j = len(p)+len(q)- 2 - n - i)
-             if (j>=0 && j<len(q)) p[i]*q[j]
-                 ])
-     ]);        
-     
 function poly_mult(p,q) = 
     is_undef(q) ?
-       len(p)==2 ? poly_mult(p[0],p[1]) 
-                 : poly_mult(p[0], poly_mult(select(p,1,-1)))
+        len(p)==2 
+        ? poly_mult(p[0],p[1]) 
+        : poly_mult(p[0], poly_mult(select(p,1,-1)))
     :
     assert( is_vector(p) && is_vector(q),"Invalid arguments to poly_mult")
-    _poly_trim( [
-                  for(n = [len(p)+len(q)-2:-1:0])
-                      sum( [for(i=[0:1:len(p)-1])
-                           let(j = len(p)+len(q)- 2 - n - i)
-                           if (j>=0 && j<len(q)) p[i]*q[j]
-                               ])
-                   ]);
+    p*p==0 || q*q==0
+    ? [0]
+    : _poly_trim(convolve(p,q));
 
     
 // Function: poly_div()
@@ -1288,19 +1352,23 @@ function poly_mult(p,q) =
 //    a list of two polynomials, [quotient, remainder].  If the division has no remainder then
 //    the zero polynomial [] is returned for the remainder.  Similarly if the quotient is zero
 //    the returned quotient will be [].  
-function poly_div(n,d,q) =
-    is_undef(q) 
-    ?   assert( is_vector(n) && is_vector(d) , "Invalid polynomials." )
-        let( d = _poly_trim(d) )
-        assert( d!=[0] , "Denominator cannot be a zero polynomial." )
-        poly_div(n,d,q=[])
-    :   len(n)<len(d) ? [q,_poly_trim(n)] : 
-        let(
-          t = n[0] / d[0], 
-          newq = concat(q,[t]),
-          newn = [for(i=[1:1:len(n)-1]) i<len(d) ? n[i] - t*d[i] : n[i]]
-        )  
-        poly_div(newn,d,newq);
+function poly_div(n,d) =
+    assert( is_vector(n) && is_vector(d) , "Invalid polynomials." )
+    let( d = _poly_trim(d), 
+         n = _poly_trim(n) )
+    assert( d!=[0] , "Denominator cannot be a zero polynomial." )
+    n==[0]
+    ? [[0],[0]]
+    : _poly_div(n,d,q=[]);
+
+function _poly_div(n,d,q) =
+    len(n)<len(d) ? [q,_poly_trim(n)] : 
+    let(
+      t = n[0] / d[0], 
+      newq = concat(q,[t]),
+      newn = [for(i=[1:1:len(n)-1]) i<len(d) ? n[i] - t*d[i] : n[i]]
+    )  
+    _poly_div(newn,d,newq);
 
 
 // Internal Function: _poly_trim()
diff --git a/partitions.scad b/partitions.scad
index 534219b..c46134c 100644
--- a/partitions.scad
+++ b/partitions.scad
@@ -31,9 +31,13 @@ _partition_cutpaths = [
 
 function _partition_cutpath(l, h, cutsize, cutpath, gap) =
     let(
+        check = assert(is_finite(l))
+            assert(is_finite(h))
+            assert(is_finite(gap))
+            assert(is_finite(cutsize) || is_vector(cutsize,2))
+            assert(is_string(cutpath) || is_path(cutpath,2)),
         cutsize = is_vector(cutsize)? cutsize : [cutsize*2, cutsize],
         cutpath = is_path(cutpath)? cutpath : (
-            assert(is_string(cutpath), "cutpath must be a 2D path or a string.")
             let(idx = search([cutpath], _partition_cutpaths))
             idx==[[]]? assert(in_list(cutpath,_partition_cutpaths,idx=0)) :
             _partition_cutpaths[idx.x][1]
@@ -79,7 +83,7 @@ function _partition_cutpath(l, h, cutsize, cutpath, gap) =
 //   partition_mask(w=20, cutpath="jigsaw");
 module partition_mask(l=100, w=100, h=100, cutsize=10, cutpath=undef, gap=0, inverse=false, spin=0, orient=UP)
 {
-    cutsize = is_vector(cutsize)? cutsize : [cutsize*2, cutsize];
+    cutsize = is_vector(cutsize)? point2d(cutsize) : [cutsize*2, cutsize];
     path = _partition_cutpath(l, h, cutsize, cutpath, gap);
     fullpath = concat(path, [[l/2,w*(inverse?-1:1)], [-l/2,w*(inverse?-1:1)]]);
     rot(from=UP,to=orient) {
diff --git a/paths.scad b/paths.scad
index 7d5f978..93d82fa 100644
--- a/paths.scad
+++ b/paths.scad
@@ -43,10 +43,12 @@ include <triangulation.scad>
 //   dim = list of allowed dimensions of the vectors in the path.  Default: [2,3]
 //   fast = set to true for fast check that only looks at first entry.  Default: false
 function is_path(list, dim=[2,3], fast=false) =
-    fast? is_list(list) && is_vector(list[0]) :
-    is_list(list) && is_list(list[0]) && len(list)>1 &&
-    (is_undef(dim) || in_list(len(list[0]), force_list(dim))) &&
-    is_list_of(list, repeat(0,len(list[0])));
+    fast
+    ?   is_list(list) && is_vector(list[0]) 
+    :   is_matrix(list) 
+        && len(list)>1 
+        && len(list[0])>0
+        && (is_undef(dim) || in_list(len(list[0]), force_list(dim)));
 
 
 // Function: is_closed_path()
@@ -105,32 +107,51 @@ function path_subselect(path, s1, u1, s2, u2, closed=false) =
 
 // Function: simplify_path()
 // Description:
-//   Takes a path and removes unnecessary collinear points.
+//   Takes a path and removes unnecessary subsequent collinear points.
 // Usage:
 //   simplify_path(path, [eps])
 // Arguments:
-//   path = A list of 2D path points.
+//   path = A list of path points of any dimension.
 //   eps = Largest positional variance allowed.  Default: `EPSILON` (1-e9)
 function simplify_path(path, eps=EPSILON) =
-    len(path)<=2? path : let(
-        indices = concat([0], [for (i=[1:1:len(path)-2]) if (!collinear_indexed(path, i-1, i, i+1, eps=eps)) i], [len(path)-1])
-    ) [for (i = indices) path[i]];
+    assert( is_path(path), "Invalid path." )
+    assert( is_undef(eps) || (is_finite(eps) && (eps>=0) ), "Invalid tolerance." )    
+    len(path)<=2 ? path 
+    :   let(
+            indices = [ 0,
+                        for (i=[1:1:len(path)-2]) 
+                            if (!collinear(path[i-1],path[i],path[i+1], eps=eps)) i, 
+                        len(path)-1 
+                      ]
+        ) 
+        [for (i = indices) path[i] ];
 
 
 // Function: simplify_path_indexed()
 // Description:
-//   Takes a list of points, and a path as a list of indices into `points`,
-//   and removes all path points that are unecessarily collinear.
+//   Takes a list of points, and a list of indices into `points`,
+//   and removes from the list all indices of subsequent indexed points that are unecessarily collinear.
+//   Returns the list of the remained indices.
 // Usage:
-//   simplify_path_indexed(path, eps)
+//   simplify_path_indexed(points,indices, eps)
 // Arguments:
 //   points = A list of points.
-//   path = A list of indices into `points` that forms a path.
+//   indices = A list of indices into `points` that forms a path.
 //   eps = Largest angle variance allowed.  Default: EPSILON (1-e9) degrees.
-function simplify_path_indexed(points, path, eps=EPSILON) =
-    len(path)<=2? path : let(
-        indices = concat([0], [for (i=[1:1:len(path)-2]) if (!collinear_indexed(points, path[i-1], path[i], path[i+1], eps=eps)) i], [len(path)-1])
-    ) [for (i = indices) path[i]];
+function simplify_path_indexed(points, indices, eps=EPSILON) =
+    len(indices)<=2? indices 
+    : let(
+        indices = concat( indices[0], 
+                          [for (i=[1:1:len(indices)-2]) 
+                              let( 
+                                  i1 = indices[i-1],
+                                  i2 = indices[i],
+                                  i3 = indices[i+1]
+                              )
+                              if (!collinear(points[i1],points[i2],points[i3], eps=eps)) indices[i]], 
+                          indices[len(indices)-1] )
+        ) 
+        indices;
 
 
 // Function: path_length()
diff --git a/polyhedra.scad b/polyhedra.scad
index b199fe9..e845acd 100644
--- a/polyhedra.scad
+++ b/polyhedra.scad
@@ -710,7 +710,7 @@ function regular_polyhedron_info(
     info == "center" ? translation :
     info == "type" ? entry[class] :
     info == "name" ? entry[pname] :
-    echo_warning(str("Unknown info type '",info,"' requested"));
+    assert(false, str("Unknown info type '",info,"' requested"));
 
 
 
diff --git a/regions.scad b/regions.scad
index 6169337..fa761f0 100644
--- a/regions.scad
+++ b/regions.scad
@@ -339,7 +339,7 @@ function region_faces(region, transform, reverse=false, vnf=EMPTY_VNF) =
 //   linear_sweep(orgn,height=20,convexity=16) show_anchors();
 module linear_sweep(region, height=1, center, twist=0, scale=1, slices, maxseg, style="default", convexity, anchor_isect=false, anchor, spin=0, orient=UP) {
     region = is_path(region)? [region] : region;
-    cp = median(flatten(region));
+    cp = mean(pointlist_bounds(flatten(region)));
     anchor = get_anchor(anchor, center, "origin", "origin");
     vnf = linear_sweep(
         region, height=height,
diff --git a/rounding.scad b/rounding.scad
index 254d035..ae2aed0 100644
--- a/rounding.scad
+++ b/rounding.scad
@@ -458,10 +458,10 @@ function smooth_path(path, tangents, size, relsize, splinesteps=10, uniform=fals
 
 
 
-// Module: offset_sweep()
+// Function&Module: offset_sweep()
 //
 // Description:
-//   Takes a 2d path as input and extrudes it upwards and/or downward.  Each layer in the extrusion is produced using `offset()` to expand or shrink the previous layer.
+//   Takes a 2d path as input and extrudes it upwards and/or downward.  Each layer in the extrusion is produced using `offset()` to expand or shrink the previous layer.  When invoked as a function returns a VNF; when invoked as a module produces geometry.  
 //   You can specify a sequence of offsets values, or you can use several built-in offset profiles that are designed to provide end treatments such as roundovers.
 //   The path is shifted by `offset()` multiple times in sequence
 //   to produce the final shape (not multiple shifts from one parent), so coarse definition of the input path will degrade
@@ -543,8 +543,12 @@ function smooth_path(path, tangents, size, relsize, splinesteps=10, uniform=fals
 //   angle = default angle for chamfers.  Default: 45
 //   joint = default joint value for smooth roundover.
 //   k = default curvature parameter value for "smooth" roundover
-//   convexity = convexity setting for use with polyhedron.  Default: 10
-//
+//   convexity = convexity setting for use with polyhedron.  (module only) Default: 10
+//   anchor = Translate so anchor point is at the origin.  (module only) Default: "origin"
+//   spin = Rotate this many degrees around Z axis after anchor.  (module only) Default: 0
+//   orient = Vector to rotate top towards after spin  (module only)
+//   extent = use extent method for computing anchors. (module only)  Default: false
+//   cp = set centerpoint for anchor computation.  (module only) Default: object centroid
 // Example: Rounding a star shaped prism with postive radius values
 //   star = star(5, r=22, ir=13);
 //   rounded_star = round_corners(star, cut=flatten(repeat([.5,0],5)), $fn=24);
@@ -650,118 +654,118 @@ function smooth_path(path, tangents, size, relsize, splinesteps=10, uniform=fals
 //       up(1)
 //         offset_sweep(offset(rhex,r=-1), height=9.5, bottom=os_circle(r=2), top=os_teardrop(r=-4));
 //     }
-module offset_sweep(
-        path, height, h, l,
-        top=[], bottom=[],
-        offset="round", r=0, steps=16,
-        quality=1, check_valid=true,
-        offset_maxstep=1, extra=0,
-        cut=undef, chamfer_width=undef, chamfer_height=undef,
-        joint=undef, k=0.75, angle=45,
-        convexity=10
-) {
-        // This function does the actual work of repeatedly calling offset() and concatenating the resulting face and vertex lists to produce
-        // the inputs for the polyhedron module.
-        function make_polyhedron(path,offsets, offset_type, flip_faces, quality, check_valid, maxstep, offsetind=0, vertexcount=0, vertices=[], faces=[] )=
-                offsetind==len(offsets)? (
-                        let(
-                                bottom = list_range(n=len(path),s=vertexcount),
-                                oriented_bottom = !flip_faces? bottom : reverse(bottom)
-                        ) [vertices, concat(faces,[oriented_bottom])]
-                ) : (
-                        let(
-                                this_offset = offsetind==0? offsets[0][0] : offsets[offsetind][0] - offsets[offsetind-1][0],
-                                delta = offset_type=="delta" || offset_type=="chamfer" ? this_offset : undef,
-                                r = offset_type=="round"? this_offset : undef,
-                                do_chamfer = offset_type == "chamfer"
+
+
+// This function does the actual work of repeatedly calling offset() and concatenating the resulting face and vertex lists to produce
+// the inputs for the polyhedron module.
+function _make_offset_polyhedron(path,offsets, offset_type, flip_faces, quality, check_valid, maxstep, offsetind=0,
+                                 vertexcount=0, vertices=[], faces=[] )=
+        offsetind==len(offsets)? (
+                let(
+                        bottom = list_range(n=len(path),s=vertexcount),
+                        oriented_bottom = !flip_faces? bottom : reverse(bottom)
+                ) [vertices, concat(faces,[oriented_bottom])]
+        ) : (
+                let(
+                        this_offset = offsetind==0? offsets[0][0] : offsets[offsetind][0] - offsets[offsetind-1][0],
+                        delta = offset_type=="delta" || offset_type=="chamfer" ? this_offset : undef,
+                        r = offset_type=="round"? this_offset : undef,
+                        do_chamfer = offset_type == "chamfer"
+                )
+                let(
+                        vertices_faces = offset(
+                                path, r=r, delta=delta, chamfer = do_chamfer, closed=true,
+                                check_valid=check_valid, quality=quality,
+                                maxstep=maxstep, return_faces=true,
+                                firstface_index=vertexcount,
+                                flip_faces=flip_faces
                         )
-                        assert(num_defined([r,delta])==1,"Must set `offset` to \"round\" or \"delta")
-                        let(
-                                vertices_faces = offset(
-                                        path, r=r, delta=delta, chamfer = do_chamfer, closed=true,
-                                        check_valid=check_valid, quality=quality,
-                                        maxstep=maxstep, return_faces=true,
-                                        firstface_index=vertexcount,
-                                        flip_faces=flip_faces
-                                )
-                        )
-                        make_polyhedron(
-                                vertices_faces[0], offsets, offset_type,
-                                flip_faces, quality, check_valid, maxstep,
-                                offsetind+1, vertexcount+len(path),
-                                vertices=concat(
-                                        vertices,
-                                        zip(vertices_faces[0],repeat(offsets[offsetind][1],len(vertices_faces[0])))
-                                ),
-                                faces=concat(faces, vertices_faces[1])
-                        )
-                );
-
-
-        argspec = [
-                ["r",r],
-                ["extra",extra],
-                ["type","circle"],
-                ["check_valid",check_valid],
-                ["quality",quality],
-                ["offset_maxstep", offset_maxstep],
-                ["steps",steps],
-                ["offset",offset],
-                ["chamfer_width",chamfer_width],
-                ["chamfer_height",chamfer_height],
-                ["angle",angle],
-                ["cut",cut],
-                ["joint",joint],
-                ["k", k],
-                ["points", []],
-        ];
-
-        path = check_and_fix_path(path, [2], closed=true);
-        clockwise = polygon_is_clockwise(path);
-
-        top = struct_set(argspec, top, grow=false);
-        bottom = struct_set(argspec, bottom, grow=false);
-
-        //  This code does not work.  It hits the error in make_polyhedron from offset being wrong
-        //  before this code executes.  Had to move the test into make_polyhedron, which is ugly since it's in the loop
-        //offsetsok = in_list(struct_val(top, "offset"),["round","delta"]) &&
-        //      in_list(struct_val(bottom, "offset"),["round","delta"]);
-        //assert(offsetsok,"Offsets must be one of \"round\" or \"delta\"");
-
-
-        offsets_bot = _rounding_offsets(bottom, -1);
-        offsets_top = _rounding_offsets(top, 1);
-
-        if (offset == "chamfer" && (len(offsets_bot)>5 || len(offsets_top)>5)) {
-          echo("WARNING: You have selected offset=\"chamfer\", which leads to exponential growth in the vertex count and requested many layers.  This can be slow or run out of recursion depth.");
-        }
-        // "Extra" height enlarges the result beyond the requested height, so subtract it
-        bottom_height = len(offsets_bot)==0 ? 0 : abs(select(offsets_bot,-1)[1]) - struct_val(bottom,"extra");
-        top_height = len(offsets_top)==0 ? 0 : abs(select(offsets_top,-1)[1]) - struct_val(top,"extra");
-
-        height = get_height(l=l,h=h,height=height,dflt=bottom_height+top_height);
-        assert(height>=0, "Height must be nonnegative");
-
-        middle = height-bottom_height-top_height;
-        assert(
-                middle>=0, str(
-                        "Specified end treatments (bottom height = ",bottom_height,
-                        " top_height = ",top_height,") are too large for extrusion height (",height,")"
+                )
+                _make_offset_polyhedron(
+                        vertices_faces[0], offsets, offset_type,
+                        flip_faces, quality, check_valid, maxstep,
+                        offsetind+1, vertexcount+len(path),
+                        vertices=concat(
+                                vertices,
+                                zip(vertices_faces[0],repeat(offsets[offsetind][1],len(vertices_faces[0])))
+                        ),
+                        faces=concat(faces, vertices_faces[1])
                 )
         );
-        initial_vertices_bot = path3d(path);
 
-        vertices_faces_bot = make_polyhedron(
+
+function offset_sweep(
+                       path, height, h, l,
+                       top=[], bottom=[],
+                       offset="round", r=0, steps=16,
+                       quality=1, check_valid=true,
+                       offset_maxstep=1, extra=0,
+                       cut=undef, chamfer_width=undef, chamfer_height=undef,
+                       joint=undef, k=0.75, angle=45
+                      ) =
+    let(
+        argspec = [
+                   ["r",r],
+                   ["extra",extra],
+                   ["type","circle"],
+                   ["check_valid",check_valid],
+                   ["quality",quality],
+                   ["offset_maxstep", offset_maxstep],
+                   ["steps",steps],
+                   ["offset",offset],
+                   ["chamfer_width",chamfer_width],
+                   ["chamfer_height",chamfer_height],
+                   ["angle",angle],
+                   ["cut",cut],
+                   ["joint",joint],
+                   ["k", k],
+                   ["points", []],
+        ],
+        path = check_and_fix_path(path, [2], closed=true),
+        clockwise = polygon_is_clockwise(path),
+        
+        top = struct_set(argspec, top, grow=false),
+        bottom = struct_set(argspec, bottom, grow=false),
+
+        //  This code does not work.  It hits the error in _make_offset_polyhedron from offset being wrong
+        //  before this code executes.  Had to move the test into _make_offset_polyhedron, which is ugly since it's in the loop
+        offsetsok = in_list(struct_val(top, "offset"),["round","delta"])
+                    && in_list(struct_val(bottom, "offset"),["round","delta"])
+    )
+    assert(offsetsok,"Offsets must be one of \"round\" or \"delta\"")
+    let( 
+        offsets_bot = _rounding_offsets(bottom, -1),
+        offsets_top = _rounding_offsets(top, 1),
+        dummy = offset == "chamfer" && (len(offsets_bot)>5 || len(offsets_top)>5)
+                ? echo("WARNING: You have selected offset=\"chamfer\", which leads to exponential growth in the vertex count and requested more than 5 layers.  This can be slow or run out of recursion depth.")
+                : 0,
+
+        // "Extra" height enlarges the result beyond the requested height, so subtract it
+        bottom_height = len(offsets_bot)==0 ? 0 : abs(select(offsets_bot,-1)[1]) - struct_val(bottom,"extra"),
+        top_height = len(offsets_top)==0 ? 0 : abs(select(offsets_top,-1)[1]) - struct_val(top,"extra"),
+
+        height = get_height(l=l,h=h,height=height,dflt=bottom_height+top_height),
+        middle = height-bottom_height-top_height
+    )
+    assert(height>=0, "Height must be nonnegative") 
+    assert(middle>=0, str("Specified end treatments (bottom height = ",bottom_height,
+                          " top_height = ",top_height,") are too large for extrusion height (",height,")"
+                         )
+    )
+    let(
+        initial_vertices_bot = path3d(path),
+
+        vertices_faces_bot = _make_offset_polyhedron(
                 path, offsets_bot, struct_val(bottom,"offset"), clockwise,
                 struct_val(bottom,"quality"),
                 struct_val(bottom,"check_valid"),
                 struct_val(bottom,"offset_maxstep"),
                 vertices=initial_vertices_bot
-        );
+        ),
 
-        top_start_ind = len(vertices_faces_bot[0]);
-        initial_vertices_top = zip(path, repeat(middle,len(path)));
-        vertices_faces_top = make_polyhedron(
+        top_start_ind = len(vertices_faces_bot[0]),
+        initial_vertices_top = zip(path, repeat(middle,len(path))),
+        vertices_faces_top = _make_offset_polyhedron(
                 path, move(p=offsets_top,[0,middle]),
                 struct_val(top,"offset"), !clockwise,
                 struct_val(top,"quality"),
@@ -769,20 +773,39 @@ module offset_sweep(
                 struct_val(top,"offset_maxstep"),
                 vertexcount=top_start_ind,
                 vertices=initial_vertices_top
-        );
+        ),
         middle_faces = middle==0 ? [] : [
                 for(i=[0:len(path)-1]) let(
                         oneface=[i, (i+1)%len(path), top_start_ind+(i+1)%len(path), top_start_ind+i]
                 ) !clockwise ? reverse(oneface) : oneface
-        ];
-        up(bottom_height) {
-                polyhedron(
-                        concat(vertices_faces_bot[0],vertices_faces_top[0]),
-                        faces=concat(vertices_faces_bot[1], vertices_faces_top[1], middle_faces),
-                        convexity=convexity
-                );
-        }
-}
+        ]
+    )
+    [up(bottom_height, concat(vertices_faces_bot[0],vertices_faces_top[0])),  // Vertices
+     concat(vertices_faces_bot[1], vertices_faces_top[1], middle_faces)];     // Faces
+
+
+module offset_sweep(path, height, h, l,
+                    top=[], bottom=[],
+                    offset="round", r=0, steps=16,
+                    quality=1, check_valid=true,
+                    offset_maxstep=1, extra=0,
+                    cut=undef, chamfer_width=undef, chamfer_height=undef,
+                    joint=undef, k=0.75, angle=45,
+                    convexity=10,anchor="origin",cp,
+                    spin=0, orient=UP, extent=false)
+{
+    vnf = offset_sweep(path=path, height=height, h=h, l=l, top=top, bottom=bottom, offset=offset, r=0, steps=steps,
+                       quality=quality, check_valid=true, offset_maxstep=1, extra=0, cut=cut, chamfer_width=chamfer_width,
+                       chamfer_height=chamfer_height, joint=joint, k=k, angle=angle);
+  
+    attachable(anchor=anchor, spin=spin, orient=orient, vnf=vnf, extent=extent, cp=is_def(cp) ? cp : vnf_centroid(vnf))
+    {
+        vnf_polyhedron(vnf,convexity=convexity);
+        children();
+    }   
+}   
+
+
 
 function os_circle(r,cut,extra,check_valid, quality,steps, offset_maxstep, offset) =
         assert(num_defined([r,cut])==1, "Must define exactly one of `r` and `cut`")
@@ -924,7 +947,6 @@ function os_profile(points, extra,check_valid, quality, offset_maxstep, offset)
 //   joint = default joint value for smooth roundover.
 //   k = default curvature parameter value for "smooth" roundover
 //   convexity = convexity setting for use with polyhedron.  Default: 10
-//
 // Example: Chamfered elliptical prism.  If you stretch a chamfered cylinder the chamfer will be uneven.
 //   convex_offset_extrude(bottom = os_chamfer(height=-2), top=os_chamfer(height=1), height=7)
 //   xscale(4)circle(r=6,$fn=64);
@@ -1364,8 +1386,6 @@ module offset_stroke(path, width=1, rounded=true, start, end, check_valid=true,
         }
 }
 
-
-
 function _rp_compute_patches(top, bot, rtop, rsides, ktop, ksides, concave) =
    let(
      N = len(top),
@@ -1395,8 +1415,8 @@ function _rp_compute_patches(top, bot, rtop, rsides, ktop, ksides, concave) =
                     let(
                        prev_corner = prev_offset + abs(rtop_in)*in_prev,
                        next_corner = next_offset + abs(rtop_in)*in_next,
-                       prev_degenerate = is_undef(ray_intersection([far_corner, far_corner+prev], [prev_offset, prev_offset+in_prev])),
-                       next_degenerate = is_undef(ray_intersection([far_corner, far_corner+next], [next_offset, next_offset+in_next]))
+                       prev_degenerate = is_undef(ray_intersection(path2d([far_corner, far_corner+prev]), path2d([prev_offset, prev_offset+in_prev]))),
+                       next_degenerate = is_undef(ray_intersection(path2d([far_corner, far_corner+next]), path2d([next_offset, next_offset+in_next])))
                     )
                     [ prev_degenerate ? far_corner : prev_corner,
                       far_corner,
@@ -1452,6 +1472,11 @@ function _rp_compute_patches(top, bot, rtop, rsides, ktop, ksides, concave) =
 //   splinesteps = number of segments to use for curved patches.  Default: 16
 //   debug = turn on debug mode which displays illegal polyhedra and shows the bezier corner patches for troubleshooting purposes.  Default: False
 //   convexity = convexity parameter for polyhedron(), only for module version.  Default: 10
+//   anchor = Translate so anchor point is at the origin.  (module only) Default: "origin"
+//   spin = Rotate this many degrees around Z axis after anchor.  (module only) Default: 0
+//   orient = Vector to rotate top towards after spin  (module only)
+//   extent = use extent method for computing anchors. (module only)  Default: false
+//   cp = set centerpoint for anchor computation.  (module only) Default: object centroid
 // Example: Uniformly rounded pentagonal prism
 //   rounded_prism(pentagon(3), height=3, joint_top=0.5, joint_bot=0.5, joint_sides=0.5);
 // Example: Maximum possible rounding.
@@ -1500,15 +1525,21 @@ function _rp_compute_patches(top, bot, rtop, rsides, ktop, ksides, concave) =
 //   rounded_prism(apply(yrot(95),path3d(hexagon(3))), apply(yrot(95), path3d(hexagon(3),3)), joint_top=2, joint_bot=1, joint_sides=1);
 
 module rounded_prism(bottom, top, joint_bot, joint_top, joint_sides, k_bot, k_top, k_sides,
-                     k=0.5, splinesteps=16, h, length, l, height, convexity=10, debug=false)
+                     k=0.5, splinesteps=16, h, length, l, height, convexity=10, debug=false,
+                     anchor="origin",cp,spin=0, orient=UP, extent=false)
 {
   result = rounded_prism(bottom=bottom, top=top, joint_bot=joint_bot, joint_top=joint_top, joint_sides=joint_sides,
                          k_bot=k_bot, k_top=k_top, k_sides=k_sides, k=k, splinesteps=splinesteps, h=h, length=length, height=height, l=l,debug=debug);
-  if (debug){
-      vnf_polyhedron(result[1], convexity=convexity);
-      trace_bezier_patches(result[0], showcps=true, splinesteps=splinesteps, $fn=16, showdots=false, showpatch=false);
+  vnf = debug ? result[1] : result;
+  attachable(anchor=anchor, spin=spin, orient=orient, vnf=vnf, extent=extent, cp=is_def(cp) ? cp : vnf_centroid(vnf))
+  {
+    if (debug){
+        vnf_polyhedron(vnf, convexity=convexity);
+        trace_bezier_patches(result[0], showcps=true, splinesteps=splinesteps, $fn=16, showdots=false, showpatch=false);
+    }
+    else vnf_polyhedron(vnf,convexity=convexity);
+    children();
   }
-  else vnf_polyhedron(result,convexity=convexity);
 }
 
 
@@ -1532,7 +1563,7 @@ function rounded_prism(bottom, top, joint_bot, joint_top, joint_sides, k_bot, k_
      // Determine which points are concave by making bottom 2d if necessary
      bot_proj = len(bottom[0])==2 ? bottom :  project_plane(bottom, select(bottom,0,2)),
      bottom_sign = polygon_is_clockwise(bot_proj) ? 1 : -1,
-     concave = [for(i=[0:N-1]) bottom_sign*sign(point_left_of_segment2d(select(bot_proj,i+1), select(bot_proj, i-1,i)))>0],
+     concave = [for(i=[0:N-1]) bottom_sign*sign(point_left_of_line2d(select(bot_proj,i+1), select(bot_proj, i-1,i)))>0],
      top = is_undef(top) ? path3d(bottom,height/2) :
            len(top[0])==2 ? path3d(top,height/2) :
            top,
@@ -1547,16 +1578,16 @@ function rounded_prism(bottom, top, joint_bot, joint_top, joint_sides, k_bot, k_
    assert(jsvecok || jssingleok,
           str("Argument joint_sides is invalid.  All entries must be nonnegative, and it must be a number, 2-vector, or a length ",N," list those."))
    assert(is_num(k_sides) || is_vector(k_sides,N), str("Curvature parameter k_sides must be a number or length ",N," vector"))
-   assert(points_are_coplanar(bottom))
-   assert(points_are_coplanar(top))
+   assert(coplanar(bottom))
+   assert(coplanar(top))
    assert(!is_num(k_sides) || (k_sides>=0 && k_sides<=1), "Curvature parameter k_sides must be in interval [0,1]")
    let(
-     non_coplanar=[for(i=[0:N-1]) if (!points_are_coplanar(concat(select(top,i,i+1), select(bottom,i,i+1)))) [i,(i+1)%N]],
+     non_coplanar=[for(i=[0:N-1]) if (!coplanar(concat(select(top,i,i+1), select(bottom,i,i+1)))) [i,(i+1)%N]],
      k_sides_vec = is_num(k_sides) ? repeat(k_sides, N) : k_sides,
      kbad = [for(i=[0:N-1]) if (k_sides_vec[i]<0 || k_sides_vec[i]>1) i],
      joint_sides_vec = jssingleok ? repeat(joint_sides,N) : joint_sides,
-     top_collinear = [for(i=[0:N-1]) if (points_are_collinear(select(top,i-1,i+1))) i],
-     bot_collinear = [for(i=[0:N-1]) if (points_are_collinear(select(bottom,i-1,i+1))) i]
+     top_collinear = [for(i=[0:N-1]) if (collinear(select(top,i-1,i+1))) i],
+     bot_collinear = [for(i=[0:N-1]) if (collinear(select(bottom,i-1,i+1))) i]
    )
    assert(non_coplanar==[], str("Side faces are non-coplanar at edges: ",non_coplanar))
    assert(top_collinear==[], str("Top has collinear or duplicated points at indices: ",top_collinear))
@@ -1622,14 +1653,14 @@ function rounded_prism(bottom, top, joint_bot, joint_top, joint_sides, k_bot, k_
                vline = concat(select(subindex(top_patch[i],j),2,4),
                               select(subindex(bot_patch[i],j),2,4))
              )
-         if (!points_are_collinear(vline)) [i,j]],
+         if (!collinear(vline)) [i,j]],
      //verify horiz edges
      verify_horiz=[for(i=[0:N-1], j=[0:4])
          let(
              hline_top = concat(select(top_patch[i][j],2,4), select(select(top_patch, i+1)[j],0,2)),
              hline_bot = concat(select(bot_patch[i][j],2,4), select(select(bot_patch, i+1)[j],0,2))
          )
-         if (!points_are_collinear(hline_top) || !points_are_collinear(hline_bot)) [i,j]]
+         if (!collinear(hline_top) || !collinear(hline_bot)) [i,j]]
     )
     assert(debug || top_intersections==[],
           "Roundovers interfere with each other on top face: either input is self intersecting or top joint length is too large")
@@ -1880,7 +1911,7 @@ function _circle_mask(r) =
 //     $fn=128;
 //     difference(){
 //       tube(or=r, wall=2, h=45);
-//       bent_cutout_mask(r-1, 2.1, apply(back(15),subdivide_path(round_corners(star(n=7,ir=5,or=10), cut=flatten(repeat([0.5,0],7))),14*15,closed=true)));
+//       bent_cutout_mask(r-1, 2.1, apply(back(15),subdivide_path(round_corners(star(n=7,ir=5,or=10), cut=flatten(repeat([0.5,0],7)),$fn=32),14*15,closed=true)));
 //     }
 //   }
 // Example(2D): Cutting a slot in a cylinder is tricky if you want rounded corners at the top.  This slot profile has slightly angled top edges to blend into the top edge of the cylinder.
@@ -1944,6 +1975,7 @@ function _circle_mask(r) =
 
 module bent_cutout_mask(r, thickness, path, convexity=10)
 {
+  no_children($children);
   assert(is_path(path,2),"Input path must be a 2d path")
   assert(r-thickness>0, "Thickness too large for radius");
   assert(thickness>0, "Thickness must be positive");
@@ -1962,4 +1994,4 @@ module bent_cutout_mask(r, thickness, path, convexity=10)
 }
 
 
-// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap
+// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap
\ No newline at end of file
diff --git a/scripts/make_all_docs.sh b/scripts/make_all_docs.sh
index ddd9937..fa05079 100755
--- a/scripts/make_all_docs.sh
+++ b/scripts/make_all_docs.sh
@@ -19,7 +19,7 @@ done
 if [[ "$FILES" != "" ]]; then
     PREVIEW_LIBS="$FILES"
 else
-    PREVIEW_LIBS="affine arrays attachments beziers bottlecaps common constants coords cubetruss debug distributors edges errors geometry hingesnaps hull involute_gears joiners knurling linear_bearings masks math metric_screws mutators nema_steppers partitions paths phillips_drive polyhedra primitives quaternions queues regions rounding screws shapes shapes2d skin sliders stacks strings structs threading torx_drive transforms triangulation vectors version vnf walls wiring"
+    PREVIEW_LIBS="affine arrays attachments beziers bottlecaps common constants coords cubetruss debug distributors edges geometry hingesnaps hull involute_gears joiners knurling linear_bearings masks math metric_screws mutators nema_steppers partitions paths phillips_drive polyhedra primitives quaternions queues regions rounding screws shapes shapes2d skin sliders stacks strings structs threading torx_drive transforms triangulation vectors version vnf walls wiring"
 fi
 
 dir="$(basename $PWD)"
diff --git a/shapes.scad b/shapes.scad
index cb6d4e5..d612b0a 100644
--- a/shapes.scad
+++ b/shapes.scad
@@ -102,8 +102,8 @@ module cuboid(
                 if (edges == EDGES_ALL && trimcorners) {
                     if (chamfer<0) {
                         cube(size, center=true) {
-                            attach(TOP) prismoid([size.x,size.y], [size.x-2*chamfer,size.y-2*chamfer], h=-chamfer, anchor=TOP);
-                            attach(BOT) prismoid([size.x,size.y], [size.x-2*chamfer,size.y-2*chamfer], h=-chamfer, anchor=TOP);
+                            attach(TOP,overlap=0) prismoid([size.x,size.y], [size.x-2*chamfer,size.y-2*chamfer], h=-chamfer, anchor=TOP);
+                            attach(BOT,overlap=0) prismoid([size.x,size.y], [size.x-2*chamfer,size.y-2*chamfer], h=-chamfer, anchor=TOP);
                         }
                     } else {
                         isize = [for (v = size) max(0.001, v-2*chamfer)];
diff --git a/shapes2d.scad b/shapes2d.scad
index 4c13d91..4c11efd 100644
--- a/shapes2d.scad
+++ b/shapes2d.scad
@@ -785,7 +785,7 @@ function rect(size=1, center, rounding=0, chamfer=0, anchor, spin=0) =
     assert(is_num(rounding) || len(rounding)==4)
     let(
         size = is_num(size)? [size,size] : point2d(size),
-        anchor = get_anchor(anchor, center, FRONT+LEFT, FRONT+LEFT),
+        anchor = point2d(get_anchor(anchor, center, FRONT+LEFT, FRONT+LEFT)),
         complex = rounding!=0 || chamfer!=0
     )
     (rounding==0 && chamfer==0)? let(
diff --git a/skin.scad b/skin.scad
index 920f4ed..87ae8d4 100644
--- a/skin.scad
+++ b/skin.scad
@@ -16,7 +16,8 @@ include <vnf.scad>
 
 // Function&Module: skin()
 // Usage: As module:
-//   skin(profiles, [slices], [refine], [method], [sampling], [caps], [closed], [z]);
+//   skin(profiles, [slices], [refine], [method], [sampling], [caps], [closed], [z], [convexity],
+//        [anchor],[cp],[spin],[orient],[extent]);
 // Usage: As function:
 //   vnf = skin(profiles, [slices], [refine], [method], [sampling], [caps], [closed], [z]);
 // Description:
@@ -117,6 +118,12 @@ include <vnf.scad>
 //   caps = true to create endcap faces when closed is false.  Can be a length 2 boolean array.  Default is true if closed is false.
 //   method = method for connecting profiles, one of "distance", "tangent", "direct" or "reindex".  Default: "direct".
 //   z = array of height values for each profile if the profiles are 2d
+//   convexity = convexity setting for use with polyhedron.  (module only) Default: 10
+//   anchor = Translate so anchor point is at the origin.  (module only) Default: "origin"
+//   spin = Rotate this many degrees around Z axis after anchor.  (module only) Default: 0
+//   orient = Vector to rotate top towards after spin  (module only)
+//   extent = use extent method for computing anchors. (module only)  Default: false
+//   cp = set centerpoint for anchor computation.  (module only) Default: object centroid
 // Example:
 //   skin([octagon(4), circle($fn=70,r=2)], z=[0,3], slices=10);
 // Example: Rotating the pentagon place the zero index at different locations, giving a twist
@@ -315,11 +322,15 @@ include <vnf.scad>
 //              stroke(zrot(30, p=yscale(0.5, p=circle(d=120))),width=10,closed=true);
 //          }
 //      }
-
-
-module skin(profiles, slices, refine=1, method="direct", sampling, caps, closed=false, z, convexity=10)
+module skin(profiles, slices, refine=1, method="direct", sampling, caps, closed=false, z, convexity=10,
+            anchor="origin",cp,spin=0, orient=UP, extent=false)
 {
-      vnf_polyhedron(skin(profiles, slices, refine, method, sampling, caps, closed, z), convexity=convexity);
+    vnf = skin(profiles, slices, refine, method, sampling, caps, closed, z);
+    attachable(anchor=anchor, spin=spin, orient=orient, vnf=vnf, extent=extent, cp=is_def(cp) ? cp : vnf_centroid(vnf))
+    {      
+        vnf_polyhedron(vnf,convexity=convexity);
+        children();
+    }
 }        
 
 
@@ -803,6 +814,12 @@ function associate_vertices(polygons, split, curpoly=0) =
 //   transformations = list of 4x4 matrices to apply
 //   closed = set to true to form a closed (torus) model.  Default: false
 //   caps = true to create endcap faces when closed is false.  Can be a singe boolean to specify endcaps at both ends, or a length 2 boolean array.  Default is true if closed is false.
+//   convexity = convexity setting for use with polyhedron.  (module only) Default: 10
+//   anchor = Translate so anchor point is at the origin.  (module only) Default: "origin"
+//   spin = Rotate this many degrees around Z axis after anchor.  (module only) Default: 0
+//   orient = Vector to rotate top towards after spin  (module only)
+//   extent = use extent method for computing anchors. (module only)  Default: false
+//   cp = set centerpoint for anchor computation.  (module only) Default: object centroid
 // Example: This is the "sweep-drop" example from list-comprehension-demos.
 //   function drop(t) = 100 * 0.5 * (1 - cos(180 * t)) * sin(180 * t) + 1;
 //   function path(t) = [0, 0, 80 + 80 * cos(180 * t)];
@@ -839,9 +856,16 @@ function sweep(shape, transformations, closed=false, caps) =
   assert(!closed || !caps, "Cannot make closed shape with caps")
   _skin_core([for(i=[0:len(transformations)-(closed?0:1)]) apply(transformations[i%len(transformations)],path3d(shape))],caps=fullcaps);
 
-module sweep(shape, transformations, closed=false, caps, convexity=10) {
-  vnf_polyhedron(sweep(shape, transformations, closed, caps), convexity=convexity);
-}  
+module sweep(shape, transformations, closed=false, caps, convexity=10,
+             anchor="origin",cp,spin=0, orient=UP, extent=false)
+{
+    vnf = sweep(shape, transformations, closed, caps);
+    attachable(anchor=anchor, spin=spin, orient=orient, vnf=vnf, extent=extent, cp=is_def(cp) ? cp : vnf_centroid(vnf))
+    {      
+        vnf_polyhedron(vnf,convexity=convexity);
+        children();
+    }
+}        
 
 
 // Function&Module: path_sweep()
@@ -906,8 +930,13 @@ module sweep(shape, transformations, closed=false, caps, convexity=10) {
 //   tangent = a list of tangent vectors in case you need more accuracy (particularly at the end points of your curve)
 //   relaxed = set to true with the "manual" method to relax the orthogonality requirement of cross sections to the path tangent.  Default: false
 //   caps = Can be a boolean or vector of two booleans.  Set to false to disable caps at the two ends.  Default: true
+//   transforms = set to true to return transforms instead of a VNF.  These transforms can be manipulated and passed to sweep().  Default: false.
 //   convexity = convexity parameter for polyhedron().  Only accepted by the module version.  Default: 10
-//   transforms = set to true to return transforms instead of a VNF.  These transforms can be manipulated and passed to sweep().  Default: false. 
+//   anchor = Translate so anchor point is at the origin.  (module only) Default: "origin"
+//   spin = Rotate this many degrees around Z axis after anchor.  (module only) Default: 0
+//   orient = Vector to rotate top towards after spin  (module only)
+//   extent = use extent method for computing anchors. (module only)  Default: false
+//   cp = set centerpoint for anchor computation.  (module only) Default: object centroid
 //
 // Example(2D): We'll use this shape in several examples
 //   ushape = [[-10, 0],[-10, 10],[ -7, 10],[ -7, 2],[  7, 2],[  7, 7],[ 10, 7],[ 10, 0]];
@@ -1121,13 +1150,19 @@ module sweep(shape, transformations, closed=false, caps, convexity=10) {
 //   outside = [for(i=[0:len(trans)-1]) trans[i]*scale(lerp(1,1.5,i/(len(trans)-1)))];
 //   inside = [for(i=[len(trans)-1:-1:0]) trans[i]*scale(lerp(1.1,1.4,i/(len(trans)-1)))];
 //   sweep(shape, concat(outside,inside),closed=true);
-
 module path_sweep(shape, path, method="incremental", normal, closed=false, twist=0, twist_by_length=true,
-                    symmetry=1, last_normal, tangent, relaxed=false, caps, convexity=10)
+                    symmetry=1, last_normal, tangent, relaxed=false, caps, convexity=10,
+                    anchor="origin",cp,spin=0, orient=UP, extent=false)
 {
-  vnf_polyhedron(path_sweep(shape, path, method, normal, closed, twist, twist_by_length,
-                    symmetry, last_normal, tangent, relaxed, caps), convexity=convexity);
-}
+    vnf = path_sweep(shape, path, method, normal, closed, twist, twist_by_length,
+                    symmetry, last_normal, tangent, relaxed, caps);
+    attachable(anchor=anchor, spin=spin, orient=orient, vnf=vnf, extent=extent, cp=is_def(cp) ? cp : vnf_centroid(vnf))
+    {      
+        vnf_polyhedron(vnf,convexity=convexity);
+        children();
+    }
+}        
+
 
 function path_sweep(shape, path, method="incremental", normal, closed=false, twist=0, twist_by_length=true,
                     symmetry=1, last_normal, tangent, relaxed=false, caps, transforms=false) = 
@@ -1202,7 +1237,7 @@ function path_sweep(shape, path, method="incremental", normal, closed=false, twi
             let (pathnormal = path_normals(path, tangents, closed))
             assert(all_defined(pathnormal),"Natural normal vanishes on your curve, select a different method")
             let( testnormals = [for(i=[0:len(pathnormal)-1-(closed?1:2)]) pathnormal[i]*select(pathnormal,i+2)],
-                 dummy = min(testnormals) < .5 ? echo_warning("abrupt change in normal direction.  Consider a different method") :0
+                 dummy = min(testnormals) < .5 ? echo("WARNING: ***** Abrupt change in normal direction.  Consider a different method *****") :0
                )
             [for(i=[0:L-(closed?0:1)]) let(
                      rotation = affine_frame_map(x=pathnormal[i%L], z=tangents[i%L])
@@ -1216,7 +1251,7 @@ function path_sweep(shape, path, method="incremental", normal, closed=false, twi
              end = reindex_polygon(start, apply(transform_list[L],path3d(shape)))
              )
              all([for(i=idx(start)) approx(start[i],end[i])]),
-      dummy = ends_match ? 0 :echo_warning("The points do not match when closing the model")
+      dummy = ends_match ? 0 : echo("WARNING: ***** The points do not match when closing the model *****")
     )
     transforms ? transform_list : sweep(shape, transform_list, closed=false, caps=fullcaps);
 
diff --git a/std.scad b/std.scad
index 8375fcc..1c2e45e 100644
--- a/std.scad
+++ b/std.scad
@@ -14,7 +14,6 @@ include <version.scad>
 include <constants.scad>
 include <edges.scad>
 include <common.scad>
-include <errors.scad>
 include <arrays.scad>
 include <strings.scad>
 include <vnf.scad>
diff --git a/tests/test_arrays.scad b/tests/test_arrays.scad
index f621cd0..8df23a8 100644
--- a/tests/test_arrays.scad
+++ b/tests/test_arrays.scad
@@ -3,14 +3,6 @@ include <../std.scad>
 
 // Section: List Query Operations
 
-module test_is_simple_list() {
-		assert(is_simple_list([1,2,3,4]));
-		assert(is_simple_list([]));
-		assert(!is_simple_list([1,2,[3,4]]));
-} 
-test_is_simple_list();
-
-
 module test_select() {
     l = [3,4,5,6,7,8,9];
     assert(select(l, 5, 6) == [8,9]);
@@ -358,12 +350,34 @@ test_add_scalar();
 module test_subindex() {
     v = [[1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]];
     assert(subindex(v,2) == [3, 7, 11, 15]);
+    assert(subindex(v,[2]) == [[3], [7], [11], [15]]);
     assert(subindex(v,[2,1]) == [[3, 2], [7, 6], [11, 10], [15, 14]]);
     assert(subindex(v,[1:3]) == [[2, 3, 4], [6, 7, 8], [10, 11, 12], [14, 15, 16]]);
 }
 test_subindex();
 
 
+// Need decision about behavior for out of bounds ranges, empty ranges
+module test_submatrix(){
+  M = [[1,2,3,4,5],
+       [6,7,8,9,10],
+       [11,12,13,14,15],
+       [16,17,18,19,20],
+       [21,22,23,24,25]];
+  assert_equal(submatrix(M,[1:2], [3:4]), [[9,10],[14,15]]);
+  assert_equal(submatrix(M,[1], [3,4]), [[9,10]]);
+  assert_equal(submatrix(M,1, [3,4]), [[9,10]]);
+  assert_equal(submatrix(M, [3,4],1), [[17],[22]]);
+  assert_equal(submatrix(M, [1,3],[2,4]), [[8,10],[18,20]]);
+  assert_equal(submatrix(M, 1,3), [[9]]);
+  A = [[true,    17, "test"],
+     [[4,2],   91, false],
+     [6,    [3,4], undef]];
+  assert_equal(submatrix(A,[0,2],[1,2]),[[17, "test"], [[3, 4], undef]]);
+}
+test_submatrix();
+
+
 module test_force_list() {
     assert_equal(force_list([3,4,5]), [3,4,5]);
     assert_equal(force_list(5), [5]);
@@ -466,6 +480,7 @@ test_array_dim();
 module test_transpose() {
     assert(transpose([[1,2,3],[4,5,6],[7,8,9]]) == [[1,4,7],[2,5,8],[3,6,9]]);
     assert(transpose([[1,2,3],[4,5,6]]) == [[1,4],[2,5],[3,6]]);
+    assert(transpose([[1,2,3],[4,5,6]],reverse=true) == [[6,3], [5,2], [4,1]]);
     assert(transpose([3,4,5]) == [3,4,5]);
 }
 test_transpose();
diff --git a/tests/test_coords.scad b/tests/test_coords.scad
index 4b89c12..9fccd65 100644
--- a/tests/test_coords.scad
+++ b/tests/test_coords.scad
@@ -29,6 +29,7 @@ test_point3d();
 
 module test_path3d() {
     assert(path3d([[1,2], [3,4], [5,6], [7,8]])==[[1,2,0],[3,4,0],[5,6,0],[7,8,0]]);
+    assert(path3d([[1,2], [3,4], [5,6], [7,8]],9)==[[1,2,9],[3,4,9],[5,6,9],[7,8,9]]);
     assert(path3d([[1,2,3], [2,3,4], [3,4,5], [4,5,6]])==[[1,2,3],[2,3,4],[3,4,5],[4,5,6]]);
     assert(path3d([[1,2,3,4], [2,3,4,5], [3,4,5,6], [4,5,6,7]])==[[1,2,3],[2,3,4],[3,4,5],[4,5,6]]);
 }
@@ -41,6 +42,9 @@ module test_point4d() {
     assert(point4d([1,2,3])==[1,2,3,0]);
     assert(point4d([2,3])==[2,3,0,0]);
     assert(point4d([1])==[1,0,0,0]);
+    assert(point4d([1,2,3],9)==[1,2,3,9]);
+    assert(point4d([2,3],9)==[2,3,9,9]);
+    assert(point4d([1],9)==[1,9,9,9]);
 }
 test_point4d();
 
diff --git a/tests/test_errors.scad b/tests/test_errors.scad
deleted file mode 100644
index fffe9eb..0000000
--- a/tests/test_errors.scad
+++ /dev/null
@@ -1,11 +0,0 @@
-include <../std.scad>
-
-
-// Can't test echo output as yet.  Include these for coverage calculations.
-module test_echo_error() {}
-module test_echo_warning() {}
-module test_deprecate() {}
-module test_deprecate_argument() {}
-
-
-// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap
diff --git a/tests/test_geometry.scad b/tests/test_geometry.scad
index 8f780f4..0950e1c 100644
--- a/tests/test_geometry.scad
+++ b/tests/test_geometry.scad
@@ -1,6 +1,146 @@
 include <../std.scad>
 
 
+//the commented lines are for tests to be written
+//the tests are ordered as they appear in geometry.scad
+
+test_point_on_segment2d();
+test_point_left_of_line2d();
+test_collinear();
+test_distance_from_line();
+test_line_normal();
+test_line_intersection();
+//test_line_ray_intersection();
+test_line_segment_intersection();
+//test_ray_intersection();
+//test_ray_segment_intersection();
+test_segment_intersection();
+test_line_closest_point();
+//test_ray_closest_point();
+test_segment_closest_point();
+test_line_from_points();
+test_tri_calc();
+//test_hyp_opp_to_adj();
+//test_hyp_ang_to_adj();
+//test_opp_ang_to_adj();
+//test_hyp_adj_to_opp();
+//test_hyp_ang_to_opp();
+//test_adj_ang_to_opp();
+//test_adj_opp_to_hyp();
+//test_adj_ang_to_hyp();
+//test_opp_ang_to_hyp();
+//test_hyp_adj_to_ang();
+//test_hyp_opp_to_ang();
+//test_adj_opp_to_ang();
+test_triangle_area();
+test_plane3pt();
+test_plane3pt_indexed();
+//test_plane_from_normal();
+test_plane_from_points();
+//test_plane_from_polygon();
+test_plane_normal();
+//test_plane_offset();
+//test_plane_transform();
+test_projection_on_plane();
+//test_plane_point_nearest_origin();
+test_distance_from_plane();
+
+test_find_circle_2tangents();
+test_find_circle_3points();
+test_circle_point_tangents();
+test_tri_functions();
+//test_closest_point_on_plane();
+//test__general_plane_line_intersection();
+//test_plane_line_angle();
+//test_plane_line_intersection();
+//test_polygon_line_intersection();
+//test_plane_intersection();
+test_coplanar();
+test_points_on_plane();
+test_in_front_of_plane();
+//test_find_circle_2tangents();
+//test_find_circle_3points();
+//test_circle_point_tangents();
+//test_circle_circle_tangents();
+test_noncollinear_triple();
+test_pointlist_bounds();
+test_closest_point();
+test_furthest_point();
+test_polygon_area();
+test_is_convex_polygon();
+test_polygon_shift();
+test_polygon_shift_to_closest_point();
+test_reindex_polygon();
+test_align_polygon();
+test_centroid();
+test_point_in_polygon();
+test_polygon_is_clockwise();
+test_clockwise_polygon();
+test_ccw_polygon();
+test_reverse_polygon();
+//test_polygon_normal();
+//test_split_polygons_at_each_x();
+//test_split_polygons_at_each_y();
+//test_split_polygons_at_each_z();
+
+//tests to migrate to other files
+test_is_path();
+test_is_closed_path();
+test_close_path();
+test_cleanup_path();
+test_simplify_path();
+test_simplify_path_indexed();
+test_is_region();
+
+// to be used when there are two alternative symmetrical outcomes 
+// from a function like a plane output.
+function standardize(v) = 
+    v==[]? [] :
+    sign([for(vi=v) if( ! approx(vi,0)) vi,0 ][0])*v;
+
+module assert_std(vc,ve) { assert(standardize(vc)==standardize(ve)); }
+
+module test_points_on_plane() {
+    pts     = [for(i=[0:40]) rands(-1,1,3) ];
+    dir     = rands(-10,10,3);
+    normal0 = unit([1,2,3]);
+    ang     = rands(0,360,1)[0];
+    normal  = rot(a=ang,p=normal0);
+    plane   = [each normal, normal*dir];
+    prj_pts = projection_on_plane(plane,pts);
+    assert(points_on_plane(prj_pts,plane));
+    assert(!points_on_plane(concat(pts,[normal-dir]),plane));
+}
+*test_points_on_plane();
+
+module test_projection_on_plane(){
+    ang     = rands(0,360,1)[0];
+    dir     = rands(-10,10,3);
+    normal0 = unit([1,2,3]);
+    normal  = rot(a=ang,p=normal0);
+    plane0  = [each normal0, 0];
+    plane   = [each normal,  0];
+    planem  = [each normal, normal*dir];
+    pts     = [for(i=[1:10]) rands(-1,1,3)];
+    assert_approx( projection_on_plane(plane,pts),
+                   projection_on_plane(plane,projection_on_plane(plane,pts)));
+    assert_approx( projection_on_plane(plane,pts),
+                   rot(a=ang,p=projection_on_plane(plane0,rot(a=-ang,p=pts))));    
+    assert_approx( move((-normal*dir)*normal,p=projection_on_plane(planem,pts)),
+                   projection_on_plane(plane,pts));
+    assert_approx( move((normal*dir)*normal,p=projection_on_plane(plane,pts)),
+                   projection_on_plane(planem,pts));
+}
+*test_projection_on_plane();
+
+module test_line_from_points() {
+    assert_approx(line_from_points([[1,0],[0,0],[-1,0]]),[[-1,0],[1,0]]);
+    assert_approx(line_from_points([[1,1],[0,1],[-1,1]]),[[-1,1],[1,1]]);
+    assert(line_from_points([[1,1],[0,1],[-1,0]])==undef);
+    assert(line_from_points([[1,1],[0,1],[-1,0]],fast=true)== [[-1,0],[1,1]]);
+}
+*test_line_from_points();
+
 module test_point_on_segment2d() {
     assert(point_on_segment2d([-15,0], [[-10,0], [10,0]]) == false);
     assert(point_on_segment2d([-10,0], [[-10,0], [10,0]]) == true);
@@ -29,42 +169,28 @@ module test_point_on_segment2d() {
     assert(point_on_segment2d([ 10, 10], [[-10,-10], [10,10]]) == true);
     assert(point_on_segment2d([ 15, 15], [[-10,-10], [10,10]]) == false);
 }
-test_point_on_segment2d();
+*test_point_on_segment2d();
 
 
-module test_point_left_of_segment() {
-    assert(point_left_of_segment2d([ -3,  0], [[-10,-10], [10,10]]) > 0);
-    assert(point_left_of_segment2d([  0,  0], [[-10,-10], [10,10]]) == 0);
-    assert(point_left_of_segment2d([  3,  0], [[-10,-10], [10,10]]) < 0);
+module test_point_left_of_line2d() {
+    assert(point_left_of_line2d([ -3,  0], [[-10,-10], [10,10]]) > 0);
+    assert(point_left_of_line2d([  0,  0], [[-10,-10], [10,10]]) == 0);
+    assert(point_left_of_line2d([  3,  0], [[-10,-10], [10,10]]) < 0);
 }
-test_point_left_of_segment();
-
+*test_point_left_of_line2d();
 
 module test_collinear() {
     assert(collinear([-10,-10], [-15, -16], [10,10]) == false);
+    assert(collinear([[-10,-10], [-15, -16], [10,10]]) == false);
     assert(collinear([-10,-10], [-15, -15], [10,10]) == true);
+    assert(collinear([[-10,-10], [-15, -15], [10,10]]) == true);
     assert(collinear([-10,-10], [ -3,   0], [10,10]) == false);
     assert(collinear([-10,-10], [  0,   0], [10,10]) == true);
     assert(collinear([-10,-10], [  3,   0], [10,10]) == false);
     assert(collinear([-10,-10], [ 15,  15], [10,10]) == true);
     assert(collinear([-10,-10], [ 15,  16], [10,10]) == false);
 }
-test_collinear();
-
-
-module test_collinear_indexed() {
-    pts = [
-        [-20,-20], [-10,-20], [0,-10], [10,0], [20,10], [20,20], [15,30]
-    ];
-    assert(collinear_indexed(pts, 0,1,2) == false);
-    assert(collinear_indexed(pts, 1,2,3) == true);
-    assert(collinear_indexed(pts, 2,3,4) == true);
-    assert(collinear_indexed(pts, 3,4,5) == false);
-    assert(collinear_indexed(pts, 4,5,6) == false);
-    assert(collinear_indexed(pts, 4,3,2) == true);
-    assert(collinear_indexed(pts, 0,5,6) == false);
-}
-test_collinear_indexed();
+*test_collinear();
 
 
 module test_distance_from_line() {
@@ -73,7 +199,7 @@ module test_distance_from_line() {
     assert(abs(distance_from_line([[-10,-10,-10], [10,10,10]], [1,-1,0]) - sqrt(2)) < EPSILON);
     assert(abs(distance_from_line([[-10,-10,-10], [10,10,10]], [8,-8,0]) - 8*sqrt(2)) < EPSILON);
 }
-test_distance_from_line();
+*test_distance_from_line();
 
 
 module test_line_normal() {
@@ -97,7 +223,7 @@ module test_line_normal() {
         assert(approx(n2, n1));
     }
 }
-test_line_normal();
+*test_line_normal();
 
 
 module test_line_intersection() {
@@ -110,7 +236,7 @@ module test_line_intersection() {
     assert(line_intersection([[-10,-10], [ 10, 10]], [[ 10,-10], [-10, 10]]) == [0,0]);
     assert(line_intersection([[ -8,  0], [ 12,  4]], [[ 12,  0], [ -8,  4]]) == [2,2]);
 }
-test_line_intersection();
+*test_line_intersection();
 
 
 module test_segment_intersection() {
@@ -126,7 +252,7 @@ module test_segment_intersection() {
     assert(segment_intersection([[-10,-10], [ 10, 10]], [[ 10,-10], [-10, 10]]) == [0,0]);
     assert(segment_intersection([[ -8,  0], [ 12,  4]], [[ 12,  0], [ -8,  4]]) == [2,2]);
 }
-test_segment_intersection();
+*test_segment_intersection();
 
 
 module test_line_segment_intersection() {
@@ -141,7 +267,7 @@ module test_line_segment_intersection() {
     assert(line_segment_intersection([[-10,-10], [ 10, 10]], [[ 10,-10], [  1, -1]]) == undef);
     assert(line_segment_intersection([[-10,-10], [ 10, 10]], [[ 10,-10], [ -1,  1]]) == [0,0]);
 }
-test_line_segment_intersection();
+*test_line_segment_intersection();
 
 
 module test_line_closest_point() {
@@ -151,7 +277,7 @@ module test_line_closest_point() {
     assert(approx(line_closest_point([[-10,-20], [10,20]], [1,2]+[2,-1]), [1,2]));
     assert(approx(line_closest_point([[-10,-20], [10,20]], [13,31]), [15,30]));
 }
-test_line_closest_point();
+*test_line_closest_point();
 
 
 module test_segment_closest_point() {
@@ -162,10 +288,10 @@ module test_segment_closest_point() {
     assert(approx(segment_closest_point([[-10,-20], [10,20]], [13,31]), [10,20]));
     assert(approx(segment_closest_point([[-10,-20], [10,20]], [15,25]), [10,20]));
 }
-test_segment_closest_point();
-
+*test_segment_closest_point();
 
 module test_find_circle_2tangents() {
+//** missing tests with arg tangent=true
     assert(approx(find_circle_2tangents([10,10],[0,0],[10,-10],r=10/sqrt(2))[0],[10,0]));
     assert(approx(find_circle_2tangents([-10,10],[0,0],[-10,-10],r=10/sqrt(2))[0],[-10,0]));
     assert(approx(find_circle_2tangents([-10,10],[0,0],[10,10],r=10/sqrt(2))[0],[0,10]));
@@ -174,9 +300,9 @@ module test_find_circle_2tangents() {
     assert(approx(find_circle_2tangents([10,0],[0,0],[0,-10],r=10)[0],[10,-10]));
     assert(approx(find_circle_2tangents([0,-10],[0,0],[-10,0],r=10)[0],[-10,-10]));
     assert(approx(find_circle_2tangents([-10,0],[0,0],[0,10],r=10)[0],[-10,10]));
-    assert(approx(find_circle_2tangents(polar_to_xy(10,60),[0,0],[10,0],r=10)[0],polar_to_xy(20,30)));
+    assert_approx(find_circle_2tangents(polar_to_xy(10,60),[0,0],[10,0],r=10)[0],polar_to_xy(20,30));
 }
-test_find_circle_2tangents();
+*test_find_circle_2tangents();
 
 
 module test_find_circle_3points() {
@@ -291,7 +417,7 @@ module test_find_circle_3points() {
         }
     }
 }
-test_find_circle_3points();
+*test_find_circle_3points();
 
 
 module test_circle_point_tangents() {
@@ -304,7 +430,7 @@ module test_circle_point_tangents() {
         assert(approx(flatten(got), flatten(expected)));
     }
 }
-test_circle_point_tangents();
+*test_circle_point_tangents();
 
 
 module test_tri_calc() {
@@ -327,23 +453,9 @@ module test_tri_calc() {
         assert(approx(tri_calc(hyp=hyp, ang2=ang2), expected));
     }
 }
-test_tri_calc();
+*test_tri_calc();
 
 
-// Dummy modules to show up in coverage check script.
-module test_hyp_opp_to_adj();
-module test_hyp_ang_to_adj();
-module test_opp_ang_to_adj();
-module test_hyp_adj_to_opp();
-module test_hyp_ang_to_opp();
-module test_adj_ang_to_opp();
-module test_adj_opp_to_hyp();
-module test_adj_ang_to_hyp();
-module test_opp_ang_to_hyp();
-module test_hyp_adj_to_ang();
-module test_hyp_opp_to_ang();
-module test_adj_opp_to_ang();
-
 module test_tri_functions() {
     sides = rands(1,100,100,seed_value=8181);
     for (p = pair_wrap(sides)) {
@@ -365,7 +477,7 @@ module test_tri_functions() {
         assert_approx(adj_opp_to_ang(adj,opp), ang);
     }
 }
-test_tri_functions();
+*test_tri_functions();
 
 
 module test_triangle_area() {
@@ -373,55 +485,53 @@ module test_triangle_area() {
     assert(abs(triangle_area([0,0], [0,10], [0,15])) < EPSILON);
     assert(abs(triangle_area([0,0], [10,0], [0,10]) - 50) < EPSILON);
 }
-test_triangle_area();
+*test_triangle_area();
 
 
 module test_plane3pt() {
-    assert(plane3pt([0,0,20], [0,10,10], [0,0,0]) == [1,0,0,0]);
-    assert(plane3pt([2,0,20], [2,10,10], [2,0,0]) == [1,0,0,2]);
-    assert(plane3pt([0,0,0], [10,0,10], [0,0,20]) == [0,1,0,0]);
-    assert(plane3pt([0,2,0], [10,2,10], [0,2,20]) == [0,1,0,2]);
-    assert(plane3pt([0,0,0], [10,10,0], [20,0,0]) == [0,0,1,0]);
-    assert(plane3pt([0,0,2], [10,10,2], [20,0,2]) == [0,0,1,2]);
+    assert_std(plane3pt([0,0,20], [0,10,10], [0,0,0]), [1,0,0,0]);
+    assert_std(plane3pt([2,0,20], [2,10,10], [2,0,0]), [1,0,0,2]);
+    assert_std(plane3pt([0,0,0], [10,0,10], [0,0,20]), [0,1,0,0]);
+    assert_std(plane3pt([0,2,0], [10,2,10], [0,2,20]), [0,1,0,2]);
+    assert_std(plane3pt([0,0,0], [10,10,0], [20,0,0]), [0,0,1,0]);
+    assert_std(plane3pt([0,0,2], [10,10,2], [20,0,2]), [0,0,1,2]);
 }
-test_plane3pt();
-
+*test_plane3pt();
 
 module test_plane3pt_indexed() {
     pts = [ [0,0,0], [10,0,0], [0,10,0], [0,0,10] ];
     s13 = sqrt(1/3);
-    assert(plane3pt_indexed(pts, 0,3,2) == [1,0,0,0]);
-    assert(plane3pt_indexed(pts, 0,2,3) == [-1,0,0,0]);
-    assert(plane3pt_indexed(pts, 0,1,3) == [0,1,0,0]);
-    assert(plane3pt_indexed(pts, 0,3,1) == [0,-1,0,0]);
-    assert(plane3pt_indexed(pts, 0,2,1) == [0,0,1,0]);
-    assert(plane3pt_indexed(pts, 0,1,2) == [0,0,-1,0]);
-    assert(plane3pt_indexed(pts, 3,2,1) == [s13,s13,s13,10*s13]);
-    assert(plane3pt_indexed(pts, 1,2,3) == [-s13,-s13,-s13,-10*s13]);
+    assert_std(plane3pt_indexed(pts, 0,3,2), [1,0,0,0]);
+    assert_std(plane3pt_indexed(pts, 0,2,3), [-1,0,0,0]);
+    assert_std(plane3pt_indexed(pts, 0,1,3), [0,1,0,0]);
+    assert_std(plane3pt_indexed(pts, 0,3,1), [0,-1,0,0]);
+    assert_std(plane3pt_indexed(pts, 0,2,1), [0,0,1,0]);
+    assert_approx(plane3pt_indexed(pts, 0,1,2), [0,0,-1,0]);
+    assert_approx(plane3pt_indexed(pts, 3,2,1), [s13,s13,s13,10*s13]);
+    assert_approx(plane3pt_indexed(pts, 1,2,3), [-s13,-s13,-s13,-10*s13]);
 }
-test_plane3pt_indexed();
-
+*test_plane3pt_indexed();
 
 module test_plane_from_points() {
-    assert(plane_from_points([[0,0,20], [0,10,10], [0,0,0], [0,5,3]]) == [1,0,0,0]);
-    assert(plane_from_points([[2,0,20], [2,10,10], [2,0,0], [2,3,4]]) == [1,0,0,2]);
-    assert(plane_from_points([[0,0,0], [10,0,10], [0,0,20], [5,0,7]]) == [0,1,0,0]);
-    assert(plane_from_points([[0,2,0], [10,2,10], [0,2,20], [4,2,3]]) == [0,1,0,2]);
-    assert(plane_from_points([[0,0,0], [10,10,0], [20,0,0], [8,3,0]]) == [0,0,1,0]);
-    assert(plane_from_points([[0,0,2], [10,10,2], [20,0,2], [3,4,2]]) == [0,0,1,2]);
+    assert_std(plane_from_points([[0,0,20], [0,10,10], [0,0,0], [0,5,3]]), [1,0,0,0]);
+    assert_std(plane_from_points([[2,0,20], [2,10,10], [2,0,0], [2,3,4]]), [1,0,0,2]);
+    assert_std(plane_from_points([[0,0,0], [10,0,10], [0,0,20], [5,0,7]]), [0,1,0,0]);
+    assert_std(plane_from_points([[0,2,0], [10,2,10], [0,2,20], [4,2,3]]), [0,1,0,2]);
+    assert_std(plane_from_points([[0,0,0], [10,10,0], [20,0,0], [8,3,0]]), [0,0,1,0]);
+    assert_std(plane_from_points([[0,0,2], [10,10,2], [20,0,2], [3,4,2]]), [0,0,1,2]);
 }
-test_plane_from_points();
+*test_plane_from_points();
 
 
 module test_plane_normal() {
-    assert(plane_normal(plane3pt([0,0,20], [0,10,10], [0,0,0])) == [1,0,0]);
-    assert(plane_normal(plane3pt([2,0,20], [2,10,10], [2,0,0])) == [1,0,0]);
-    assert(plane_normal(plane3pt([0,0,0], [10,0,10], [0,0,20])) == [0,1,0]);
-    assert(plane_normal(plane3pt([0,2,0], [10,2,10], [0,2,20])) == [0,1,0]);
-    assert(plane_normal(plane3pt([0,0,0], [10,10,0], [20,0,0])) == [0,0,1]);
-    assert(plane_normal(plane3pt([0,0,2], [10,10,2], [20,0,2])) == [0,0,1]);
+    assert_std(plane_normal(plane3pt([0,0,20], [0,10,10], [0,0,0])), [1,0,0]);
+    assert_std(plane_normal(plane3pt([2,0,20], [2,10,10], [2,0,0])), [1,0,0]);
+    assert_std(plane_normal(plane3pt([0,0,0], [10,0,10], [0,0,20])), [0,1,0]);
+    assert_std(plane_normal(plane3pt([0,2,0], [10,2,10], [0,2,20])), [0,1,0]);
+    assert_std(plane_normal(plane3pt([0,0,0], [10,10,0], [20,0,0])), [0,0,1]);
+    assert_std(plane_normal(plane3pt([0,0,2], [10,10,2], [20,0,2])), [0,0,1]);
 }
-test_plane_normal();
+*test_plane_normal();
 
 
 module test_distance_from_plane() {
@@ -429,20 +539,16 @@ module test_distance_from_plane() {
     assert(distance_from_plane(plane1, [0,0,5]) == 5);
     assert(distance_from_plane(plane1, [5,5,8]) == 8);
 }
-test_distance_from_plane();
+*test_distance_from_plane();
 
 
 module test_coplanar() {
-    plane = plane3pt([0,0,0], [0,10,10], [10,0,10]);
-    assert(coplanar(plane, [5,5,10]) == true);
-    assert(coplanar(plane, [10/3,10/3,20/3]) == true);
-    assert(coplanar(plane, [0,0,0]) == true);
-    assert(coplanar(plane, [1,1,0]) == false);
-    assert(coplanar(plane, [-1,1,0]) == true);
-    assert(coplanar(plane, [1,-1,0]) == true);
-    assert(coplanar(plane, [5,5,5]) == false);
+    assert(coplanar([ [5,5,1],[0,0,1],[-1,-1,1] ]) == false);
+    assert(coplanar([ [5,5,1],[0,0,0],[-1,-1,1] ]) == true);
+    assert(coplanar([ [0,0,0],[1,0,1],[1,1,1], [0,1,2] ]) == false);
+    assert(coplanar([ [0,0,0],[1,0,1],[1,1,2], [0,1,1] ]) == true);
 }
-test_coplanar();
+*test_coplanar();
 
 
 module test_in_front_of_plane() {
@@ -455,7 +561,7 @@ module test_in_front_of_plane() {
     assert(in_front_of_plane(plane, [0,0,5]) == true);
     assert(in_front_of_plane(plane, [0,0,-5]) == false);
 }
-test_in_front_of_plane();
+*test_in_front_of_plane();
 
 
 module test_is_path() {
@@ -470,35 +576,43 @@ module test_is_path() {
     assert(is_path([[1,2,3],[4,5,6]]));
     assert(is_path([[1,2,3],[4,5,6],[7,8,9]]));
 }
-test_is_path();
+*test_is_path();
 
 
 module test_is_closed_path() {
     assert(!is_closed_path([[1,2,3],[4,5,6],[1,8,9]]));
     assert(is_closed_path([[1,2,3],[4,5,6],[1,8,9],[1,2,3]]));
 }
-test_is_closed_path();
+*test_is_closed_path();
 
 
 module test_close_path() {
     assert(close_path([[1,2,3],[4,5,6],[1,8,9]]) == [[1,2,3],[4,5,6],[1,8,9],[1,2,3]]);
     assert(close_path([[1,2,3],[4,5,6],[1,8,9],[1,2,3]]) == [[1,2,3],[4,5,6],[1,8,9],[1,2,3]]);
 }
-test_close_path();
+*test_close_path();
 
 
 module test_cleanup_path() {
     assert(cleanup_path([[1,2,3],[4,5,6],[1,8,9]]) == [[1,2,3],[4,5,6],[1,8,9]]);
     assert(cleanup_path([[1,2,3],[4,5,6],[1,8,9],[1,2,3]]) == [[1,2,3],[4,5,6],[1,8,9]]);
 }
-test_cleanup_path();
+*test_cleanup_path();
 
 
 module test_polygon_area() {
     assert(approx(polygon_area([[1,1],[-1,1],[-1,-1],[1,-1]]), 4));
     assert(approx(polygon_area(circle(r=50,$fn=1000)), -PI*50*50, eps=0.1));
 }
-test_polygon_area();
+*test_polygon_area();
+
+
+module test_is_convex_polygon() {
+    assert(is_convex_polygon([[1,1],[-1,1],[-1,-1],[1,-1]]));
+    assert(is_convex_polygon(circle(r=50,$fn=1000)));
+    assert(!is_convex_polygon([[1,1],[0,0],[-1,1],[-1,-1],[1,-1]]));
+}
+*test_is_convex_polygon();
 
 
 module test_polygon_shift() {
@@ -506,7 +620,7 @@ module test_polygon_shift() {
     assert(polygon_shift(path,1) == [[-1,1],[-1,-1],[1,-1],[1,1]]);
     assert(polygon_shift(path,2) == [[-1,-1],[1,-1],[1,1],[-1,1]]);
 }
-test_polygon_shift();
+*test_polygon_shift();
 
 
 module test_polygon_shift_to_closest_point() {
@@ -516,56 +630,45 @@ module test_polygon_shift_to_closest_point() {
     assert(polygon_shift_to_closest_point(path,[-1.1,-1.1]) == [[-1,-1],[1,-1],[1,1],[-1,1]]);
     assert(polygon_shift_to_closest_point(path,[1.1,-1.1]) == [[1,-1],[1,1],[-1,1],[-1,-1]]);
 }
-test_polygon_shift_to_closest_point();
+*test_polygon_shift_to_closest_point();
 
 
-/*
-module test_first_noncollinear(){
-    pts = [
-        [1,1], [2,2], [3,3], [4,4], [4,5], [5,6]
-    ];
-    assert(first_noncollinear(0,1,pts) == 4);
-    assert(first_noncollinear(1,0,pts) == 4);
-    assert(first_noncollinear(0,2,pts) == 4);
-    assert(first_noncollinear(2,0,pts) == 4);
-    assert(first_noncollinear(1,2,pts) == 4);
-    assert(first_noncollinear(2,1,pts) == 4);
-    assert(first_noncollinear(0,3,pts) == 4);
-    assert(first_noncollinear(3,0,pts) == 4);
-    assert(first_noncollinear(1,3,pts) == 4);
-    assert(first_noncollinear(3,1,pts) == 4);
-    assert(first_noncollinear(2,3,pts) == 4);
-    assert(first_noncollinear(3,2,pts) == 4);
-    assert(first_noncollinear(0,4,pts) == 1);
-    assert(first_noncollinear(4,0,pts) == 1);
-    assert(first_noncollinear(1,4,pts) == 0);
-    assert(first_noncollinear(4,1,pts) == 0);
-    assert(first_noncollinear(2,4,pts) == 0);
-    assert(first_noncollinear(4,2,pts) == 0);
-    assert(first_noncollinear(3,4,pts) == 0);
-    assert(first_noncollinear(4,3,pts) == 0);
-    assert(first_noncollinear(0,5,pts) == 1);
-    assert(first_noncollinear(5,0,pts) == 1);
-    assert(first_noncollinear(1,5,pts) == 0);
-    assert(first_noncollinear(5,1,pts) == 0);
-    assert(first_noncollinear(2,5,pts) == 0);
-    assert(first_noncollinear(5,2,pts) == 0);
-    assert(first_noncollinear(3,5,pts) == 0);
-    assert(first_noncollinear(5,3,pts) == 0);
-    assert(first_noncollinear(4,5,pts) == 0);
-    assert(first_noncollinear(5,4,pts) == 0);
+module test_reindex_polygon() {
+   pent = subdivide_path([for(i=[0:4])[sin(72*i),cos(72*i)]],5);
+   circ = circle($fn=5,r=2.2);
+   assert_approx(reindex_polygon(circ,pent), [[0.951056516295,0.309016994375],[0.587785252292,-0.809016994375],[-0.587785252292,-0.809016994375],[-0.951056516295,0.309016994375],[0,1]]);
+   poly = [[-1,1],[-1,-1],[1,-1],[1,1],[0,0]];
+   ref  = [for(i=[0:4])[sin(72*i),cos(72*i)]];
+   assert_approx(reindex_polygon(ref,poly),[[0,0],[1,1],[1,-1],[-1,-1],[-1,1]]);
 }
-test_first_noncollinear();
-*/
+*test_reindex_polygon();
 
 
-module test_find_noncollinear_points() {
-    assert(find_noncollinear_points([[1,1],[2,2],[3,3],[4,4],[4,5],[5,6]]) == [0,5,3]);
-    assert(find_noncollinear_points([[1,1],[2,2],[8,3],[4,4],[4,5],[5,6]]) == [0,2,5]);
+module test_align_polygon() {
+   pentagon = subdivide_path(pentagon(side=2),10);
+   hexagon  = subdivide_path(hexagon(side=2.7),10);
+   aligned =  [[2.7,0],[2.025,-1.16913429511],[1.35,-2.33826859022],
+               [-1.35,-2.33826859022],[-2.025,-1.16913429511],[-2.7,0],
+               [-2.025,1.16913429511],[-1.35,2.33826859022],[1.35,2.33826859022],
+               [2.025,1.16913429511]];
+   assert_approx(align_polygon(pentagon,hexagon,[0:10:359]), aligned);
+   aligned2 = [[1.37638192047,0],[1.37638192047,-1],[0.425325404176,-1.30901699437],
+               [-0.525731112119,-1.61803398875],[-1.11351636441,-0.809016994375],
+               [-1.7013016167,0],[-1.11351636441,0.809016994375],
+               [-0.525731112119,1.61803398875],[0.425325404176,1.30901699437],
+               [1.37638192047,1]];
+   assert_approx(align_polygon(hexagon,pentagon,[0:10:359]), aligned2);
+}
+*test_align_polygon();
+
+
+module test_noncollinear_triple() {
+    assert(noncollinear_triple([[1,1],[2,2],[3,3],[4,4],[4,5],[5,6]]) == [0,5,3]);
+    assert(noncollinear_triple([[1,1],[2,2],[8,3],[4,4],[4,5],[5,6]]) == [0,2,5]);
     u = unit([5,3]);
-    assert_equal(find_noncollinear_points([for(i = [2,3,4,5,7,12,15]) i * u], error=false),[]);
+    assert_equal(noncollinear_triple([for(i = [2,3,4,5,7,12,15]) i * u], error=false),[]);
 }
-test_find_noncollinear_points();
+*test_noncollinear_triple();
 
 
 module test_centroid() {
@@ -573,15 +676,18 @@ module test_centroid() {
     assert_approx(centroid(circle(d=100)), [0,0]);
     assert_approx(centroid(rect([40,60],rounding=10,anchor=LEFT)), [20,0]);
     assert_approx(centroid(rect([40,60],rounding=10,anchor=FWD)), [0,30]);
+    poly = [for(a=[0:90:360]) 
+              move([1,2.5,3.1], rot(p=[cos(a),sin(a),0],from=[0,0,1],to=[1,1,1])) ];
+    assert_approx(centroid(poly), [1,2.5,3.1]);
 }
-test_centroid();
+*test_centroid();
 
 
 module test_simplify_path() {
     path = [[-20,-20], [-10,-20], [0,-10], [10,0], [20,10], [20,20], [15,30]];
     assert(simplify_path(path) == [[-20,-20], [-10,-20], [20,10], [20,20], [15,30]]);
 }
-test_simplify_path();
+*test_simplify_path();
 
 
 module test_simplify_path_indexed() {
@@ -589,23 +695,29 @@ module test_simplify_path_indexed() {
     path = [4,6,1,0,3,2,5];
     assert(simplify_path_indexed(pts, path) == [4,6,3,2,5]);
 }
-test_simplify_path_indexed();
+*test_simplify_path_indexed();
 
 
 module test_point_in_polygon() {
     poly = [for (a=[0:30:359]) 10*[cos(a),sin(a)]];
+    poly2 = [ [-3,-3],[2,-3],[2,1],[-1,1],[-1,-1],[1,-1],[1,2],[-3,2] ];
     assert(point_in_polygon([0,0], poly) == 1);
     assert(point_in_polygon([20,0], poly) == -1);
+    assert(point_in_polygon([20,0], poly,EPSILON,nonzero=false) == -1);
     assert(point_in_polygon([5,5], poly) == 1);
     assert(point_in_polygon([-5,5], poly) == 1);
     assert(point_in_polygon([-5,-5], poly) == 1);
     assert(point_in_polygon([5,-5], poly) == 1);
+    assert(point_in_polygon([5,-5], poly,EPSILON,nonzero=false) == 1);
     assert(point_in_polygon([-10,-10], poly) == -1);
     assert(point_in_polygon([10,0], poly) == 0);
     assert(point_in_polygon([0,10], poly) == 0);
     assert(point_in_polygon([0,-10], poly) == 0);
+    assert(point_in_polygon([0,-10], poly,EPSILON,nonzero=false) == 0);
+    assert(point_in_polygon([0,0], poly2,EPSILON,nonzero=true) == 1);
+    assert(point_in_polygon([0,0], poly2,EPSILON,nonzero=false) == -1);
 }
-test_point_in_polygon();
+*test_point_in_polygon();
 
 
 module test_pointlist_bounds() {
@@ -626,16 +738,16 @@ module test_pointlist_bounds() {
     ];
     assert(pointlist_bounds(pts2d) == [[-63,-42],[84,42]]);
     pts5d = [
-        [-53,27,12,-53,12],
-        [-63,97,36,-63,36],
-        [84,-32,-5,84,-5], 
-        [63,-24,42,63,42], 
-        [23,57,-42,23,-42]
+        [-53, 27, 12,-53, 12],
+        [-63, 97, 36,-63, 36],
+        [ 84,-32, -5, 84, -5], 
+        [ 63,-24, 42, 63, 42], 
+        [ 23, 57,-42, 23,-42]
     ];
     assert(pointlist_bounds(pts5d) == [[-63,-32,-42,-63,-42],[84,97,42,84,42]]);
     assert(pointlist_bounds([[3,4,5,6]]), [[3,4,5,6],[3,4,5,6]]);
 }
-test_pointlist_bounds();
+*test_pointlist_bounds();
 
 
 module test_closest_point() {
@@ -648,7 +760,7 @@ module test_closest_point() {
         assert(mindist == dists[pidx]);
     }
 }
-test_closest_point();
+*test_closest_point();
 
 
 module test_furthest_point() {
@@ -661,7 +773,7 @@ module test_furthest_point() {
         assert(mindist == dists[pidx]);
     }
 }
-test_furthest_point();
+*test_furthest_point();
 
 
 module test_polygon_is_clockwise() {
@@ -670,7 +782,7 @@ module test_polygon_is_clockwise() {
     assert(polygon_is_clockwise(circle(d=100)));
     assert(polygon_is_clockwise(square(100)));
 }
-test_polygon_is_clockwise();
+*test_polygon_is_clockwise();
 
 
 module test_clockwise_polygon() {
@@ -679,7 +791,7 @@ module test_clockwise_polygon() {
     assert(clockwise_polygon(path) == path);
     assert(clockwise_polygon(rpath) == path);
 }
-test_clockwise_polygon();
+*test_clockwise_polygon();
 
 
 module test_ccw_polygon() {
@@ -688,7 +800,7 @@ module test_ccw_polygon() {
     assert(ccw_polygon(path) == rpath);
     assert(ccw_polygon(rpath) == rpath);
 }
-test_ccw_polygon();
+*test_ccw_polygon();
 
 
 module test_reverse_polygon() {
@@ -697,7 +809,7 @@ module test_reverse_polygon() {
     assert(reverse_polygon(path) == rpath);
     assert(reverse_polygon(rpath) == path);
 }
-test_reverse_polygon();
+*test_reverse_polygon();
 
 
 module test_is_region() {
@@ -709,7 +821,7 @@ module test_is_region() {
     assert(!is_region(true));
     assert(!is_region("foo"));
 }
-test_is_region();
+*test_is_region();
 
 
 
diff --git a/tests/test_math.scad b/tests/test_math.scad
index f1b36b9..c3eb601 100644
--- a/tests/test_math.scad
+++ b/tests/test_math.scad
@@ -100,6 +100,106 @@ module test_is_matrix() {
 test_is_matrix();
 
 
+module test_is_zero() {
+    assert(is_zero(0));
+    assert(is_zero([0,0,0]));
+    assert(is_zero([[0,0,0],[0,0]]));
+    assert(is_zero([EPSILON/2,EPSILON/2,EPSILON/2]));
+    assert(!is_zero(1e-3));
+    assert(!is_zero([0,0,1e-3]));
+    assert(!is_zero([EPSILON*10,0,0]));
+    assert(!is_zero([0,EPSILON*10,0]));
+    assert(!is_zero([0,0,EPSILON*10]));
+    assert(!is_zero(true));
+    assert(!is_zero(false));
+    assert(!is_zero(INF));
+    assert(!is_zero(-INF));
+    assert(!is_zero(NAN));
+    assert(!is_zero("foo"));
+    assert(!is_zero([]));
+    assert(!is_zero([0:1:2]));
+}
+test_is_zero();
+
+
+module test_is_positive() {
+    assert(!is_positive(-2));
+    assert(!is_positive(0));
+    assert(is_positive(2));
+    assert(!is_positive([0,0,0]));
+    assert(!is_positive([0,1,2]));
+    assert(is_positive([3,1,2]));
+    assert(!is_positive([3,-1,2]));
+    assert(!is_positive([]));
+    assert(!is_positive(true));
+    assert(!is_positive(false));
+    assert(!is_positive("foo"));
+    assert(!is_positive([0:1:2]));
+}
+test_is_positive();
+
+
+module test_is_negative() {
+    assert(is_negative(-2));
+    assert(!is_negative(0));
+    assert(!is_negative(2));
+    assert(!is_negative([0,0,0]));
+    assert(!is_negative([0,1,2]));
+    assert(!is_negative([3,1,2]));
+    assert(!is_negative([3,-1,2]));
+    assert(is_negative([-3,-1,-2]));
+    assert(!is_negative([-3,1,-2]));
+    assert(is_negative([[-5,-7],[-3,-1,-2]]));
+    assert(!is_negative([[-5,-7],[-3,1,-2]]));
+    assert(!is_negative([]));
+    assert(!is_negative(true));
+    assert(!is_negative(false));
+    assert(!is_negative("foo"));
+    assert(!is_negative([0:1:2]));
+}
+test_is_negative();
+
+
+module test_is_nonpositive() {
+    assert(is_nonpositive(-2));
+    assert(is_nonpositive(0));
+    assert(!is_nonpositive(2));
+    assert(is_nonpositive([0,0,0]));
+    assert(!is_nonpositive([0,1,2]));
+    assert(is_nonpositive([0,-1,-2]));
+    assert(!is_nonpositive([3,1,2]));
+    assert(!is_nonpositive([3,-1,2]));
+    assert(!is_nonpositive([]));
+    assert(!is_nonpositive(true));
+    assert(!is_nonpositive(false));
+    assert(!is_nonpositive("foo"));
+    assert(!is_nonpositive([0:1:2]));
+}
+test_is_nonpositive();
+
+
+module test_is_nonnegative() {
+    assert(!is_nonnegative(-2));
+    assert(is_nonnegative(0));
+    assert(is_nonnegative(2));
+    assert(is_nonnegative([0,0,0]));
+    assert(is_nonnegative([0,1,2]));
+    assert(is_nonnegative([3,1,2]));
+    assert(!is_nonnegative([3,-1,2]));
+    assert(!is_nonnegative([-3,-1,-2]));
+    assert(!is_nonnegative([[-5,-7],[-3,-1,-2]]));
+    assert(!is_nonnegative([[-5,-7],[-3,1,-2]]));
+    assert(!is_nonnegative([[5,7],[3,-1,2]]));
+    assert(is_nonnegative([[5,7],[3,1,2]]));
+    assert(!is_nonnegative([]));
+    assert(!is_nonnegative(true));
+    assert(!is_nonnegative(false));
+    assert(!is_nonnegative("foo"));
+    assert(!is_nonnegative([0:1:2]));
+}
+test_is_nonnegative();
+
+
 module test_approx() {
     assert_equal(approx(PI, 3.141592653589793236), true);
     assert_equal(approx(PI, 3.1415926), false);
@@ -391,11 +491,13 @@ module test_mean() {
 }
 test_mean();
 
+/*
 module test_median() {
     assert_equal(median([2,3,7]), 4.5);
     assert_equal(median([[1,2,3], [3,4,5], [8,9,10]]), [4.5,5.5,6.5]);
 }
 test_median();
+*/
 
 
 module test_convolve() {
@@ -851,22 +953,6 @@ module test_real_roots(){
 }
 test_real_roots();
 
-// Need decision about behavior for out of bounds ranges, empty ranges
-module test_submatrix(){
-  M = [[1,2,3,4,5],
-       [6,7,8,9,10],
-       [11,12,13,14,15],
-       [16,17,18,19,20],
-       [21,22,23,24,25]];
-  assert_equal(submatrix(M,[1:2], [3:4]), [[9,10],[14,15]]);
-  assert_equal(submatrix(M,[1], [3,4]), [[9,10]]);
-  assert_equal(submatrix(M,1, [3,4]), [[9,10]]);
-  assert_equal(submatrix(M, [3,4],1), [[17],[22]]);
-  assert_equal(submatrix(M, [1,3],[2,4]), [[8,10],[18,20]]);
-}
-test_submatrix();
-
-
 
 module test_qr_factor() {
   // Check that R is upper triangular
@@ -911,23 +997,21 @@ test_qr_factor();
 
 module test_poly_mult(){
   assert_equal(poly_mult([3,2,1],[4,5,6,7]),[12,23,32,38,20,7]);
-  assert_equal(poly_mult([3,2,1],[0]),[0]);
-//  assert_equal(poly_mult([3,2,1],[]),[]);
   assert_equal(poly_mult([[1,2],[3,4],[5,6]]), [15,68,100,48]);
+  assert_equal(poly_mult([3,2,1],[0]),[0]);
   assert_equal(poly_mult([[1,2],[0],[5,6]]), [0]);
-//  assert_equal(poly_mult([[1,2],[],[5,6]]), []);
-  assert_equal(poly_mult([[3,4,5],[0,0,0]]),[0]);
-//  assert_equal(poly_mult([[3,4,5],[0,0,0]]),[]);
+  assert_equal(poly_mult([[3,4,5],[0,0,0]]), [0]);
+  assert_equal(poly_mult([[0],[0,0,0]]),[0]);
 }
 test_poly_mult();
 
 
 module test_poly_div(){
   assert_equal(poly_div(poly_mult([4,3,3,2],[2,1,3]), [2,1,3]),[[4,3,3,2],[0]]);
-//  assert_equal(poly_div(poly_mult([4,3,3,2],[2,1,3]), [2,1,3]),[[4,3,3,2],[]]);
   assert_equal(poly_div([1,2,3,4],[1,2,3,4,5]), [[], [1,2,3,4]]);
   assert_equal(poly_div(poly_add(poly_mult([1,2,3,4],[2,0,2]), [1,1,2]), [1,2,3,4]), [[2,0,2],[1,1,2]]);
   assert_equal(poly_div([1,2,3,4], [1,-3]), [[1,5,18],[58]]);
+  assert_equal(poly_div([0], [1,-3]), [[0],[0]]);
 }
 test_poly_div();
 
diff --git a/transforms.scad b/transforms.scad
index be73fbc..e9cfafd 100644
--- a/transforms.scad
+++ b/transforms.scad
@@ -306,11 +306,11 @@ function up(z=0,p=undef) = move([0,0,z],p=p);
 //   * Called as a function with a `p` argument containing a list of points, returns the list of rotated points.
 //   * Called as a function with a [bezier patch](beziers.scad) in the `p` argument, returns the rotated patch.
 //   * Called as a function with a [VNF structure](vnf.scad) in the `p` argument, returns the rotated VNF.
-//   * Called as a function without a `p` argument, and `planar` is true, returns the affine2d rotational matrix.
+//   * Called as a function without a `p` argument, and `planar` is true, returns the affine2d rotational matrix.  Requires that `a` is a finite scalar.
 //   * Called as a function without a `p` argument, and `planar` is false, returns the affine3d rotational matrix.
 //
 // Arguments:
-//   a = Scalar angle or vector of XYZ rotation angles to rotate by, in degrees.
+//   a = Scalar angle or vector of XYZ rotation angles to rotate by, in degrees.  If `planar` is true and `p` is not given, then `a` must be a finite scalar.  Default: `0`
 //   v = vector for the axis of rotation.  Default: [0,0,1] or UP
 //   cp = centerpoint to rotate around. Default: [0,0,0]
 //   from = Starting vector for vector-based rotations.
@@ -343,16 +343,21 @@ module rot(a=0, v=undef, cp=undef, from=undef, to=undef, reverse=false)
 
 function rot(a=0, v, cp, from, to, reverse=false, planar=false, p, _m) =
     assert(is_undef(from)==is_undef(to), "from and to must be specified together.")
+    assert(is_undef(from) || is_vector(from, zero=false), "'from' must be a non-zero vector.")
+    assert(is_undef(to) || is_vector(to, zero=false), "'to' must be a non-zero vector.")
+    assert(is_undef(v) || is_vector(v, zero=false), "'v' must be a non-zero vector.")
+    assert(is_undef(cp) || is_vector(cp), "'cp' must be a vector.")
+    assert(is_finite(a) || is_vector(a), "'a' must be a finite scalar or a vector.")
+    assert(is_bool(reverse))
+    assert(is_bool(planar))
     is_undef(p)? (
         planar? let(
+            check = assert(is_num(a)),
             cp = is_undef(cp)? cp : point2d(cp),
             m1 = is_undef(from)? affine2d_zrot(a) :
-                assert(is_vector(from))
-                assert(!approx(norm(from),0))
-                assert(approx(point3d(from).z, 0))
-                assert(is_vector(to))
-                assert(!approx(norm(to),0))
-                assert(approx(point3d(to).z, 0))
+                assert(a==0, "'from' and 'to' cannot be used with 'a' when 'planar' is true.")
+                assert(approx(point3d(from).z, 0), "'from' must be a 2D vector when 'planar' is true.")
+                assert(approx(point3d(to).z, 0), "'to' must be a 2D vector when 'planar' is true.")
                 affine2d_zrot(
                     vang(point2d(to)) -
                     vang(point2d(from))
@@ -364,13 +369,10 @@ function rot(a=0, v, cp, from, to, reverse=false, planar=false, p, _m) =
             to = is_undef(to)? undef : point3d(to),
             cp = is_undef(cp)? undef : point3d(cp),
             m1 = !is_undef(from)? (
-                    assert(is_vector(from))
-                    assert(!approx(norm(from),0))
-                    assert(is_vector(to))
-                    assert(!approx(norm(to),0))
+                    assert(is_num(a))
                     affine3d_rot_from_to(from,to) * affine3d_zrot(a)
                 ) :
-                !is_undef(v)? affine3d_rot_by_axis(v,a) :
+                !is_undef(v)? assert(is_num(a)) affine3d_rot_by_axis(v,a) :
                 is_num(a)? affine3d_zrot(a) :
                 affine3d_zrot(a.z) * affine3d_yrot(a.y) * affine3d_xrot(a.x),
             m2 = is_undef(cp)? m1 : (move(cp) * m1 * move(-cp)),
@@ -558,12 +560,12 @@ function scale(v=1, p=undef) =
         len(v)==2? affine2d_scale(v) : affine3d_scale(point3d(v))
     ) : (
         assert(is_list(p))
-        is_num(p.x)? vmul(p,v) :
+        is_vector(p)? ( len(p)==2? vmul(p,point2d(v)) : vmul(p,point3d(v,1)) ) :
         is_vnf(p)? let(inv=product([for (x=v) x<0? -1 : 1])) [
-            scale(v=v,p=p.x),
-            inv>=0? p.y : [for (l=p.y) reverse(l)]
+            scale(v=v, p=p[0]),
+            inv>=0? p[1] : [for (l=p[1]) reverse(l)]
         ] :
-        [for (l=p) is_vector(l)? vmul(l,v) : scale(v=v, p=l)]
+        [ for (pp=p) scale(v=v, p=pp) ]
     );
 
 
diff --git a/vectors.scad b/vectors.scad
index 80909b2..1749981 100644
--- a/vectors.scad
+++ b/vectors.scad
@@ -64,13 +64,10 @@ function vang(v) =
 // Example:
 //   vmul([3,4,5], [8,7,6]);  // Returns [24, 28, 30]
 function vmul(v1, v2) = 
-// this thighter check can be done yet because it would break other codes in the library
-//    assert( is_vector(v1) && is_vector(v2,len(v1)), "Incompatible vectors")
-    assert( is_vector(v1) && is_vector(v2), "Invalid vector(s)")
+    assert( is_vector(v1) && is_vector(v2,len(v1)), "Incompatible vectors")
     [for (i = [0:1:len(v1)-1]) v1[i]*v2[i]];
     
 
-
 // Function: vdiv()
 // Description:
 //   Element-wise vector division.  Divides each element of vector `v1` by
diff --git a/version.scad b/version.scad
index 6e03aa9..c2e6cda 100644
--- a/version.scad
+++ b/version.scad
@@ -8,7 +8,7 @@
 //////////////////////////////////////////////////////////////////////
 
 
-BOSL_VERSION = [2,0,401];
+BOSL_VERSION = [2,0,409];
 
 
 // Section: BOSL Library Version Functions
diff --git a/vnf.scad b/vnf.scad
index b72902e..05fd82f 100644
--- a/vnf.scad
+++ b/vnf.scad
@@ -403,6 +403,16 @@ function _triangulate_planar_convex_polygons(polys) =
         outtris = concat(tris, newtris, newtris2)
     ) outtris;
 
+//**
+// this function may produce degenerate triangles:
+//    _triangulate_planar_convex_polygons([ [for(i=[0:1]) [i,i],
+//                                           [1,-1], [-1,-1],
+//                                           for(i=[-1:0]) [i,i] ] ] )
+//    == [[[-1, -1], [ 0,  0], [0,  0]]
+//        [[-1, -1], [-1, -1], [0,  0]]
+//        [[ 1, -1], [-1, -1], [0,  0]]
+//        [[ 0,  0], [ 1,  1], [1, -1]] ]
+//
 
 // Function: vnf_bend()
 // Usage:
@@ -647,7 +657,7 @@ function vnf_validate(vnf, show_warns=true, check_isects=false) =
         nonplanars = unique([
             for (face = faces) let(
                 faceverts = [for (k=face) varr[k]]
-            ) if (!points_are_coplanar(faceverts)) [
+            ) if (!coplanar(faceverts)) [
                 "ERROR",
                 "NONPLANAR",
                 "Face vertices are not coplanar",