diff --git a/.github/workflows/docsgen.yml b/.github/workflows/docsgen.yml
index 8ab135c..911fac2 100644
--- a/.github/workflows/docsgen.yml
+++ b/.github/workflows/docsgen.yml
@@ -41,6 +41,11 @@ jobs:
         cd $GITHUB_WORKSPACE
         ./scripts/genindex.sh
 
+    - name: Generate Cheat Sheet
+      run: |
+        cd $GITHUB_WORKSPACE
+        ./scripts/gencheat.sh
+
     - name: Generating Docs
       env:
         GH_PAT: ${{ secrets.GH_PAT }}
diff --git a/arrays.scad b/arrays.scad
index 835fd65..8f44fe8 100644
--- a/arrays.scad
+++ b/arrays.scad
@@ -19,6 +19,32 @@
 // Section: List Query Operations
 
 
+// Function: is_homogeneous()
+// Usage:
+//   is_homogeneous(list,depth)
+// Description:
+//   Returns true when the list have elements of same type up to the depth `depth`.
+//   Booleans and numbers are not distinguinshed as of distinct types. 
+// Arguments:
+//   list = the list to check
+//   depth = the lowest level the check is done
+// Example:
+//   is_homogeneous( [[1,["a"]], [2,["b"]]] )     // Returns true
+//   is_homogeneous( [[1,["a"]], [2,[true]]] )    // Returns false
+//   is_homogeneous( [[1,["a"]], [2,[true]]], 1 ) // Returns true
+//   is_homogeneous( [[1,["a"]], [2,[true]]], 2 ) // Returns false
+//   is_homogeneous( [[1,["a"]], [true,["b"]]] )  // Returns true
+function is_homogeneous(l, depth) =
+    !is_list(l) || l==[] ? false :
+    let( l0=l[0] )
+    [] == [for(i=[1:len(l)-1]) if( ! _same_type(l[i],l0, depth+1) )  0 ];
+                 
+function _same_type(a,b, depth) = 
+    (depth==0) || (a>=b) || (a==b) || (a<=b)  
+    ||  ( is_list(a) && is_list(b) && len(a)==len(b) 
+          && []==[for(i=idx(a)) if( ! _same_type(a[i],b[i],depth-1) )0] ); 
+  
+
 // Function: select()
 // Description:
 //   Returns a portion of a list, wrapping around past the beginning, if end<start. 
@@ -597,6 +623,23 @@ function list_fit(array, length, fill) =
 
 // Section: List Shuffling and Sorting
 
+
+// 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
+// this allows imax=INF as a bound to numerical lists
+function _valid_idx(idx,imin,imax) =
+    is_undef(idx) 
+    || ( is_finite(idx)  
+         && ( is_undef(imin) || idx>=imin ) 
+         && ( is_undef(imax) || idx< imax ) )
+    || ( is_list(idx)  
+         && ( is_undef(imin) || min(idx)>=imin ) 
+         && ( is_undef(imax) || max(idx)< imax ) )
+    || ( is_range(idx) 
+         && ( is_undef(imin) || (idx[1]>0 && idx[0]>=imin ) || (idx[1]<0 && idx[0]<=imax ) )
+         && ( is_undef(imax) || (idx[1]>0 && idx[2]<=imax ) || (idx[1]<0 && idx[2]>=imin ) ) );
+    
+
 // Function: shuffle()
 // Description:
 //   Shuffles the input list into random order.
@@ -611,7 +654,8 @@ function shuffle(list) =
     concat(shuffle(left), shuffle(right));
 
 
-// Sort a vector of scalar values
+// Sort a vector of scalar values with the native comparison operator
+// all elements should have the same type.
 function _sort_scalars(arr) =
     len(arr)<=1 ? arr : 
     let(
@@ -623,105 +667,67 @@ function _sort_scalars(arr) =
     concat( _sort_scalars(lesser), equal, _sort_scalars(greater) );
 
 
-// Sort a vector of vectors based on the first entry only of each vector
-function _sort_vectors1(arr) =
-    len(arr)<=1 ? arr :
-    !(len(arr)>0) ? [] : 
+// lexical sort of a homogeneous list of vectors 
+// uses native comparison operator
+function _sort_vectors(arr, _i=0) =
+    len(arr)<=1 || _i>=len(arr[0]) ? arr :
     let(
-        pivot   = arr[floor(len(arr)/2)],
-        lesser  = [ for (y = arr) if (y[0]  < pivot[0]) y ],
-        equal   = [ for (y = arr) if (y[0] == pivot[0]) y ],
-        greater = [ for (y = arr) if (y[0]  > pivot[0]) y ]
-    ) 
-    concat( _sort_vectors1(lesser), equal, _sort_vectors1(greater) );
+        pivot   = arr[floor(len(arr)/2)][_i],
+        lesser  = [ for (entry=arr) if (entry[_i]  < pivot ) entry ],
+        equal   = [ for (entry=arr) if (entry[_i] == pivot ) entry ],
+        greater = [ for (entry=arr) if (entry[_i]  > pivot ) entry ]
+      )
+    concat(
+        _sort_vectors(lesser,  _i   ), 
+        _sort_vectors(equal,   _i+1 ), 
+        _sort_vectors(greater, _i ) );
+        
 
-
-// Sort a vector of vectors based on the first two entries of each vector
-// Lexicographic order, remaining entries of vector ignored
-function _sort_vectors2(arr) =
-    len(arr)<=1 ? arr :
-    !(len(arr)>0) ? [] : 
+// lexical sort of a homogeneous list of vectors by the vector components with indices in idxlist
+// all idxlist indices should be in the range of the vector dimensions
+// idxlist must be undef or a simple list of numbers
+// uses native comparison operator
+function _sort_vectors(arr, idxlist, _i=0) =
+    len(arr)<=1 || ( is_list(idxlist) && _i>=len(idxlist) ) || _i>=len(arr[0])  ? arr :
     let(
-        pivot   = arr[floor(len(arr)/2)],
-        lesser  = [ for (y = arr) if (y[0] < pivot[0] || (y[0]==pivot[0] && y[1]<pivot[1])) y ],
-        equal   = [ for (y = arr) if (y[0] == pivot[0] && y[1]==pivot[1]) y ],
-        greater = [ for (y = arr) if (y[0] > pivot[0] || (y[0]==pivot[0] && y[1]>pivot[1])) y ]
-    ) 
-    concat( _sort_vectors2(lesser), equal, _sort_vectors2(greater) );
+        k = is_list(idxlist) ? idxlist[_i] : _i,
+        pivot   = arr[floor(len(arr)/2)][k],
+        lesser  = [ for (entry=arr) if (entry[k]  < pivot ) entry ],
+        equal   = [ for (entry=arr) if (entry[k] == pivot ) entry ],
+        greater = [ for (entry=arr) if (entry[k]  > pivot ) entry ]
+      )
+    concat(
+        _sort_vectors(lesser,  idxlist, _i  ), 
+        _sort_vectors(equal,   idxlist, _i+1), 
+        _sort_vectors(greater, idxlist, _i  ) );
+        
 
-
-// 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) =
-    len(arr)<=1 ? arr : let(
-        pivot   = arr[floor(len(arr)/2)],
-        lesser  = [ for (y = arr) 
-                      if ( y[0] < pivot[0] 
-                           || ( y[0]==pivot[0] 
-                                && ( y[1]<pivot[1] 
-                                     || ( y[1]==pivot[1] 
-                                          && y[2]<pivot[2] ))))
-                    y ],
-        equal = [ for (y = arr) 
-                    if ( y[0] == pivot[0] 
-                         && y[1]== pivot[1] 
-                         && y[2]==pivot[2] )
-                  y ],
-        greater = [ for (y = arr) 
-                      if ( y[0] > pivot[0] 
-                           || ( y[0]==pivot[0] 
-                                && ( y[1] > pivot[1] 
-                                     || ( y[1]==pivot[1] 
-                                          && y[2] > pivot[2] )))) 
-                    y ]
-    ) 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) =
-    len(arr)<=1 ? arr : let(
-        pivot = arr[floor(len(arr)/2)],
-        lesser = [  for (y = arr) 
-                      if ( y[0] < pivot[0] 
-                           || ( y[0]==pivot[0] 
-                                && ( y[1]<pivot[1] 
-                                     || ( y[1]==pivot[1] 
-                                         && ( y[2]<pivot[2] 
-                                              || ( y[2]==pivot[2] 
-                                                   && y[3]<pivot[3] ))))))
-                      y ],
-        equal = [ for (y = arr) 
-                    if (  y[0] == pivot[0] 
-                          && y[1] == pivot[1] 
-                          && y[2] == pivot[2] 
-                          && y[3] == pivot[3] ) 
-                  y  ],
-        greater = [ for (y = arr) 
-                      if ( y[0] > pivot[0] 
-                           || ( y[0]==pivot[0] 
-                                && ( y[1]>pivot[1] 
-                                     || ( y[1]==pivot[1] 
-                                          && ( y[2]>pivot[2] 
-                                               || ( y[2]==pivot[2] 
-                                                    && 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) =
+// sorting using compare_vals(); returns indexed list when `indexed==true`
+function _sort_general(arr, idx=undef, indexed=false) =
     (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);
+    ! indexed && is_undef(idx)
+    ? _lexical_sort(arr)
+    : let( arrind = _indexed_sort(enumerate(arr,idx)) )
+      indexed 
+      ? arrind
+      : [for(i=arrind) arr[i]];
       
+// lexical sort using compare_vals()
+function _lexical_sort(arr) = 
+    arr==[] ? [] : len(arr)==1? arr : 
+    let( pivot = arr[floor(len(arr)/2)] )
+    let(
+        lesser  = [ for (entry=arr) if (compare_vals(entry, pivot) <0 ) entry ],
+        equal   = [ for (entry=arr) if (compare_vals(entry, pivot)==0 ) entry ],
+        greater = [ for (entry=arr) if (compare_vals(entry, pivot) >0 ) entry ]
+      )
+    concat(_lexical_sort(lesser), equal, _lexical_sort(greater));
+
+
 // 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]] :
+    arrind==[] ? [] : len(arrind)==1? [arrind[0][0]] : 
     let( pivot = arrind[floor(len(arrind)/2)][1] )
     let(
         lesser  = [ for (entry=arrind) if (compare_vals(entry[1], pivot) <0 ) entry ],
@@ -731,52 +737,46 @@ function _indexed_sort(arrind) =
     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()
 // Usage:
 //   sort(list, [idx])
 // Description:
-//   Sorts the given list using `compare_vals()`, sorting in lexicographic order, with types ordered according to
+//   Sorts the given list in lexicographic order. If the input is a homogeneous simple list or a homogeneous 
+//   list of vectors (see function is_homogeneous), the sorting method uses the native comparison operator and is faster. 
+//   When sorting non homogeneous list the elements are compared with `compare_vals`, with types ordered according to
 //   `undef < boolean < number < string < list`.  Comparison of lists is recursive. 
-//   If the list is a list of vectors whose length is from 1 to 4 and the `idx` parameter is not passed, then 
-//   `sort` uses a much more efficient method for comparisons and will run much faster.  In this case, all entries
-//   in the data are compared using the native comparison operator, so comparisons between types will fail.  
+//   When comparing vectors, homogeneous or not, the parameter `idx` may be used to select the components to compare.
+//   Note that homogeneous lists of vectors may contain mixed types provided that for any two list elements
+//   list[i] and list[j] satisfies  type(list[i][k])==type(list[j][k]) for all k. 
+//   Strings are allowed as any list element and are compared with the native operators although no substring
+//   comparison is possible.  
 // Arguments:
 //   list = The list to sort.
 //   idx = If given, do the comparison based just on the specified index, range or list of indices.  
-// Example:
-//   l = [45,2,16,37,8,3,9,23,89,12,34];
-//   sorted = sort(l);  // Returns [2,3,8,9,12,16,23,34,37,45,89]
-function sort(list, idx=undef) =
+// Example: 
+//   // Homogeneous lists
+//   l1 = [45,2,16,37,8,3,9,23,89,12,34];
+//   sorted1 = sort(l1);  // Returns [2,3,8,9,12,16,23,34,37,45,89]
+//   l2 = [["oat",0], ["cat",1], ["bat",3], ["bat",2], ["fat",3]];
+//   sorted2 = sort(l2); // Returns: [["bat",2],["bat",3],["cat",1],["fat",3],["oat",0]]
+//   // Non-homegenous list
+//   l3 = [[4,0],[7],[3,9],20,[4],[3,1],[8]];
+//   sorted3 = sort(l3); // Returns: [20,[3,1],[3,9],[4],[4,0],[7],[8]]
+function sort(list, idx=undef) = 
     !is_list(list) || len(list)<=1 ? 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);
-
+    is_homogeneous(list,1)
+    ?   let(size = array_dim(list[0]))
+        size==0 ?         _sort_scalars(list)
+        : len(size)!=1 ?  _sort_general(list,idx)  
+        : is_undef(idx) ? _sort_vectors(list)
+        : assert( _valid_idx(idx) , "Invalid indices.")
+          _sort_vectors(list,[for(i=idx) i])        
+    : _sort_general(list,idx);
+        
 
 // Function: sortidx()
 // Description:
-//   Given a list, calculates the sort order of the list, and returns
+//   Given a list, sort it as function `sort()`, and returns
 //   a list of indexes into the original list in that sorted order.
 //   If you iterate the returned list in order, and use the list items
 //   to index into the original list, you will be iterating the original
@@ -795,31 +795,28 @@ function sort(list, idx=undef) =
 //   idxs1 = sortidx(lst, idx=1); // Returns: [3,0,2,1]
 //   idxs2 = sortidx(lst, idx=0); // Returns: [1,2,0,3]
 //   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( _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)))
-              :  0
-    )
-    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
-        _sort_general(list,idx);
-
-
-// sort() does not accept strings but sortidx does; isn't inconsistent ?
-
+function sortidx(list, idx=undef) = 
+    !is_list(list) || len(list)<=1 ? list :
+    is_homogeneous(list,1)
+    ?   let( 
+            size = array_dim(list[0]),
+            aug  = ! (size==0 || len(size)==1) ? 0 // for general sorting
+                   : [for(i=[0:len(list)-1]) concat(i,list[i])], // for scalar or vector sorting
+            lidx = size==0? [1] :                                // scalar sorting
+                   len(size)==1 
+                   ? is_undef(idx) ? [for(i=[0:len(list[0])-1]) i+1] // vector sorting
+                                   : [for(i=idx) i+1]                // vector sorting
+                   : 0   // just to signal
+            )
+        assert( ! ( size==0 && is_def(idx) ), 
+                "The specification of `idx` is incompatible with scalar sorting." ) 
+        assert( _valid_idx(idx) , "Invalid indices." ) 
+        lidx!=0
+        ?   let( lsort = _sort_vectors(aug,lidx) )
+            [for(li=lsort) li[0] ]
+        :   _sort_general(list,idx,indexed=true)
+    : _sort_general(list,idx,indexed=true);
+        
 
 // Function: unique()
 // Usage:
@@ -1202,6 +1199,52 @@ function zip(vecs, v2, v3, fit=false, fill=undef) =
     :   [for(i=[0:1:minlen-1]) [for(v=vecs) for(x=v[i]) x] ];
 
 
+// Function: block_matrix()
+// Usage:
+//    block_matrix([[M11, M12,...],[M21, M22,...], ... ])
+// Description:
+//    Create a block matrix by supplying a matrix of matrices, which will
+//    be combined into one unified matrix.  Every matrix in one row
+//    must have the same height, and the combined width of the matrices
+//    in each row must be equal.  
+function block_matrix(M) =
+    let(
+        bigM = [for(bigrow = M) each zip(bigrow)],
+        len0=len(bigM[0]),
+        badrows = [for(row=bigM) if (len(row)!=len0) 1]
+    )
+    assert(badrows==[], "Inconsistent or invalid input")
+    bigM;
+
+// Function: diagonal_matrix()
+// Usage:
+//   diagonal_matrix(diag, [offdiag])
+// Description:
+//   Creates a square matrix with the items in the list `diag` on
+//   its diagonal.  The off diagonal entries are set to offdiag,
+//   which is zero by default. 
+function diagonal_matrix(diag,offdiag=0) =
+  [for(i=[0:1:len(diag)-1]) [for(j=[0:len(diag)-1]) i==j?diag[i] : offdiag]];
+
+
+// Function: submatrix_set()
+// Usage: submatrix_set(M,A,[m],[n])
+// Description:
+//    Sets a submatrix of M equal to the matrix A.  By default the top left corner of M is set to A, but
+//    you can specify offset coordinates m and n.  If A (as adjusted by m and n) extends beyond the bounds
+//    of M then the extra entries are ignored.  You can pass in A=[[]], a null matrix, and M will be
+//    returned unchanged.  Note that the input M need not be rectangular in shape.  
+function submatrix_set(M,A,m=0,n=0) =
+    assert(is_list(M))
+    assert(is_list(A))
+    let( badrows = [for(i=idx(A)) if (!is_list(A[i])) i])
+    assert(badrows==[], str("Input submatrix malformed rows: ",badrows))
+    [for(i=[0:1:len(M)-1])
+        assert(is_list(M[i]), str("Row ",i," of input matrix is not a list"))
+        [for(j=[0:1:len(M[i])-1]) 
+            i>=m && i <len(A)+m && j>=n && j<len(A[0])+n ? A[i-m][j-n] : M[i][j]]];
+
+
 // Function: array_group()
 // Description:
 //   Takes a flat array of values, and groups items in sets of `cnt` length.
@@ -1241,16 +1284,22 @@ function full_flatten(l) = [for(a=l) if(is_list(a)) (each full_flatten(a)) else
 // Internal.  Not exposed.
 function _array_dim_recurse(v) =
     !is_list(v[0])
-    ?   sum( [for(entry=v) is_list(entry) ? 1 : 0] ) == 0 ? [] : [undef]
+    ?   len( [for(entry=v) if(!is_list(entry)) 0] ) == 0 ? [] : [undef]
     :   let(
-          firstlen = len(v[0]),
-          first = sum( [for(entry = v) len(entry) == firstlen  ? 0 : 1]   ) == 0 ? firstlen : undef,
+          firstlen = is_list(v[0]) ? len(v[0]): undef,
+          first = len( [for(entry = v) if(! is_list(entry) || (len(entry) != firstlen)) 0  ]   ) == 0 ? firstlen : undef,
           leveldown = flatten(v)
         ) 
         is_list(leveldown[0])
         ?  concat([first],_array_dim_recurse(leveldown))
         : [first];
 
+function _array_dim_recurse(v) =
+    let( alen = [for(vi=v) is_list(vi) ? len(vi): -1] )
+    v==[] || max(alen)==-1 ? [] :
+    let( add = max(alen)!=min(alen) ? undef : alen[0] ) 
+    concat( add, _array_dim_recurse(flatten(v)));
+
 
 // Function: array_dim()
 // Usage:
@@ -1281,15 +1330,19 @@ function array_dim(v, depth=undef) =
         ?  len(v)
         :  let( dimlist = _array_dim_recurse(v))
            (depth > len(dimlist))? 0 : dimlist[depth-1] ;
+           
+           
 
 // This function may return 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)
+// Usage:
+//    transpose(arr, [reverse])
+// Description:
+//    Returns the transpose of the given input array.  The input should be a list of lists that are
+//    all the same length.  If you give a vector then transpose returns it unchanged.  
+//    When reverse=true, the transpose is done across to the secondary diagonal.  (See example below.)
 //    By default, reverse=false.
 // Example:
 //   arr = [
@@ -1329,19 +1382,19 @@ function array_dim(v, depth=undef) =
 //   //  ["h", "e", "b"],
 //   //  ["g", "d", "a"]
 //   // ]
-// Example:
+// Example: Transpose on a list of numbers returns the list unchanged
 //   transpose([3,4,5]);  // Returns: [3,4,5]
 function transpose(arr, reverse=false) =
-    assert( is_list(arr) && len(arr)>0, "The array is not a vector neither a matrix." )
+    assert( is_list(arr) && len(arr)>0, "Input to transpose must be a nonempty list.")
     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." )
+    ?   let( len0 = len(arr[0]) )
+        assert([for(a=arr) if(!is_list(a) || len(a)!=len0) 1 ]==[], "Input to transpose has inconsistent row lengths." )
         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 (i=[0:1:len0-1]) 
+              [ for (j=[0:1:len(arr)-1]) arr[len(arr)-1-j][len0-1-i] ] ] 
+        : [for (i=[0:1:len0-1]) 
               [ for (j=[0:1:len(arr)-1]) arr[j][i] ] ] 
-    :  assert( is_vector(arr), "The array is not a vector neither a matrix." )
+    :  assert( is_vector(arr), "Input to transpose must be a vector or list of lists.")
            arr;
 
 
diff --git a/attachments.scad b/attachments.scad
index 48ed522..ae93412 100644
--- a/attachments.scad
+++ b/attachments.scad
@@ -461,13 +461,16 @@ function find_anchor(anchor, geom) =
             rpts = apply(rot(from=anchor, to=RIGHT) * move(point3d(-cp)), points),
             hits = [
                 for (face = faces) let(
-                    verts = select(rpts, face)
+                    verts = select(rpts, face),
+                    xs = subindex(verts,0),
+                    ys = subindex(verts,1),
+                    zs = subindex(verts,2)
                 ) if (
-                    max(subindex(verts,0)) >= -eps &&
-                    max(subindex(verts,1)) >= -eps &&
-                    max(subindex(verts,2)) >= -eps &&
-                    min(subindex(verts,1)) <=  eps &&
-                    min(subindex(verts,2)) <=  eps
+                    max(xs) >= -eps &&
+                    max(ys) >= -eps &&
+                    max(zs) >= -eps &&
+                    min(ys) <=  eps &&
+                    min(zs) <=  eps
                 ) let(
                     poly = select(points, face),
                     pt = polygon_line_intersection(poly, [cp,cp+anchor], bounded=[true,false], eps=eps)
diff --git a/beziers.scad b/beziers.scad
index 13a74ab..a8484b6 100644
--- a/beziers.scad
+++ b/beziers.scad
@@ -356,7 +356,7 @@ function bezier_segment_length(curve, start_u=0, end_u=1, max_deflect=0.01) =
 
 // Function: fillet3pts()
 // Usage:
-//   fillet3pts(p0, p1, p2, r);
+//   fillet3pts(p0, p1, p2, r|d);
 // Description:
 //   Takes three points, defining two line segments, and works out the
 //   cubic (degree 3) bezier segment (and surrounding control points)
@@ -368,7 +368,8 @@ function bezier_segment_length(curve, start_u=0, end_u=1, max_deflect=0.01) =
 //   p1 = The middle point.
 //   p2 = The ending point.
 //   r = The radius of the fillet/rounding.
-//   maxerr = Max amount bezier curve should diverge from actual radius curve.  Default: 0.1
+//   d = The diameter of the fillet/rounding.
+//   maxerr = Max amount bezier curve should diverge from actual curve.  Default: 0.1
 // Example(2D):
 //   p0 = [40, 0];
 //   p1 = [0, 0];
@@ -376,7 +377,8 @@ function bezier_segment_length(curve, start_u=0, end_u=1, max_deflect=0.01) =
 //   trace_polyline([p0,p1,p2], showpts=true, size=0.5, color="green");
 //   fbez = fillet3pts(p0,p1,p2, 10);
 //   trace_bezier(slice(fbez, 1, -2), size=1);
-function fillet3pts(p0, p1, p2, r, maxerr=0.1, w=0.5, dw=0.25) = let(
+function fillet3pts(p0, p1, p2, r, d, maxerr=0.1, w=0.5, dw=0.25) = let(
+        r = get_radius(r=r,d=d),
         v0 = unit(p0-p1),
         v1 = unit(p2-p1),
         midv = unit((v0+v1)/2),
@@ -391,8 +393,8 @@ function fillet3pts(p0, p1, p2, r, maxerr=0.1, w=0.5, dw=0.25) = let(
         bp = bezier_points([tp0, cp0, cp1, tp1], 0.5),
         tdist = norm(cp-bp)
     ) (abs(tdist-cpr) <= maxerr)? [tp0, tp0, cp0, cp1, tp1, tp1] :
-        (tdist<cpr)? fillet3pts(p0, p1, p2, r, maxerr=maxerr, w=w+dw, dw=dw/2) :
-        fillet3pts(p0, p1, p2, r, maxerr=maxerr, w=w-dw, dw=dw/2);
+        (tdist<cpr)? fillet3pts(p0, p1, p2, r=r, maxerr=maxerr, w=w+dw, dw=dw/2) :
+        fillet3pts(p0, p1, p2, r=r, maxerr=maxerr, w=w-dw, dw=dw/2);
 
 
 
@@ -613,7 +615,7 @@ function fillet_path(pts, fillet, maxerr=0.1) = concat(
             p1 = pts[p],
             p0 = (pts[p-1]+p1)/2,
             p2 = (pts[p+1]+p1)/2
-        ) for (pt = fillet3pts(p0, p1, p2, fillet, maxerr=maxerr)) pt
+        ) for (pt = fillet3pts(p0, p1, p2, r=fillet, maxerr=maxerr)) pt
     ],
     [pts[len(pts)-1], pts[len(pts)-1]]
 );
diff --git a/common.scad b/common.scad
index db4f3bb..570ff8e 100644
--- a/common.scad
+++ b/common.scad
@@ -353,6 +353,7 @@ function segs(r) =
 // Arguments:
 //   $children = number of children the module has.  
 module no_children(count) {
+  assert($children==0, "Module no_children() does not support child modules");
   assert(count==0, str("Module ",parent_module(1),"() does not support child modules"));
 }
 
@@ -377,6 +378,7 @@ function _valstr(x) =
 //   expected = The value that was expected.
 //   info = Extra info to print out to make the error clearer.
 module assert_approx(got, expected, info) {
+    no_children($children);
     if (!approx(got, expected)) {
         echo();
         echo(str("EXPECT: ", _valstr(expected)));
@@ -404,6 +406,7 @@ module assert_approx(got, expected, info) {
 //   expected = The value that was expected.
 //   info = Extra info to print out to make the error clearer.
 module assert_equal(got, expected, info) {
+    no_children($children);
     if (got != expected || (is_nan(got) && is_nan(expected))) {
         echo();
         echo(str("EXPECT: ", _valstr(expected)));
diff --git a/distributors.scad b/distributors.scad
index 787693e..6176165 100644
--- a/distributors.scad
+++ b/distributors.scad
@@ -686,7 +686,7 @@ module rot_copies(rots=[], v=undef, cp=[0,0,0], n=undef, sa=0, offset=0, delta=[
 //   cp = Centerpoint to rotate around.
 //   n = Optional number of evenly distributed copies to be rotated around the ring.
 //   sa = Starting angle, in degrees.  For use with `n`.  Angle is in degrees counter-clockwise from Y+, when facing the origin from X+.  First unrotated copy is placed at that angle.
-//   r = Radius to move children back, away from cp, before rotating.  Makes rings of copies.
+//   r = Radius to move children back (Y+), away from cp, before rotating.  Makes rings of copies.
 //   subrot = If false, don't sub-rotate children as they are copied around the ring.
 //
 // Side Effects:
@@ -743,7 +743,7 @@ module xrot_copies(rots=[], cp=[0,0,0], n=undef, sa=0, r=0, subrot=true)
 //   cp = Centerpoint to rotate around.
 //   n = Optional number of evenly distributed copies to be rotated around the ring.
 //   sa = Starting angle, in degrees.  For use with `n`.  Angle is in degrees counter-clockwise from X-, when facing the origin from Y+.
-//   r = Radius to move children left, away from cp, before rotating.  Makes rings of copies.
+//   r = Radius to move children left (X-), away from cp, before rotating.  Makes rings of copies.
 //   subrot = If false, don't sub-rotate children as they are copied around the ring.
 //
 // Side Effects:
@@ -800,7 +800,7 @@ module yrot_copies(rots=[], cp=[0,0,0], n=undef, sa=0, r=0, subrot=true)
 //   cp = Centerpoint to rotate around.  Default: [0,0,0]
 //   n = Optional number of evenly distributed copies to be rotated around the ring.
 //   sa = Starting angle, in degrees.  For use with `n`.  Angle is in degrees counter-clockwise from X+, when facing the origin from Z+.  Default: 0
-//   r = Radius to move children right, away from cp, before rotating.  Makes rings of copies.  Default: 0
+//   r = Radius to move children right (X+), away from cp, before rotating.  Makes rings of copies.  Default: 0
 //   subrot = If false, don't sub-rotate children as they are copied around the ring.  Default: true
 //
 // Side Effects:
diff --git a/geometry.scad b/geometry.scad
index dc86187..0ab9f73 100644
--- a/geometry.scad
+++ b/geometry.scad
@@ -1007,11 +1007,20 @@ function closest_point_on_plane(plane, 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( 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 ];
+    let(
+        l0 = line[0],  // Ray start point
+        u = line[1] - l0,  // Ray direction vector
+        n = plane_normal(plane),
+        p0 = n * plane[3],  // A point on the plane
+        w = l0 - p0  // Vector from plane point to ray start
+    ) approx(n*u, 0, eps=eps) ? (
+        // Line is parallel to plane.
+        approx(n*w, 0, eps=eps)
+          ? [line, undef] // Line is on the plane.
+          : undef // Line never intersects the plane.
+    ) : let(
+        t = (-n * w) / (n * u)  // Distance ratio along ray
+    ) [ l0 + u*t, t ];
 
 
 // Function: plane_line_angle()
@@ -1098,8 +1107,8 @@ function polygon_line_intersection(poly, line, bounded=false, eps=EPSILON) =
             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),
+            poly2d = clockwise_polygon(project_plane(poly, plane)),
+            line2d = project_plane([lp1,lp2], plane),
             parts = split_path_at_region_crossings(line2d, [poly2d], closed=false),
             inside = [for (part = parts)
                           if (point_in_polygon(mean(part), poly2d)>0) part
@@ -1107,7 +1116,7 @@ function polygon_line_intersection(poly, line, bounded=false, eps=EPSILON) =
         ) 
         !inside? undef :
         let(
-            isegs = [for (seg = inside) lift_plane(seg, p1, p2, p3) ]
+            isegs = [for (seg = inside) lift_plane(seg, plane) ]
         ) 
         isegs
     ) 
@@ -1264,7 +1273,6 @@ function find_circle_2tangents(pt1, pt2, pt3, r, d, tangents=false) =
             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)
             ) 
diff --git a/masks.scad b/masks.scad
index 894e50d..b421e7d 100644
--- a/masks.scad
+++ b/masks.scad
@@ -368,8 +368,8 @@ module chamfer_hole_mask(r=undef, d=undef, chamfer=0.25, ang=45, from_end=false,
 
 // Module: rounding_mask()
 // Usage:
-//   rounding_mask(l|h, r)
-//   rounding_mask(l|h, r1, r2)
+//   rounding_mask(l|h, r|d)
+//   rounding_mask(l|h, r1|d1, r2|d2)
 // Description:
 //   Creates a shape that can be used to round a vertical 90 degree edge.
 //   Difference it from the object to be rounded.  The center of the mask
@@ -379,6 +379,9 @@ module chamfer_hole_mask(r=undef, d=undef, chamfer=0.25, ang=45, from_end=false,
 //   r = Radius of the rounding.
 //   r1 = Bottom radius of rounding.
 //   r2 = Top radius of rounding.
+//   d = Diameter of the rounding.
+//   d1 = Bottom diameter of rounding.
+//   d2 = Top diameter of rounding.
 //   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`
@@ -409,11 +412,11 @@ module chamfer_hole_mask(r=undef, d=undef, chamfer=0.25, ang=45, from_end=false,
 //               rounding_mask(l=p.x, r=25, spin=45, orient=RIGHT);
 //       }
 //   }
-module rounding_mask(l=undef, r=undef, r1=undef, r2=undef, anchor=CENTER, spin=0, orient=UP, h=undef)
+module rounding_mask(l, r, r1, r2, d, d1, d2, anchor=CENTER, spin=0, orient=UP, h=undef)
 {
     l = first_defined([l, h, 1]);
-    r1 = get_radius(r1=r1, r=r, dflt=1);
-    r2 = get_radius(r1=r2, r=r, dflt=1);
+    r1 = get_radius(r1=r1, r=r, d1=d1, d=d, dflt=1);
+    r2 = get_radius(r1=r2, r=r, d1=d2, d=d, dflt=1);
     sides = quantup(segs(max(r1,r2)),4);
     attachable(anchor,spin,orient, size=[2*r1,2*r1,l], size2=[2*r2,2*r2]) {
         if (r1<r2) {
@@ -440,8 +443,8 @@ module rounding_mask(l=undef, r=undef, r1=undef, r2=undef, anchor=CENTER, spin=0
 
 // Module: rounding_mask_x()
 // Usage:
-//   rounding_mask_x(l, r, [anchor])
-//   rounding_mask_x(l, r1, r2, [anchor])
+//   rounding_mask_x(l, r|d, [anchor])
+//   rounding_mask_x(l, r1|d1, r2|d2, [anchor])
 // Description:
 //   Creates a shape that can be used to round a 90 degree edge oriented
 //   along the X axis.  Difference it from the object to be rounded.
@@ -452,6 +455,9 @@ module rounding_mask(l=undef, r=undef, r1=undef, r2=undef, anchor=CENTER, spin=0
 //   r = Radius of the rounding.
 //   r1 = Left end radius of rounding.
 //   r2 = Right end radius of rounding.
+//   d = Diameter of the rounding.
+//   d1 = Left end diameter of rounding.
+//   d2 = Right end diameter of rounding.
 //   spin = Rotate this many degrees around the Z axis after anchor.  See [spin](attachments.scad#spin).  Default: `0`
 // Example:
 //   difference() {
@@ -463,10 +469,10 @@ module rounding_mask(l=undef, r=undef, r1=undef, r2=undef, anchor=CENTER, spin=0
 //       cube(size=100, center=false);
 //       #rounding_mask_x(l=100, r1=10, r2=30, anchor=LEFT);
 //   }
-module rounding_mask_x(l=1.0, r=undef, r1=undef, r2=undef, anchor=CENTER, spin=0)
+module rounding_mask_x(l=1.0, r, r1, r2, d, d1, d2, anchor=CENTER, spin=0)
 {
     anchor = rot(p=anchor, from=RIGHT, to=TOP);
-    rounding_mask(l=l, r=r, r1=r1, r2=r2, anchor=anchor, spin=spin, orient=RIGHT) {
+    rounding_mask(l=l, r=r, r1=r1, r2=r2, d=d, d1=d1, d2=d2, anchor=anchor, spin=spin, orient=RIGHT) {
         for (i=[0:1:$children-2]) children(i);
         if ($children) children($children-1);
     }
@@ -475,8 +481,8 @@ module rounding_mask_x(l=1.0, r=undef, r1=undef, r2=undef, anchor=CENTER, spin=0
 
 // Module: rounding_mask_y()
 // Usage:
-//   rounding_mask_y(l, r, [anchor])
-//   rounding_mask_y(l, r1, r2, [anchor])
+//   rounding_mask_y(l, r|d, [anchor])
+//   rounding_mask_y(l, r1|d1, r2|d2, [anchor])
 // Description:
 //   Creates a shape that can be used to round a 90 degree edge oriented
 //   along the Y axis.  Difference it from the object to be rounded.
@@ -487,6 +493,9 @@ module rounding_mask_x(l=1.0, r=undef, r1=undef, r2=undef, anchor=CENTER, spin=0
 //   r = Radius of the rounding.
 //   r1 = Front end radius of rounding.
 //   r2 = Back end radius of rounding.
+//   d = Diameter of the rounding.
+//   d1 = Front end diameter of rounding.
+//   d2 = Back end diameter of rounding.
 //   spin = Rotate this many degrees around the Z axis after anchor.  See [spin](attachments.scad#spin).  Default: `0`
 // Example:
 //   difference() {
@@ -498,10 +507,10 @@ module rounding_mask_x(l=1.0, r=undef, r1=undef, r2=undef, anchor=CENTER, spin=0
 //       cube(size=100, center=false);
 //       right(100) #rounding_mask_y(l=100, r1=10, r2=30, anchor=FRONT);
 //   }
-module rounding_mask_y(l=1.0, r=undef, r1=undef, r2=undef, anchor=CENTER, spin=0)
+module rounding_mask_y(l=1.0, r, r1, r2, d, d1, d2, anchor=CENTER, spin=0)
 {
     anchor = rot(p=anchor, from=BACK, to=TOP);
-    rounding_mask(l=l, r=r, r1=r1, r2=r2, anchor=anchor, spin=spin, orient=BACK) {
+    rounding_mask(l=l, r=r, r1=r1, r2=r2, d=d, d1=d1, d2=d2, anchor=anchor, spin=spin, orient=BACK) {
         for (i=[0:1:$children-2]) children(i);
         if ($children) children($children-1);
     }
@@ -510,8 +519,8 @@ module rounding_mask_y(l=1.0, r=undef, r1=undef, r2=undef, anchor=CENTER, spin=0
 
 // Module: rounding_mask_z()
 // Usage:
-//   rounding_mask_z(l, r, [anchor])
-//   rounding_mask_z(l, r1, r2, [anchor])
+//   rounding_mask_z(l, r|d, [anchor])
+//   rounding_mask_z(l, r1|d1, r2|d2, [anchor])
 // Description:
 //   Creates a shape that can be used to round a 90 degree edge oriented
 //   along the Z axis.  Difference it from the object to be rounded.
@@ -522,6 +531,9 @@ module rounding_mask_y(l=1.0, r=undef, r1=undef, r2=undef, anchor=CENTER, spin=0
 //   r = Radius of the rounding.
 //   r1 = Bottom radius of rounding.
 //   r2 = Top radius of rounding.
+//   d = Diameter of the rounding.
+//   d1 = Bottom diameter of rounding.
+//   d2 = Top diameter of rounding.
 //   spin = Rotate this many degrees around the Z axis after anchor.  See [spin](attachments.scad#spin).  Default: `0`
 // Example:
 //   difference() {
@@ -533,9 +545,9 @@ module rounding_mask_y(l=1.0, r=undef, r1=undef, r2=undef, anchor=CENTER, spin=0
 //       cube(size=100, center=false);
 //       #rounding_mask_z(l=100, r1=10, r2=30, anchor=BOTTOM);
 //   }
-module rounding_mask_z(l=1.0, r=undef, r1=undef, r2=undef, anchor=CENTER, spin=0)
+module rounding_mask_z(l=1.0, r, r1, r2, d, d1, d2, anchor=CENTER, spin=0)
 {
-    rounding_mask(l=l, r=r, r1=r1, r2=r2, anchor=anchor, spin=spin, orient=UP) {
+    rounding_mask(l=l, r=r, r1=r1, r2=r2, d=d, d1=d1, d2=d2, anchor=anchor, spin=spin, orient=UP) {
         for (i=[0:1:$children-2]) children(i);
         if ($children) children($children-1);
     }
@@ -544,11 +556,12 @@ module rounding_mask_z(l=1.0, r=undef, r1=undef, r2=undef, anchor=CENTER, spin=0
 
 // Module: rounding()
 // Usage:
-//   rounding(r, size, [edges]) ...
+//   rounding(r|d, size, [edges]) ...
 // Description:
 //   Rounds the edges of a cuboid region containing the given children.
 // Arguments:
 //   r = Radius of the rounding. (Default: 1)
+//   d = Diameter of the rounding. (Default: 1)
 //   size = The size of the rectangular cuboid we want to chamfer.
 //   edges = Edges to round.  See the docs for [`edges()`](edges.scad#edges) to see acceptable values.  Default: All edges.
 //   except_edges = Edges to explicitly NOT round.  See the docs for [`edges()`](edges.scad#edges) to see acceptable values.  Default: No edges.
@@ -560,8 +573,9 @@ module rounding_mask_z(l=1.0, r=undef, r1=undef, r2=undef, anchor=CENTER, spin=0
 //   rounding(r=10, size=[50,50,75], edges=[TOP,FRONT+RIGHT], except_edges=TOP+LEFT, $fn=24) {
 //     cube(size=[50,50,75], center=true);
 //   }
-module rounding(r=1, size=[1,1,1], edges=EDGES_ALL, except_edges=[])
+module rounding(r, size=[1,1,1], d, edges=EDGES_ALL, except_edges=[])
 {
+    r = get_radius(r=r, d=d, dflt=1);
     difference() {
         children();
         difference() {
@@ -574,8 +588,8 @@ module rounding(r=1, size=[1,1,1], edges=EDGES_ALL, except_edges=[])
 
 // Module: rounding_angled_edge_mask()
 // Usage:
-//   rounding_angled_edge_mask(h, r, [ang]);
-//   rounding_angled_edge_mask(h, r1, r2, [ang]);
+//   rounding_angled_edge_mask(h, r|d, [ang]);
+//   rounding_angled_edge_mask(h, r1|d1, r2|d2, [ang]);
 // Description:
 //   Creates a vertical mask that can be used to round the edge where two face meet, at any arbitrary
 //   angle.  Difference it from the object to be rounded.  The center of the mask should align exactly
@@ -585,6 +599,9 @@ module rounding(r=1, size=[1,1,1], edges=EDGES_ALL, except_edges=[])
 //   r = Radius of the rounding.
 //   r1 = Bottom radius of rounding.
 //   r2 = Top radius of rounding.
+//   d = Diameter of the rounding.
+//   d1 = Bottom diameter of rounding.
+//   d2 = Top diameter of rounding.
 //   ang = Angle that the planes meet at.
 //   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`
@@ -599,7 +616,7 @@ module rounding(r=1, size=[1,1,1], edges=EDGES_ALL, except_edges=[])
 //       angle_pie_mask(ang=70, h=50, d=100);
 //       #rounding_angled_edge_mask(h=51, r1=10, r2=25, ang=70, $fn=32);
 //   }
-module rounding_angled_edge_mask(h=1.0, r=undef, r1=undef, r2=undef, ang=90, anchor=CENTER, spin=0, orient=UP)
+module rounding_angled_edge_mask(h=1.0, r, r1, r2, d, d1, d2, ang=90, anchor=CENTER, spin=0, orient=UP)
 {
     function _mask_shape(r) = [
         for (i = [0:1:n]) let (a=90+ang+i*sweep/n) [r*cos(a)+x, r*sin(a)+r],
@@ -609,8 +626,8 @@ module rounding_angled_edge_mask(h=1.0, r=undef, r1=undef, r2=undef, ang=90, anc
     ];
 
     sweep = 180-ang;
-    r1 = get_radius(r1=r1, r=r, dflt=1);
-    r2 = get_radius(r1=r2, r=r, dflt=1);
+    r1 = get_radius(r1=r1, r=r, d1=d1, d=d, dflt=1);
+    r2 = get_radius(r1=r2, r=r, d1=d2, d=d, dflt=1);
     n = ceil(segs(max(r1,r2))*sweep/360);
     x = sin(90-(ang/2))/sin(ang/2) * (r1<r2? r2 : r1);
     if(r1<r2) {
@@ -635,13 +652,14 @@ module rounding_angled_edge_mask(h=1.0, r=undef, r1=undef, r2=undef, ang=90, anc
 
 // Module: rounding_angled_corner_mask()
 // Usage:
-//   rounding_angled_corner_mask(r, ang);
+//   rounding_angled_corner_mask(r|d, ang);
 // Description:
 //   Creates a shape that can be used to round the corner of an angle.
 //   Difference it from the object to be rounded.  The center of the mask
 //   object should align exactly with the point of the corner to be rounded.
 // Arguments:
 //   r = Radius of the rounding.
+//   d = Diameter of the rounding.
 //   ang = Angle between planes that you need to round the corner of.
 //   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`
@@ -656,8 +674,9 @@ module rounding_angled_edge_mask(h=1.0, r=undef, r1=undef, r2=undef, ang=90, anc
 //       }
 //       rounding_angled_edge_mask(h=51, r=20, ang=ang);
 //   }
-module rounding_angled_corner_mask(r=1.0, ang=90, anchor=CENTER, spin=0, orient=UP)
+module rounding_angled_corner_mask(r, ang=90, d, anchor=CENTER, spin=0, orient=UP)
 {
+    r = get_radius(r=r, d=d, dflt=1);
     dx = r / tan(ang/2);
     dx2 = dx / cos(ang/2) + 1;
     fn = quantup(segs(r), 4);
@@ -683,13 +702,14 @@ module rounding_angled_corner_mask(r=1.0, ang=90, anchor=CENTER, spin=0, orient=
 
 // Module: rounding_corner_mask()
 // Usage:
-//   rounding_corner_mask(r, [anchor]);
+//   rounding_corner_mask(r|d, [anchor]);
 // Description:
 //   Creates a shape that you can use to round 90 degree corners.
 //   Difference it from the object to be rounded.  The center of the mask
 //   object should align exactly with the corner to be rounded.
 // Arguments:
 //   r = Radius of corner rounding.
+//   d = Diameter of corner rounding.
 //   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`
@@ -703,8 +723,9 @@ module rounding_angled_corner_mask(r=1.0, ang=90, anchor=CENTER, spin=0, orient=
 //     translate([15, 25, 0]) rounding_mask_z(l=81, r=15);
 //     translate([15, 25, 40]) #rounding_corner_mask(r=15);
 //   }
-module rounding_corner_mask(r=1.0, anchor=CENTER, spin=0, orient=UP)
+module rounding_corner_mask(r, d, anchor=CENTER, spin=0, orient=UP)
 {
+    r = get_radius(r=r, d=d, dflt=1);
     attachable(anchor,spin,orient, size=[2,2,2]*r) {
         difference() {
             cube(size=r*2, center=true);
@@ -719,7 +740,7 @@ module rounding_corner_mask(r=1.0, anchor=CENTER, spin=0, orient=UP)
 
 // Module: rounding_cylinder_mask()
 // Usage:
-//   rounding_cylinder_mask(r, rounding);
+//   rounding_cylinder_mask(r|d, rounding);
 // Description:
 //   Create a mask that can be used to round the end of a cylinder.
 //   Difference it from the cylinder to be rounded.  The center of the
@@ -727,6 +748,7 @@ module rounding_corner_mask(r=1.0, anchor=CENTER, spin=0, orient=UP)
 //   cylinder to be rounded.
 // Arguments:
 //   r = Radius of cylinder. (Default: 1.0)
+//   d = Diameter of cylinder. (Default: 1.0)
 //   rounding = Radius of the edge rounding. (Default: 0.25)
 // Example:
 //   difference() {
@@ -738,8 +760,9 @@ module rounding_corner_mask(r=1.0, anchor=CENTER, spin=0, orient=UP)
 //     cylinder(r=50, h=50, center=false);
 //     up(50) rounding_cylinder_mask(r=50, rounding=10);
 //   }
-module rounding_cylinder_mask(r=1.0, rounding=0.25)
+module rounding_cylinder_mask(r, rounding=0.25, d)
 {
+    r = get_radius(r=r, d=d, dflt=1);
     cylinder_mask(l=rounding*3, r=r, rounding2=rounding, excess=rounding, ends_only=true, anchor=TOP);
 }
 
@@ -775,7 +798,7 @@ 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, excess=0.1, anchor=CENTER, spin=0, orient=UP)
+module rounding_hole_mask(r, rounding=0.25, excess=0.1, d, 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) {
@@ -813,7 +836,7 @@ module rounding_hole_mask(r=undef, d=undef, rounding=0.25, excess=0.1, anchor=CE
 //          corner_profile(BOT,r=10)
 //              mask2d_teardrop(r=10, angle=40);
 //   }
-module teardrop_corner_mask(r, d, angle, excess=0.1, anchor=CENTER, spin=0, orient=UP) {
+module teardrop_corner_mask(r, angle, excess=0.1, d, anchor=CENTER, spin=0, orient=UP) {
     assert(is_num(angle));
     assert(is_num(excess));
     assert(angle>0 && angle<90);
diff --git a/math.scad b/math.scad
index 90182f1..599a0c4 100644
--- a/math.scad
+++ b/math.scad
@@ -680,7 +680,7 @@ function convolve(p,q) =
 //   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
 //   transpose the returned value.  
-function linear_solve(A,b) =
+function linear_solve(A,b,pivot=true) =
     assert(is_matrix(A), "Input should be a matrix.")
     let(
         m = len(A),
@@ -688,19 +688,17 @@ function linear_solve(A,b) =
     )
     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),
+        qr = m<n? qr_factor(transpose(A),pivot) : qr_factor(A,pivot),
         maxdim = max(n,m),
         mindim = min(n,m),
         Q = submatrix(qr[0],[0:maxdim-1], [0:mindim-1]),
         R = submatrix(qr[1],[0:mindim-1], [0:mindim-1]),
+        P = qr[2],
         zeros = [for(i=[0:mindim-1]) if (approx(R[i][i],0)) i]
     )
     zeros != [] ? [] :
-    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);
+    m<n ? Q*back_substitute(R,transpose(P)*b,transpose=true) // Too messy to avoid input checks here
+        : P*_back_substitute(R, transpose(Q)*b);             // Calling internal version skips input checks
 
 // Function: matrix_inverse()
 // Usage:
@@ -714,20 +712,40 @@ function matrix_inverse(A) =
     assert(is_matrix(A,square=true),"Input to matrix_inverse() must be a square matrix")
     linear_solve(A,ident(len(A)));
 
+// Function: null_space()
+// Usage:
+//   null_space(A)
+// Description:
+//   Returns an orthonormal basis for the null space of A, namely the vectors {x} such that Ax=0.  If the null space
+//   is just the origin then returns an empty list. 
+function null_space(A,eps=1e-12) =
+    assert(is_matrix(A))
+    let(
+      Q_R=qr_factor(transpose(A),pivot=true),
+      R=Q_R[1],
+      zrow = [for(i=idx(R)) if (is_zero(R[i],eps)) i]
+    )
+    len(zrow)==0
+      ? []
+      : transpose(subindex(Q_R[0],zrow));
+
 
 // Function: qr_factor()
-// Usage: qr = qr_factor(A)
+// Usage: qr = qr_factor(A,[pivot])
 // Description:
-//   Calculates the QR factorization of the input matrix A and returns it as the list [Q,R].  This factorization can be
-//   used to solve linear systems of equations.  
-function qr_factor(A) =
+//   Calculates the QR factorization of the input matrix A and returns it as the list [Q,R,P].  This factorization can be
+//   used to solve linear systems of equations.  The factorization is A = Q*R*transpose(P).  If pivot is false (the default)
+//   then P is the identity matrix and A = Q*R.  If pivot is true then column pivoting results in an R matrix where the diagonal
+//   is non-decreasing.  The use of pivoting is supposed to increase accuracy for poorly conditioned problems, and is necessary
+//   for rank estimation or computation of the null space, but it may be slower.  
+function qr_factor(A, pivot=false) =
     assert(is_matrix(A), "Input must be a matrix." )
     let(
       m = len(A),
       n = len(A[0])
     )
     let(
-        qr = _qr_factor(A, Q=ident(m), column=0, m = m, n=n),
+        qr =_qr_factor(A, Q=ident(m),P=ident(n), pivot=pivot, column=0, m = m, n=n),
         Rzero = 
           let( R = qr[1] )
           [ for(i=[0:m-1]) [
@@ -735,25 +753,31 @@ function qr_factor(A) =
               for(j=[0:n-1]) i>j ? 0 : ri[j]
             ]
           ]
-    ) [qr[0],Rzero];
+    ) [qr[0],Rzero,qr[2]];
 
-function _qr_factor(A,Q, column, m, n) =
-    column >= min(m-1,n) ? [Q,A] :
+function _qr_factor(A,Q,P, pivot, column, m, n) =
+    column >= min(m-1,n) ? [Q,A,P] :
     let(
+        swap = !pivot ? 1
+             : _swap_matrix(n,column,column+max_index([for(i=[column:n-1]) sum_of_squares([for(j=[column:m-1]) A[j][i]])])),
+        A = pivot ? A*swap : A,
         x = [for(i=[column:1:m-1]) A[i][column]],
         alpha = (x[0]<=0 ? 1 : -1) * norm(x),
         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);
+    _qr_factor(Qf*A, Q*Qf, P*swap, pivot, column+1, m, n);
+
+// Produces an n x n matrix that swaps column i and j (when multiplied on the right)
+function _swap_matrix(n,i,j) =
+  assert(i<n && j<n && i>=0 && j>=0, "Swap indices out of bounds")
+  [for(y=[0:n-1]) [for (x=[0:n-1])
+     x==i ? (y==j ? 1 : 0)
+   : x==j ? (y==i ? 1 : 0)
+   : x==y ? 1 : 0]];
+
 
 
 // Function: back_substitute()
@@ -862,6 +886,17 @@ function is_matrix(A,m,n,square=false) =
     && ( !square || len(A)==len(A[0]));
 
 
+// Function: norm_fro()
+// Usage:
+//    norm_fro(A)
+// Description:
+//    Computes frobenius norm of input matrix.  The frobenius norm is the square root of the sum of the
+//    squares of all of the entries of the matrix.  On vectors it is the same as the usual 2-norm.
+//    This is an easily computed norm that is convenient for comparing two matrices.  
+function norm_fro(A) =
+    sqrt(sum([for(entry=A) sum_of_squares(entry)]));
+
+
 // Section: Comparisons and Logic
 
 // Function: is_zero()
@@ -1292,7 +1327,7 @@ function deriv3(data, h=1, closed=false) =
 // Description:
 //   Multiplies two complex numbers represented by 2D vectors.  
 function C_times(z1,z2) = 
-    assert( is_vector(z1+z2,2), "Complex numbers should be represented by 2D vectors." )
+    assert( is_matrix([z1,z2],2,2), "Complex numbers should be represented by 2D vectors" )
     [ z1.x*z2.x - z1.y*z2.y, z1.x*z2.y + z1.y*z2.x ];
 
 // Function: C_div()
@@ -1309,6 +1344,41 @@ function C_div(z1,z2) =
 
 // Section: Polynomials
 
+// Function: quadratic_roots()
+// Usage:
+//    roots = quadratic_roots(a,b,c,[real])
+// Description:
+//    Computes roots of the quadratic equation a*x^2+b*x+c==0, where the
+//    coefficients are real numbers.  If real is true then returns only the
+//    real roots.  Otherwise returns a pair of complex values.  This method
+//    may be more reliable than the general root finder at distinguishing
+//    real roots from complex roots.  
+
+// https://people.csail.mit.edu/bkph/articles/Quadratics.pdf
+
+function quadratic_roots(a,b,c,real=false) =
+  real ? [for(root = quadratic_roots(a,b,c,real=false)) if (root.y==0) root.x]
+  :
+  is_undef(b) && is_undef(c) && is_vector(a,3) ? quadratic_roots(a[0],a[1],a[2]) :
+  assert(is_num(a) && is_num(b) && is_num(c))
+  assert(a!=0 || b!=0 || c!=0, "Quadratic must have a nonzero coefficient")
+  a==0 && b==0 ? [] :     // No solutions
+  a==0 ? [[-c/b,0]] : 
+  let(
+      descrim = b*b-4*a*c,
+      sqrt_des = sqrt(abs(descrim))
+  )
+  descrim < 0 ?             // Complex case
+     [[-b, sqrt_des],
+      [-b, -sqrt_des]]/2/a :
+  b<0 ?                     // b positive
+     [[2*c/(-b+sqrt_des),0],
+      [(-b+sqrt_des)/a/2,0]]
+      :                     // b negative
+     [[(-b-sqrt_des)/2/a, 0],
+      [2*c/(-b-sqrt_des),0]];
+
+
 // Function: polynomial() 
 // Usage:
 //   polynomial(p, z)
diff --git a/mutators.scad b/mutators.scad
index b801e48..4db899b 100644
--- a/mutators.scad
+++ b/mutators.scad
@@ -321,7 +321,7 @@ module chain_hull()
 // Usage:
 //   cylindrical_extrude(size, ir|id, or|od, [convexity]) ...
 // Description:
-//   Cylindrically extrudes all 2D children, curved around a cylidrical shape.
+//   Extrudes all 2D children outwards, curved around a cylindrical shape.
 // Arguments:
 //   or = The outer radius to extrude to.
 //   od = The outer diameter to extrude to.
diff --git a/paths.scad b/paths.scad
index 93d82fa..60bc34e 100644
--- a/paths.scad
+++ b/paths.scad
@@ -418,7 +418,7 @@ function path_torsion(path, closed=false) =
 //   scale = [X,Y] scaling factors for each axis.  Default: `[1,1]`
 // Example(3D):
 //   trace_polyline(path3d_spiral(turns=2.5, h=100, n=24, r=50), N=1, showpts=true);
-function path3d_spiral(turns=3, h=100, n=12, r=undef, d=undef, cp=[0,0], scale=[1,1]) = let(
+function path3d_spiral(turns=3, h=100, n=12, r, d, cp=[0,0], scale=[1,1]) = let(
         rr=get_radius(r=r, d=d, dflt=100),
         cnt=floor(turns*n),
         dz=h/cnt
@@ -774,15 +774,19 @@ function assemble_path_fragments(fragments, eps=EPSILON, _finished=[]) =
 
 
 // Module: modulated_circle()
+// Usage:
+//   modulated_circle(r|d, sines);
 // Description:
 //   Creates a 2D polygon circle, modulated by one or more superimposed sine waves.
 // Arguments:
-//   r = radius of the base circle.
+//   r = Radius of the base circle. Default: 40
+//   d = Diameter of the base circle.
 //   sines = array of [amplitude, frequency] pairs, where the frequency is the number of times the cycle repeats around the circle.
 // Example(2D):
 //   modulated_circle(r=40, sines=[[3, 11], [1, 31]], $fn=6);
-module modulated_circle(r=40, sines=[10])
+module modulated_circle(r, sines=[10], d)
 {
+    r = get_radius(r=r, d=d, dflt=40);
     freqs = len(sines)>0? [for (i=sines) i[1]] : [5];
     points = [
         for (a = [0 : (360/segs(r)/max(freqs)) : 360])
@@ -829,7 +833,8 @@ module extrude_from_to(pt1, pt2, convexity=undef, twist=undef, scale=undef, slic
 // Arguments:
 //   polyline = Array of points of a polyline path, to be extruded.
 //   h = height of the spiral to extrude along.
-//   r = radius of the spiral to extrude along.
+//   r = Radius of the spiral to extrude along. Default: 50
+//   d = Diameter of the spiral to extrude along.
 //   twist = number of degrees of rotation to spiral up along height.
 //   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`
@@ -838,7 +843,8 @@ module extrude_from_to(pt1, pt2, convexity=undef, twist=undef, scale=undef, slic
 // Example:
 //   poly = [[-10,0], [-3,-5], [3,-5], [10,0], [0,-30]];
 //   spiral_sweep(poly, h=200, r=50, twist=1080, $fn=36);
-module spiral_sweep(polyline, h, r, twist=360, center, anchor, spin=0, orient=UP) {
+module spiral_sweep(polyline, h, r, twist=360, center, d, anchor, spin=0, orient=UP) {
+    r = get_radius(r=r, d=d, dflt=50);
     polyline = path3d(polyline);
     pline_count = len(polyline);
     steps = ceil(segs(r)*(twist/360));
diff --git a/polyhedra.scad b/polyhedra.scad
index e845acd..13cbe54 100644
--- a/polyhedra.scad
+++ b/polyhedra.scad
@@ -652,7 +652,7 @@ function regular_polyhedron_info(
     let(
         entry = (
             name == "trapezohedron"? (
-                trapezohedron(faces=faces, side=side, longside=longside, h=h, r=r)
+                _trapezohedron(faces=faces, side=side, longside=longside, h=h, r=r)
             ) : (
                 _polyhedra_[!is_undef(index)?
                     indexlist[index] :
@@ -671,7 +671,7 @@ function regular_polyhedron_info(
             ) / entry[edgelen]
         ),
         face_triangles = hull(entry[vertices]),
-        faces_normals_vertices = stellate_faces(
+        faces_normals_vertices = _stellate_faces(
             entry[edgelen], stellate, entry[vertices],
             entry[facevertices]==[3]?
                 [face_triangles, [for(face=face_triangles) _facenormal(entry[vertices],face)]] :
@@ -713,11 +713,7 @@ function regular_polyhedron_info(
     assert(false, str("Unknown info type '",info,"' requested"));
 
 
-
-/// hull solution fails due to roundoff
-/// either cross product or just rotate to
-///
-function stellate_faces(scalefactor,stellate,vertices,faces_normals) =
+function _stellate_faces(scalefactor,stellate,vertices,faces_normals) =
     (stellate == false || stellate == 0)? concat(faces_normals,[vertices]) :
     let(
         faces = [for(face=faces_normals[0]) select(face,hull(select(vertices,face)))],
@@ -730,9 +726,10 @@ function stellate_faces(scalefactor,stellate,vertices,faces_normals) =
     ) [newfaces, normals, allpts];
 
 
-function trapezohedron(faces, r, side, longside, h) =
-    assert(faces%2==0, "Number of faces must be even")
+function _trapezohedron(faces, r, side, longside, h, d) =
+    assert(faces%2==0, "Must set 'faces' to an even number for trapezohedron")
     let(
+        r = get_radius(r=r, d=d, dflt=1),
         N = faces/2,
         parmcount = num_defined([r,side,longside,h])
     )
diff --git a/rounding.scad b/rounding.scad
index ae2aed0..c57a76f 100644
--- a/rounding.scad
+++ b/rounding.scad
@@ -488,6 +488,7 @@ function smooth_path(path, tangents, size, relsize, splinesteps=10, uniform=fals
 //   - smooth: os_smooth(cut|joint).  Define continuous curvature rounding, with `cut` and `joint` as for round_corners.
 //   - teardrop: os_teardrop(r|cut).  Rounding using a 1/8 circle that then changes to a 45 degree chamfer.  The chamfer is at the end, and enables the object to be 3d printed without support.  The radius gives the radius of the circular part.
 //   - chamfer: os_chamfer([height], [width], [cut], [angle]).  Chamfer the edge at desired angle or with desired height and width.  You can specify height and width together and the angle will be ignored, or specify just one of height and width and the angle is used to determine the shape.  Alternatively, specify "cut" along with angle to specify the cut back distance of the chamfer.
+//   - mask: os_mask(mask, [out]).  Create a profile from one of the [2d masking shapes](shapes2d.scad#5-2d-masking-shapes).  The `out` parameter specifies that the mask should flare outward (like crown molding or baseboard).  This is set false by default.  
 //   .
 //   The general settings that you can use with all of the helper functions are mostly used to control how offset_sweep() calls the offset() function.
 //   - extra: Add an extra vertical step of the specified height, to be used for intersections or differences.  This extra step will extend the resulting object beyond the height you specify.  Default: 0
@@ -654,6 +655,15 @@ 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));
 //     }
+// Example: Using os_mask to create ogee profiles:
+//   ogee = mask2d_ogee([
+//       "xstep",1,  "ystep",1,  // Starting shoulder.
+//       "fillet",5, "round",5,  // S-curve.
+//       "ystep",1,              // Ending shoulder.
+//   ]);
+//   star = star(5, r=220, ir=130);
+//   rounded_star = round_corners(star, cut=flatten(repeat([5,0],5)), $fn=24);
+//   offset_sweep(rounded_star, height=100, top=os_mask(ogee), bottom=os_mask(ogee,out=true));
 
 
 // This function does the actual work of repeatedly calling offset() and concatenating the resulting face and vertex lists to produce
@@ -880,6 +890,18 @@ function os_profile(points, extra,check_valid, quality, offset_maxstep, offset)
         ]);
 
 
+function os_mask(mask, out=false, extra,check_valid, quality, offset_maxstep, offset) =
+  let(
+      origin_index = [for(i=idx(mask)) if (mask[i].x<0 && mask[i].y<0) i],
+      xfactor = out ? -1 : 1
+  )
+  assert(len(origin_index)==1,"Cannot find origin in the mask")
+  let(
+      points = ([for(pt=polygon_shift(mask,origin_index[0])) [xfactor*max(pt.x,0),-max(pt.y,0)]])
+  )
+  os_profile(deduplicate(move(-points[1],p=select(points,1,-1))), extra,check_valid,quality,offset_maxstep,offset);
+
+
 // Module: convex_offset_extrude()
 //
 // Description:
@@ -1994,4 +2016,4 @@ module bent_cutout_mask(r, thickness, path, convexity=10)
 }
 
 
-// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap
\ No newline at end of file
+// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap
diff --git a/scripts/gencheat.sh b/scripts/gencheat.sh
new file mode 100755
index 0000000..e07a909
--- /dev/null
+++ b/scripts/gencheat.sh
@@ -0,0 +1,98 @@
+#!/bin/bash
+
+
+function ucase
+{
+    echo "$1" | tr '[:lower:]' '[:upper:]'
+}
+
+function lcase
+{
+    echo "$1" | tr '[:upper:]' '[:lower:]'
+}
+
+function columnize
+{
+    cols=$2
+    TMPFILE=$(mktemp -t $(basename $0).XXXXXX) || exit 1
+    cat >>$TMPFILE
+    if [[ $(wc -l $TMPFILE | awk '{print $1}') -gt 1 ]] ; then
+        totcnt=$(wc -l $TMPFILE | awk '{print $1}')
+        maxrows=$((($totcnt+$cols-1)/$cols))
+        maxcols=$cols
+        if [[ $maxcols -gt $totcnt ]] ; then
+            maxcols=$totcnt
+        fi
+        cnt=0
+        hdrln1="| $(ucase $1)  "
+        hdrln2='|:-----'
+        n=1
+        while [[ $n -lt $maxcols ]] ; do
+            hdrln1+=' | &nbsp;'
+            hdrln2+=' |:------'
+            n=$(($n+1))
+        done
+        hdrln1+=' |'
+        hdrln2+=' |'
+        n=0
+        while [[ $n -lt $maxrows ]] ; do
+            lines[$n]=""
+            n=$(($n+1))
+        done
+        col=0
+        while IFS= read -r line; do
+            if [[ $col != 0 ]] ; then
+                lines[$cnt]+=" | "
+            fi
+            lines[$cnt]+="$line"
+            cnt=$(($cnt+1))
+            if [[ $cnt = $maxrows ]] ; then
+                cnt=0
+                col=$(($col+1))
+            fi
+        done <$TMPFILE
+        rm -f $TMPFILE
+
+        echo
+        echo $hdrln1
+        echo $hdrln2
+        n=0
+        while [[ $n -lt $maxrows ]] ; do
+            echo "| ${lines[$n]} |"
+            n=$(($n+1))
+        done
+    fi
+}
+
+function mkconstindex
+{
+    sed 's/([^)]*)//g' | sed 's/[^a-zA-Z0-9_.:$]//g' | awk -F ':' '{printf "[%s](%s#%s)\n", $3, $1, $3}'
+}
+
+function mkconstindex2
+{
+    sed 's/ *=.*$//' | sed 's/[^a-zA-Z0-9_.:$]//g' | awk -F ':' '{printf "[%s](%s#%s)\n", $2, $1, $2}'
+}
+
+function mkotherindex
+{
+    sed 's/([^)]*)//g' | sed 's/[^a-zA-Z0-9_.:$]//g' | awk -F ':' '{printf "[%s()](%s#%s)\n", $3, $1, $3}'
+}
+
+CHEAT_FILES=$(grep '^include' std.scad | sed 's/^.*<\([a-zA-Z0-9.]*\)>/\1/' | grep -v 'version.scad' | grep -v 'primitives.scad')
+
+(
+    echo '## Belfry OpenScad Library Cheat Sheet'
+    echo
+    echo '( [Alphabetic Index](Index) )'
+    echo
+    (
+        grep -H '// Constant: ' $CHEAT_FILES | mkconstindex
+        grep -H '^[A-Z$][A-Z0-9_]* *=.*//' $CHEAT_FILES | mkconstindex2
+    ) | sort -u | columnize 'Constants' 6
+    for f in $CHEAT_FILES ; do
+        egrep -H '// Function: |// Function&Module: |// Module: ' $f | mkotherindex | columnize "[$f]($f)" 4
+        echo
+    done
+) > BOSL2.wiki/CheatSheet.md
+
diff --git a/scripts/genindex.sh b/scripts/genindex.sh
index 9d3e8cd..51a9075 100755
--- a/scripts/genindex.sh
+++ b/scripts/genindex.sh
@@ -6,16 +6,18 @@ function ucase
     echo "$1" | tr '[:lower:]' '[:upper:]'
 }
 
+
 function lcase
 {
     echo "$1" | tr '[:upper:]' '[:lower:]'
 }
 
-function mkindex
+
+function alphaindex
 {
-    TMPFILE=$(mktemp -t $(basename $0).XXXXXX) || exit 1
-    sed 's/([^)]*)//g' | sed 's/[^a-zA-Z0-9_.:$]//g' | awk -F ':' '{printf "- [%s](%s#%s)\n", $3, $1, $3}' | sort -d -f -u >> $TMPFILE
     alpha="A B C D E F G H I J K L M N O P Q R S T U V W X Y Z"
+    TMPFILE=$(mktemp -t $(basename $0).XXXXXX) || exit 1
+    sort -d -f >> $TMPFILE
     for a in $alpha; do
         echo -n "[$a](#$(lcase "$a")) "
     done
@@ -33,13 +35,31 @@ function mkindex
 }
 
 
+function constlist
+{
+    sed 's/([^)]*)//g' | sed 's/[^a-zA-Z0-9_.:$]//g' | awk -F ':' '{printf "- [%s](%s#%s) (in %s)\n", $3, $1, $3, $1}'
+}
+ 
+function constlist2
+{
+    sed 's/ *=.*$//' | sed 's/[^a-zA-Z0-9_.:$]//g' | awk -F ':' '{printf "- [%s](%s#%s) (in %s)\n", $2, $1, $2, $1}'
+}
+ 
+
+function funclist
+{
+    sed 's/([^)]*)//g' | sed 's/[^a-zA-Z0-9_.:$]//g' | awk -F ':' '{printf "- [%s()](%s#%s) (in %s)\n", $3, $1, $3, $1}'
+}
+
+
 (
     echo "## Belfry OpenScad Library Index"
     (
-        grep 'Constant: ' *.scad
-        grep 'Function: ' *.scad
-        grep 'Function&Module: ' *.scad
-        grep 'Module: ' *.scad
-    ) | mkindex
+    	(
+	    grep 'Constant: ' *.scad | constlist
+	    grep '^[A-Z]* *=.*//' *.scad | constlist2
+	) | sort -u
+        egrep 'Function: |Function&Module: |Module: ' *.scad | sort -u | funclist
+    ) | sort | alphaindex
 ) > BOSL2.wiki/Index.md
 
diff --git a/shapes.scad b/shapes.scad
index d612b0a..b010121 100644
--- a/shapes.scad
+++ b/shapes.scad
@@ -1498,13 +1498,14 @@ module pie_slice(
 //   Center this part along the concave edge to be chamfered and union it in.
 //
 // Usage:
-//   interior_fillet(l, r, [ang], [overlap]);
+//   interior_fillet(l, r|d, [ang], [overlap]);
 //
 // Arguments:
-//   l = length of edge to fillet.
-//   r = radius of fillet.
-//   ang = angle between faces to fillet.
-//   overlap = overlap size for unioning with faces.
+//   l = Length of edge to fillet.
+//   r = Radius of fillet.
+//   d = Diameter of fillet.
+//   ang = Angle between faces to fillet.
+//   overlap = Overlap size for unioning with faces.
 //   anchor = Translate so anchor point is at origin (0,0,0).  See [anchor](attachments.scad#anchor).  Default: `FRONT+LEFT`
 //   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`
@@ -1526,7 +1527,8 @@ module pie_slice(
 //     position(BOT+FRONT)
 //       interior_fillet(l=50, r=10, spin=180, orient=RIGHT);
 //   }
-module interior_fillet(l=1.0, r=1.0, ang=90, overlap=0.01, anchor=FRONT+LEFT, spin=0, orient=UP) {
+module interior_fillet(l=1.0, r, ang=90, overlap=0.01, d, anchor=FRONT+LEFT, spin=0, orient=UP) {
+    r = get_radius(r=r, d=d, dflt=1);
     dy = r/tan(ang/2);
     steps = ceil(segs(r)*ang/360);
     step = ang/steps;
diff --git a/std.scad b/std.scad
index 1c2e45e..685b6ab 100644
--- a/std.scad
+++ b/std.scad
@@ -12,21 +12,6 @@ assert(version_num()>=20190301, "BOSL2 requires OpenSCAD version 2019.03.01 or l
 include <version.scad>
 
 include <constants.scad>
-include <edges.scad>
-include <common.scad>
-include <arrays.scad>
-include <strings.scad>
-include <vnf.scad>
-include <debug.scad>
-
-include <math.scad>
-include <vectors.scad>
-include <quaternions.scad>
-include <affine.scad>
-include <coords.scad>
-include <geometry.scad>
-include <regions.scad>
-
 include <transforms.scad>
 include <distributors.scad>
 include <mutators.scad>
@@ -36,6 +21,19 @@ include <shapes.scad>
 include <shapes2d.scad>
 include <masks.scad>
 include <paths.scad>
+include <edges.scad>
+include <arrays.scad>
+include <math.scad>
+include <vectors.scad>
+include <quaternions.scad>
+include <affine.scad>
+include <coords.scad>
+include <geometry.scad>
+include <regions.scad>
+include <strings.scad>
+include <vnf.scad>
+include <common.scad>
+include <debug.scad>
 
 
 // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap
diff --git a/strings.scad b/strings.scad
index cc0713b..5e34fd9 100644
--- a/strings.scad
+++ b/strings.scad
@@ -701,7 +701,9 @@ function str_format(fmt, vals, use_nbsp=false) =
 //   echofmt("{:-10s}{:.3f}", ["plecostamus",27.43982]);  // ECHO: "plecostamus27.440"
 //   echofmt("{:-10.9s}{:.3f}", ["plecostamus",27.43982]);  // ECHO: "plecostam 27.440"
 function echofmt(fmt, vals, use_nbsp=false) = echo(str_format(fmt,vals,use_nbsp));
-module echofmt(fmt, vals, use_nbsp=false) echo(str_format(fmt,vals,use_nbsp));
-
+module echofmt(fmt, vals, use_nbsp=false) {
+   no_children($children);
+   echo(str_format(fmt,vals,use_nbsp));
+}
 
 // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap
diff --git a/structs.scad b/structs.scad
index e7effb4..1e9f321 100644
--- a/structs.scad
+++ b/structs.scad
@@ -101,6 +101,7 @@ function struct_echo(struct,name="") =
     undef;
 
 module struct_echo(struct,name="") {
+    no_children($children);
     dummy = struct_echo(struct,name);
 }
 
diff --git a/tests/test_arrays.scad b/tests/test_arrays.scad
index 8df23a8..1dfe421 100644
--- a/tests/test_arrays.scad
+++ b/tests/test_arrays.scad
@@ -3,6 +3,16 @@ include <../std.scad>
 
 // Section: List Query Operations
 
+module test_is_homogeneous(){
+    assert(is_homogeneous([[1,["a"]], [2,["b"]]])==true);
+    assert(is_homogeneous([[1,["a"]], [2,[true]]])==false);
+    assert(is_homogeneous([[1,["a"]], [2,[true]]],1)==true);
+    assert(is_homogeneous([[1,["a"]], [2,[true]]],2)==false);
+    assert(is_homogeneous([[1,["a"]], [true,["b"]]])==true);
+}
+test_is_homogeneous();
+
+
 module test_select() {
     l = [3,4,5,6,7,8,9];
     assert(select(l, 5, 6) == [8,9]);
@@ -46,7 +56,6 @@ test_in_list();
 module test_min_index() {
     assert(min_index([5,3,9,6,2,7,8,2,1])==8);
     assert(min_index([5,3,9,6,2,7,8,2,7],all=true)==[4,7]);
-//    assert(min_index([],all=true)==[]);
 }
 test_min_index();
 
@@ -54,7 +63,6 @@ test_min_index();
 module test_max_index() {
     assert(max_index([5,3,9,6,2,7,8,9,1])==2);
     assert(max_index([5,3,9,6,2,7,8,9,7],all=true)==[2,7]);
-//    assert(max_index([],all=true)==[]);
 }
 test_max_index();
 
@@ -74,6 +82,7 @@ module test_list_decreasing() {
 }
 test_list_decreasing();
 
+
 // Section: Basic List Generation
 
 module test_repeat() {
@@ -155,7 +164,6 @@ module test_list_remove() {
 }
 test_list_remove();
 
-
 module test_list_remove_values() {
     animals = ["bat", "cat", "rat", "dog", "bat", "rat"];
     assert(list_remove_values(animals, "rat") == ["bat","cat","dog","bat","rat"]);
@@ -258,15 +266,24 @@ test_shuffle();
 
 module test_sort() {
     assert(sort([7,3,9,4,3,1,8]) == [1,3,3,4,7,8,9]);
-    assert(sort(["cat", "oat", "sat", "bat", "vat", "rat", "pat", "mat", "fat", "hat", "eat"]) == ["bat", "cat", "eat", "fat", "hat", "mat", "oat", "pat", "rat", "sat", "vat"]);
+    assert(sort([[4,0],[7],[3,9],20,[4],[3,1],[8]]) == [20,[3,1],[3,9],[4],[4,0],[7],[8]]);
+    assert(sort([[4,0],[7],[3,9],20,[4],[3,1],[8]],idx=1) == [[7],20,[4],[8],[4,0],[3,1],[3,9]]);
+    assert(sort([[8,6],[3,1],[9,2],[4,3],[3,4],[1,5],[8,0]]) == [[1,5],[3,1],[3,4],[4,3],[8,0],[8,6],[9,2]]);
+    assert(sort([[8,0],[3,1],[9,2],[4,3],[3,4],[1,5],[8,6]],idx=1) == [[8,0],[3,1],[9,2],[4,3],[3,4],[1,5],[8,6]]);
+    assert(sort(["cat", "oat", "sat", "bat", "vat", "rat", "pat", "mat", "fat", "hat", "eat"]) 
+           == ["bat", "cat", "eat", "fat", "hat", "mat", "oat", "pat", "rat", "sat", "vat"]);
     assert(sort(enumerate([[2,3,4],[1,2,3],[2,4,3]]),idx=1)==[[1,[1,2,3]], [0,[2,3,4]], [2,[2,4,3]]]);
+    assert(sort([0,"1",[1,0],2,"a",[1]])== [0,2,"1","a",[1],[1,0]]);
+    assert(sort([["oat",0], ["cat",1], ["bat",3], ["bat",2], ["fat",3]])==  [["bat",2],["bat",3],["cat",1],["fat",3],["oat",0]]);
 }
 test_sort();
 
 
 module test_sortidx() {
-    lst1 = ["d","b","e","c"];
+    lst1 = ["da","bax","eaw","cav"];
     assert(sortidx(lst1) == [1,3,0,2]);
+    lst5 = [3,5,1,7];
+    assert(sortidx(lst5) == [2,0,1,3]);
     lst2 = [
         ["foo", 88, [0,0,1], false],
         ["bar", 90, [0,1,0], true],
@@ -276,11 +293,22 @@ module test_sortidx() {
     assert(sortidx(lst2, idx=1) == [3,0,2,1]);
     assert(sortidx(lst2, idx=0) == [1,2,0,3]);
     assert(sortidx(lst2, idx=[1,3]) == [3,0,2,1]);
-    lst3 = [[-4, 0, 0], [0, 0, -4], [0, -4, 0], [-4, 0, 0], [0, -4, 0], [0, 0, 4], [0, 0, -4], [0, 4, 0], [4, 0, 0], [0, 0, 4], [0, 4, 0], [4, 0, 0]];
+    lst3 = [[-4,0,0],[0,0,-4],[0,-4,0],[-4,0,0],[0,-4,0],[0,0,4],
+            [0,0,-4],[0,4,0],[4,0,0],[0,0,4],[0,4,0],[4,0,0]];
     assert(sortidx(lst3)==[0,3,2,4,1,6,5,9,7,10,8,11]);
+    assert(sortidx([[4,0],[7],[3,9],20,[4],[3,1],[8]]) == [3,5,2,4,0,1,6]);
+    assert(sortidx([[4,0],[7],[3,9],20,[4],[3,1],[8]],idx=1) ==  [1,3,4,6,0,5,2]);
+    lst4=[0,"1",[1,0],2,"a",[1]];
+    assert(sortidx(lst4)== [0,3,1,4,5,2]);
+    assert(sortidx(["cat","oat","sat","bat","vat","rat","pat","mat","fat","hat","eat"]) 
+             == [3,0,10,8,9,7,1,6,5,2,4]);
+    assert(sortidx([["oat",0], ["cat",1], ["bat",3], ["bat",2], ["fat",3]])==  [3,2,1,4,0]);
+    assert(sortidx(["Belfry", "OpenScad", "Library", "Documentation"])==[0,3,2,1]);
+    assert(sortidx(["x",1,[],0,"abc",true])==[5,3,1,4,0,2]);
 }
 test_sortidx();
 
+
 module test_unique() {
     assert(unique([]) == []);
     assert(unique([8]) == [8]);
@@ -336,10 +364,8 @@ module test_set_intersection() {
 test_set_intersection();
 
 
-
 // Arrays
 
-
 module test_add_scalar() {
     assert(add_scalar([1,2,3],3) == [4,5,6]);
     assert(add_scalar([[1,2,3],[3,4,5]],3) == [[4,5,6],[6,7,8]]);
@@ -444,6 +470,40 @@ module test_zip() {
 }
 test_zip();
 
+module test_block_matrix() {
+    A = [[1,2],[3,4]];
+    B = ident(2);
+    assert_equal(block_matrix([[A,B],[B,A],[A,B]]), [[1,2,1,0],[3,4,0,1],[1,0,1,2],[0,1,3,4],[1,2,1,0],[3,4,0,1]]);
+    assert_equal(block_matrix([[A,B],ident(4)]), [[1,2,1,0],[3,4,0,1],[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]]);
+    text = [["a","b"],["c","d"]];
+    assert_equal(block_matrix([[text,B]]), [["a","b",1,0],["c","d",0,1]]);
+}
+test_block_matrix();
+
+
+module test_diagonal_matrix() {
+    assert_equal(diagonal_matrix([1,2,3]), [[1,0,0],[0,2,0],[0,0,3]]);
+    assert_equal(diagonal_matrix([1,"c",2]), [[1,0,0],[0,"c",0],[0,0,2]]);
+    assert_equal(diagonal_matrix([1,"c",2],"X"), [[1,"X","X"],["X","c","X"],["X","X",2]]);
+    assert_equal(diagonal_matrix([[1,1],[2,2],[3,3]], [0,0]), [[ [1,1],[0,0],[0,0]], [[0,0],[2,2],[0,0]], [[0,0],[0,0],[3,3]]]);
+}
+test_diagonal_matrix();
+
+module test_submatrix_set() {
+    test = [[1,2,3,4,5],[6,7,8,9,10],[11,12,13,14,15], [16,17,18,19,20]];
+    ragged = [[1,2,3,4,5],[6,7,8,9,10],[11,12], [16,17]];
+    assert_equal(submatrix_set(test,[[9,8],[7,6]]), [[9,8,3,4,5],[7,6,8,9,10],[11,12,13,14,15], [16,17,18,19,20]]);
+    assert_equal(submatrix_set(test,[[9,7],[8,6]],1),[[1,2,3,4,5],[9,7,8,9,10],[8,6,13,14,15], [16,17,18,19,20]]);
+    assert_equal(submatrix_set(test,[[9,8],[7,6]],n=1), [[1,9,8,4,5],[6,7,6,9,10],[11,12,13,14,15], [16,17,18,19,20]]);
+    assert_equal(submatrix_set(test,[[9,8],[7,6]],1,2), [[1,2,3,4,5],[6,7,9,8,10],[11,12,7,6,15], [16,17,18,19,20]]);
+    assert_equal(submatrix_set(test,[[9,8],[7,6]],-1,-1), [[6,2,3,4,5],[6,7,8,9,10],[11,12,13,14,15], [16,17,18,19,20]]);
+    assert_equal(submatrix_set(test,[[9,8],[7,6]],n=4), [[1,2,3,4,9],[6,7,8,9,7],[11,12,13,14,15], [16,17,18,19,20]]);
+    assert_equal(submatrix_set(test,[[9,8],[7,6]],7,7), [[1,2,3,4,5],[6,7,8,9,10],[11,12,13,14,15], [16,17,18,19,20]]);
+    assert_equal(submatrix_set(ragged, [["a","b"],["c","d"]], 1, 1), [[1,2,3,4,5],[6,"a","b",9,10],[11,"c"], [16,17]]);
+    assert_equal(submatrix_set(test, [[]]), test);
+}
+test_submatrix_set();
+
 
 module test_array_group() {
     v = [1,2,3,4,5,6];
@@ -473,6 +533,12 @@ module test_array_dim() {
     assert(array_dim([[[1,2,3],[4,5,6]],[[7,8,9],[10,11,12]]], 0) == 2);
     assert(array_dim([[[1,2,3],[4,5,6]],[[7,8,9],[10,11,12]]], 2) == 3);
     assert(array_dim([[[1,2,3],[4,5,6]],[[7,8,9]]]) == [2,undef,3]);
+    assert(array_dim([1,2,3,4,5,6,7,8,9]) == [9]);
+    assert(array_dim([[1],[2],[3],[4],[5],[6],[7],[8],[9]]) == [9,1]);
+    assert(array_dim([]) == [0]);
+    assert(array_dim([[]]) == [1,0]);
+    assert(array_dim([[],[]]) == [2,0]);
+    assert(array_dim([[],[1]]) == [2,undef]);
 }
 test_array_dim();
 
@@ -486,7 +552,4 @@ module test_transpose() {
 test_transpose();
 
 
-cube();
-
-
 // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap
diff --git a/tests/test_geometry.scad b/tests/test_geometry.scad
index 0950e1c..e69a8ac 100644
--- a/tests/test_geometry.scad
+++ b/tests/test_geometry.scad
@@ -53,7 +53,7 @@ test_tri_functions();
 //test__general_plane_line_intersection();
 //test_plane_line_angle();
 //test_plane_line_intersection();
-//test_polygon_line_intersection();
+test_polygon_line_intersection();
 //test_plane_intersection();
 test_coplanar();
 test_points_on_plane();
@@ -542,6 +542,17 @@ module test_distance_from_plane() {
 *test_distance_from_plane();
 
 
+module test_polygon_line_intersection() {
+    poly1 = [[50,50,50], [50,-50,50], [-50,-50,50]];
+    assert_approx(polygon_line_intersection(poly1, [CENTER, UP]), [0,0,50]);
+    assert_approx(polygon_line_intersection(poly1, [CENTER, UP+RIGHT]), [50,0,50]);
+    assert_approx(polygon_line_intersection(poly1, [CENTER, UP+BACK+RIGHT]), [50,50,50]);
+    assert_approx(polygon_line_intersection(poly1, [[0,0,50], [1,0,50]]), [[[0,0,50], [50,0,50]]]);
+    assert_approx(polygon_line_intersection(poly1, [[0,0,0], [1,0,0]]), undef);
+}
+*test_polygon_line_intersection();
+
+
 module test_coplanar() {
     assert(coplanar([ [5,5,1],[0,0,1],[-1,-1,1] ]) == false);
     assert(coplanar([ [5,5,1],[0,0,0],[-1,-1,1] ]) == true);
diff --git a/tests/test_math.scad b/tests/test_math.scad
index c3eb601..9a25cec 100644
--- a/tests/test_math.scad
+++ b/tests/test_math.scad
@@ -781,6 +781,12 @@ test_back_substitute();
 
 
 
+module test_norm_fro(){
+  assert_approx(norm_fro([[2,3,4],[4,5,6]]), 10.29563014098700);
+
+} test_norm_fro();  
+
+
 module test_linear_solve(){
   M = [[-2,-5,-1,3],
        [3,7,6,2],
@@ -954,6 +960,38 @@ module test_real_roots(){
 test_real_roots();
 
 
+
+module test_quadratic_roots(){
+    assert_approx(quadratic_roots([1,4,4]),[[-2,0],[-2,0]]);
+    assert_approx(quadratic_roots([1,4,4],real=true),[-2,-2]);
+    assert_approx(quadratic_roots([1,-5,6],real=true), [2,3]);
+    assert_approx(quadratic_roots([1,-5,6]), [[2,0],[3,0]]);
+}
+test_quadratic_roots();
+
+
+module test_null_space(){
+    assert_equal(null_space([[3,2,1],[3,6,3],[3,9,-3]]),[]);
+
+    function nullcheck(A,dim) =
+      let(v=null_space(A))
+        len(v)==dim && is_zero(A*transpose(v),eps=1e-12);
+    
+   A = [[-1, 2, -5, 2],[-3,-1,3,-3],[5,0,5,0],[3,-4,11,-4]];
+   assert(nullcheck(A,1));
+
+   B = [
+        [  4,    1,    8,    6,   -2,    3],
+        [ 10,    5,   10,   10,    0,    5],
+        [  8,    1,    8,    8,   -6,    1],
+        [ -8,   -8,    6,   -1,   -8,   -1],
+        [  2,    2,    0,    1,    2,    1],
+        [  2,   -3,   10,    6,   -8,    1],
+       ];
+   assert(nullcheck(B,3));
+}
+test_null_space();
+
 module test_qr_factor() {
   // Check that R is upper triangular
   function is_ut(R) =
@@ -962,7 +1000,15 @@ module test_qr_factor() {
 
   // Test the R is upper trianglar, Q is orthogonal and qr=M
   function qrok(qr,M) =
-     is_ut(qr[1]) && approx(qr[0]*transpose(qr[0]), ident(len(qr[0]))) && approx(qr[0]*qr[1],M);
+     is_ut(qr[1]) && approx(qr[0]*transpose(qr[0]), ident(len(qr[0]))) && approx(qr[0]*qr[1],M) && qr[2]==ident(len(qr[2]));
+
+  // Test the R is upper trianglar, Q is orthogonal, R diagonal non-increasing and qrp=M
+  function qrokpiv(qr,M) =
+       is_ut(qr[1])
+    && approx(qr[0]*transpose(qr[0]), ident(len(qr[0])))
+    && approx(qr[0]*qr[1]*transpose(qr[2]),M)
+    && list_decreasing([for(i=[0:1:min(len(qr[1]),len(qr[1][0]))-1]) abs(qr[1][i][i])]);
+
   
   M = [[1,2,9,4,5],
        [6,7,8,19,10],
@@ -991,6 +1037,15 @@ module test_qr_factor() {
   assert(qrok(qr_factor([[7]]), [[7]]));
   assert(qrok(qr_factor([[1,2,3]]), [[1,2,3]]));
   assert(qrok(qr_factor([[1],[2],[3]]), [[1],[2],[3]]));
+
+
+  assert(qrokpiv(qr_factor(M,pivot=true),M));
+  assert(qrokpiv(qr_factor(select(M,0,3),pivot=true),select(M,0,3)));
+  assert(qrokpiv(qr_factor(transpose(select(M,0,3)),pivot=true),transpose(select(M,0,3))));
+  assert(qrokpiv(qr_factor(B,pivot=true),B));
+  assert(qrokpiv(qr_factor([[7]],pivot=true), [[7]]));
+  assert(qrokpiv(qr_factor([[1,2,3]],pivot=true), [[1,2,3]]));
+  assert(qrokpiv(qr_factor([[1],[2],[3]],pivot=true), [[1],[2],[3]]));
 }
 test_qr_factor();
 
diff --git a/tests/test_transforms.scad b/tests/test_transforms.scad
index a44e47c..3fb044b 100644
--- a/tests/test_transforms.scad
+++ b/tests/test_transforms.scad
@@ -232,7 +232,7 @@ module test_rot() {
         assert_equal(rot(a,p=pts2d), pts2d, info=str("rot(",a,",p=...), 2D"));
         assert_equal(rot(a,p=pts3d), pts3d, info=str("rot(",a,",p=...), 3D"));
     }
-    assert_equal(rot(90), [[0,-1,0,0],[1,0,0,0],[0,0,1,0],[0,0,0,1]])
+    assert_equal(rot(90), [[0,-1,0,0],[1,0,0,0],[0,0,1,0],[0,0,0,1]]);
     for (a=angs) {
         assert_equal(rot(a), affine3d_zrot(a), info=str("Z angle (only) = ",a));
         assert_equal(rot([a,0,0]), affine3d_xrot(a), info=str("X angle = ",a));
diff --git a/version.scad b/version.scad
index c2e6cda..262b891 100644
--- a/version.scad
+++ b/version.scad
@@ -8,7 +8,7 @@
 //////////////////////////////////////////////////////////////////////
 
 
-BOSL_VERSION = [2,0,409];
+BOSL_VERSION = [2,0,419];
 
 
 // Section: BOSL Library Version Functions
@@ -49,6 +49,7 @@ function bosl_version_str() = version_to_str(BOSL_VERSION);
 // Description:
 //   Given a version as a list, number, or string, asserts that the currently installed BOSL library is at least the given version.
 module bosl_required(target) {
+    no_children($children);
     assert(
         version_cmp(bosl_version(), target) >= 0,
         str(
diff --git a/vnf.scad b/vnf.scad
index 05fd82f..085c281 100644
--- a/vnf.scad
+++ b/vnf.scad
@@ -341,9 +341,18 @@ function vnf_vertex_array(
 // Arguments:
 //   vnf = A VNF structure, or list of VNF structures.
 //   convexity = Max number of times a line could intersect a wall of the shape.
-module vnf_polyhedron(vnf, convexity=2) {
+//   extent = If true, calculate anchors by extents, rather than intersection.  Default: true.
+//   cp = Centerpoint of VNF to use for anchoring when `extent` is false.  Default: `[0, 0, 0]`
+//   anchor = Translate so anchor point is at origin (0,0,0).  See [anchor](attachments.scad#anchor).  Default: `"origin"`
+//   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`
+module vnf_polyhedron(vnf, convexity=2, extent=true, cp=[0,0,0], anchor="origin", spin=0, orient=UP) {
     vnf = is_vnf_list(vnf)? vnf_merge(vnf) : vnf;
-    polyhedron(vnf[0], vnf[1], convexity=convexity);
+    cp = is_def(cp) ? cp : vnf_centroid(vnf);
+    attachable(anchor,spin,orient, vnf=vnf, extent=extent, cp=cp) {
+        polyhedron(vnf[0], vnf[1], convexity=convexity);
+        children();
+    }
 }