diff --git a/math/MathUtil.cs b/math/MathUtil.cs
index ef708a78..6c5723a6 100644
--- a/math/MathUtil.cs
+++ b/math/MathUtil.cs
@@ -10,6 +10,7 @@ public static class MathUtil
         public const double Deg2Rad = (Math.PI / 180.0);
         public const double Rad2Deg = (180.0 / Math.PI);
         public const double TwoPI = 2.0 * Math.PI;
+        public const double FourPI = 4.0 * Math.PI;
         public const double HalfPI = 0.5 * Math.PI;
         public const double ZeroTolerance = 1e-08;
         public const double Epsilon = 2.2204460492503131e-016;
diff --git a/math/Matrix3d.cs b/math/Matrix3d.cs
index 7fc052f2..f0fd6395 100644
--- a/math/Matrix3d.cs
+++ b/math/Matrix3d.cs
@@ -85,6 +85,18 @@ public Matrix3d(double m00, double m01, double m02, double m10, double m11, doub
         }
 
 
+        /// <summary>
+        /// Construct outer-product of u*transpose(v) of u and v
+        /// result is that Mij = u_i * v_j
+        /// </summary>
+        public Matrix3d(ref Vector3d u, ref Vector3d v)
+        {
+            Row0 = new Vector3d(u.x*v.x, u.x*v.y, u.x*v.z);
+            Row1 = new Vector3d(u.y*v.x, u.y*v.y, u.y*v.z);
+            Row2 = new Vector3d(u.z*v.x, u.z*v.y, u.z*v.z);
+        }
+
+
         public static readonly Matrix3d Identity = new Matrix3d(true);
         public static readonly Matrix3d Zero = new Matrix3d(false);
 
@@ -201,6 +213,14 @@ public void Multiply(ref Vector3d v, ref Vector3d vOut) {
         }
 
 
+
+        public double InnerProduct(ref Matrix3d m2)
+        {
+            return Row0.Dot(ref m2.Row0) + Row1.Dot(ref m2.Row1) + Row2.Dot(ref m2.Row2);
+        }
+
+
+
         public double Determinant {
 			get {
 				double a11 = Row0.x, a12 = Row0.y, a13 = Row0.z, a21 = Row1.x, a22 = Row1.y, a23 = Row1.z, a31 = Row2.x, a32 = Row2.y, a33 = Row2.z;
diff --git a/spatial/DMeshAABBTree.cs b/spatial/DMeshAABBTree.cs
index 20107b4f..6220377c 100644
--- a/spatial/DMeshAABBTree.cs
+++ b/spatial/DMeshAABBTree.cs
@@ -30,8 +30,8 @@ namespace g3
     /// </summary>
     public class DMeshAABBTree3 : ISpatial
     {
-        DMesh3 mesh;
-        int mesh_timestamp;
+        protected DMesh3 mesh;
+        protected int mesh_timestamp;
 
         public DMeshAABBTree3(DMesh3 m, bool autoBuild = false)
         {
@@ -1311,9 +1311,9 @@ public double TotalExtentSum()
         // storage for box nodes. 
         //   - box_to_index is a pointer into index_list
         //   - box_centers and box_extents are the centers/extents of the bounding boxes
-        DVector<int> box_to_index;
-        DVector<Vector3f> box_centers;
-        DVector<Vector3f> box_extents;
+        protected DVector<int> box_to_index;
+        protected DVector<Vector3f> box_centers;
+        protected DVector<Vector3f> box_extents;
 
         // list of indices for a given box. There is *no* marker/sentinel between
         // boxes, you have to get the starting index from box_to_index[]
@@ -1325,13 +1325,13 @@ public double TotalExtentSum()
         //       internal box, with index (-index_list[i])-1     (shift-by-one in case actual value is 0!)
         //   - if i > triangles_end and index_list[i] > 0, this is a two-child
         //       internal box, with indices index_list[i]-1 and index_list[i+1]-1
-        DVector<int> index_list;
+        protected DVector<int> index_list;
 
         // index_list[i] for i < triangles_end is a triangle-index list, otherwise box-index pair/single
-        int triangles_end = -1;
+        protected int triangles_end = -1;
 
         // box_to_index[root_index] is the root node of the tree
-        int root_index = -1;
+        protected int root_index = -1;
 
 
 
@@ -2040,7 +2040,7 @@ double box_distance_sqr(int iBox, Vector3d p)
         }
 
 
-        bool box_contains(int iBox, Vector3d p)
+        protected bool box_contains(int iBox, Vector3d p)
         {
             // [TODO] this could be way faster...
             Vector3d c = (Vector3d)box_centers[iBox];