-- COGSpace.mesa: Basic geometrical operations in cartesian three-space
-- last modified by Stolfi - September 29, 1982 8:09 pm
-- compile COGSpace
DIRECTORY
COGRandom USING [Toss],
Real USING [SqRt];
COGSpace: CEDAR DEFINITIONS
IMPORTS COGRandom, Real =
BEGIN OPEN Rand: COGRandom, Real;
-- BASIC DATA TYPES - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Point: TYPE = RECORD [x, y, z: REAL];
Vector: TYPE = Point;
Box: TYPE = RECORD [min, max: Point];
-- PROCEDURES FOR POINTS AND VECTORS - - - - - - - - - - - - - - - - - - - - - - - - - -
Add: PROC [a, b: Vector] RETURNS [Vector] = INLINE
-- Adds two vectors (or points).
{RETURN [[a.x+b.x, a.y+b.y, a.z+b.z]]};
Sub: PROC [a, b: Vector] RETURNS [Vector] = INLINE
-- Subtracts vector b from vector a.
{RETURN [[a.x-b.x, a.y-b.y, a.z-b.z]]};
Mul: PROC [a: REAL, u: Vector] RETURNS [w: Vector] = INLINE
-- Multiplies the vector u by a real number a.
{RETURN [[u.x*a, u.y*a, u.z*a]]};
Way: PROC [alpha: REAL, p, q: Point] RETURNS [Point] = INLINE
-- Computes the point alpha the way from u to v.
{beta: REAL = (1.0-alpha);
RETURN [[beta*p.x+alpha*q.x, beta*p.y+alpha*q.y, beta*p.z+alpha*q.z]]};
Dot: PROC [u, v: Vector] RETURNS [REAL] = INLINE
-- Dot product of u and v.
{RETURN [u.x*v.x+u.y*v.y+u.z*v.z]};
XYRot: PROC [u: Vector, c, s: REAL] RETURNS [ur: Vector] = INLINE
-- Rotates the vector u around the z-axis by an angle whose cosine is c and whose sine is s.
{RETURN [[c*u.x-s*u.y, c*u.y+s*u.x, u.z]]};
YZRot: PROC [u: Vector, c, s: REAL] RETURNS [ur: Vector] = INLINE
-- Rotates the vector u around the x-axis by an angle whose cosine is c and whose sine is s.
{RETURN [[u.x, c*u.y-s*u.z, c*u.z+s*u.y]]};
ZXRot: PROC [u: Vector, c, s: REAL] RETURNS [ur: Vector] = INLINE
-- Rotates the vector u around the y-axis by an angle whose cosine is c and whose sine is s.
{RETURN [[c*u.x+s*u.z, u.y, c*u.z-s*u.x]]};
DoubleRot: PROC [u, v: Vector, c, s: REAL] RETURNS [ur, vr: Vector];
-- Rotates the vectors u in their own plane by an angle whose cosine is c and whose sine is s. Note: If the lengths of u and v are different, the angle is measured "elliptically".
NMirror: PROC [u, n: Vector] RETURNS [um: Vector] = INLINE
-- Mirrors the vector u about the plane whose normal is n (which should have nonzero length).
{RETURN [Sub [u, Mul [2.0*Dot[u, n]/Dot[n, n], n]]]};
Transf: PROC [u, xb, yb, zb: Vector] RETURNS [ut: Vector];
-- Post-multiplies the vector u by the 3x3 array whose lines are xv, yv and zv. That is, returns the vector u.x*xb+u.y*yb+u.z*zb.

Length
: PROC [v: Vector] RETURNS [REAL] = TRUSTED INLINE
-- Length of vector v.
{RETURN [SqRt[v.x*v.x+v.y*v.y+v.z*v.z]]};
LengthSq: PROC [v: Vector] RETURNS [REAL] = INLINE
-- Length of vector v, squared (saves a square root extraction).
{RETURN [v.x*v.x+v.y*v.y+v.z*v.z]};
Dist: PROC [p, q: Point] RETURNS [REAL] = INLINE
-- Distance between points p and q. Returns LargestNumber if one point is infinite, and TrappingNaN if both are infinite or one is indeterminate.
{RETURN [Length[Sub[p,q]]]};
Random: PROC [box: Box] RETURNS [v: Vector] = INLINE
-- Returns a random vector in the given box.
{RETURN [
[Rand.Toss [box.min.x, box.max.x],
Rand.Toss [box.min.y, box.max.y],
Rand.Toss [box.min.z, box.max.z]] ]};
END...