From cc644c7a77b16448541d31903e4d530f363e8ada Mon Sep 17 00:00:00 2001
From: Adrian Mariano <avm4@cornell.edu>
Date: Sat, 7 Mar 2020 14:22:53 -0500
Subject: [PATCH] Added plane intersection and plane_from_normal and fixed docs
 for plane definition.

---
 geometry.scad | 42 +++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 41 insertions(+), 1 deletion(-)

diff --git a/geometry.scad b/geometry.scad
index f8ab499..cf315e6 100644
--- a/geometry.scad
+++ b/geometry.scad
@@ -548,7 +548,7 @@ function triangle_area(a,b,c) =
 //   plane3pt(p1, p2, p3);
 // Description:
 //   Generates the cartesian equation of a plane from three non-collinear points on the plane.
-//   Returns [A,B,C,D] where Ax+By+Cz+D=0 is the equation of a plane.
+//   Returns [A,B,C,D] where Ax + By + Cz = D is the equation of a plane.
 // Arguments:
 //   p1 = The first point on the plane.
 //   p2 = The second point on the plane.
@@ -583,6 +583,46 @@ function plane3pt_indexed(points, i1, i2, i3) =
 	) plane3pt(p1,p2,p3);
 
 
+// Function: plane_intersection()
+// Usage:
+//   plane_intersection(plane1, plane2, [plane3])
+// 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.  
+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(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),
+      dd=echo(point=point, normal=normal)
+   )
+   [point, point+normal];
+
+
+// Function: plane_from_normal()
+// Usage:
+//   plane_from_normal(normal, pt)
+// Description:
+//   Returns a plane defined by a normal vector and a point.
+// 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) =
+  concat(normal, [normal*pt]);
+
+
 // Function: plane_from_pointslist()
 // Usage:
 //   plane_from_pointslist(points);