// Formulae for geometry. // Find the length of the hypotenuse of a right triangle given sides a and b. hypotenuseLength[a, b] := sqrt[a^2 + b^2] // Calculate area of arbitrary triangle given 3 sides of known length. // Calculate using Heron's formula. triangleArea[a, b, c] := { // Semiperimeter s s = 1/2 (a + b + c) area = sqrt[s (s-a) (s-b) (s-c)] return area } // Find the distance between a point (x0,y0) and a line specified by // two points (x1, y1), (x2, y2) pointToLineDistance[x0, y0, x1, y1, x2, y2] := { abs[ (x2-x1)(y1-y0) - (x1-x0)(y2-y1)] / sqrt[(x2-x1)^2+(y2-y1)^2 ] } // Find the intersection between two lines specified by two points on each // line. The first line is specified by the points // (x1, y1) and (x2, y2) // and the second line is specified by the points // (x3, y3) and (x4, y4) // // This returns the intersection point as [x, y] or undef if the lines // are paralllel. lineIntersection[x1, y1, x2, y2, x3, y3, x4, y4] := { det = (x1 - x2)(y3 - y4) - (y1 - y2)(x3 - x4) if det == 0 return undef t1 = (x1 y2 - y1 x2) t2 = (x3 y4 - y3 x4) px = (t1 (x3 - x4) - t2 (x1 - x2)) / det py = (t1 (y3 - y4) - t2 (y1 - y2)) / det return [px, py] } // Find the intersection between a plane and a line. // The plane is defined by three points 1,2,3 and the line is defined // by two points 4,5. Each point's coordinates are given by (x1,y1,z1), // (x2,y2,z2)... // See: // https://mathworld.wolfram.com/Line-PlaneIntersection.html // // returns [x,y,z,t] of which [x,y,z] is the intersection point (if one // exists) and t is a parameter indicating where on the parameterized line // equation it intersects with the plane: // // x(t) = x4 + (x5-x4) t // y(t) = y4 + (y5-y4) t // z(t) = z4 + (z5-z4) t planeLineIntersection[x1,y1,z1, x2,y2,z2, x3,y3,z3, x4,y4,z4, x5,y5,z5] := { num = x3 y2 z1 - x4 y2 z1 - x2 y3 z1 + x4 y3 z1 + x2 y4 z1 - x3 y4 z1 - x3 y1 z2 + x4 y1 z2 + x1 y3 z2 - x4 y3 z2 - x1 y4 z2 + x3 y4 z2 + x2 y1 z3 - x4 y1 z3 - x1 y2 z3 + x4 y2 z3 + x1 y4 z3 - x2 y4 z3 - x2 y1 z4 + x3 y1 z4 + x1 y2 z4 - x3 y2 z4 - x1 y3 z4 + x2 y3 z4 det = (z5 - z4) (x1 y2 - x1 y3 - x2 y1 + x2 y3 + x3 y1 - x3 y2) - (y5 - y4) (x1 z2 - x1 z3 - x2 z1 + x2 z3 + x3 z1 - x3 z2) + (x5 - x4) (y1 z2 - y1 z3 - y2 z1 + y2 z3 + y3 z1 - y3 z2) // println["num=\$num, det=\$det"] if det == 0 if num == 0 // Line is coplanar with plane (embedded in plane, // infinite solutions along entire line) return [undef, undef, undef, 0] else return [undef, undef, undef, undef] // Does not intersect, no solutions t = - num/det x = x4 + (x5-x4) t y = y4 + (y5-y4) t z = z4 + (z5-z4) t return [x,y,z,t] } // Returns true if a point is inside a polygon. // Adapted from // http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html // // points is an array of two-dimensional arrays of points. // x,y is the point to test. // // Returns: // true if [x,y] is inside the polygon, false otherwise. pointInPolygon[points, x, y] := { c = false i = 0 npol = length[points] j = npol-1 while i < npol { xi = points@i@0 xj = points@j@0 yi = points@i@1 yj = points@j@1 if ((((yi<=y) and (y inside Ax + By + Cz + D > 0 ==> outside The normal to the plane (pointing outward) is [A,B,C]. This does not test that the three points are non-colinear. */ planeCoefficients[p1, p2 ,p3] := { [A, B, C] = crossProduct[p1, p2, p3] D = (-A * p1@0) + (-B * p1@1) + (-C * p1@2) return [A, B, C, D] } "Geometry.frink included OK"