Fun With a Matrix Math Library

[WayneHill]

I found a matrix math library in Delphi at this web site:

http://www.efg2.com/Lab/Library/Delp...ibrary.PAS.TXT

It is written by Earl F. Glynn and his copyright notice allows reprinting of his code for educational purposes.

http://www.efg2.com/Lab/Copyright.htm

// Graphics Math Library
//
// Copyright (C) 1982, 1985, 1992, 1995-1998 Earl F. Glynn, Overland Park, KS.
// All Rights Reserved.


UNIT GraphicsMathLibrary;  // Matrix/Vector Operations for 2D/3D Graphics

INTERFACE

  USES
    SysUtils;  {Exception}

  CONST
    sizeUndefined = 1;
    size2D        = 3;   // 'size' of 2D homogeneous vector or transform matrix
    size3D        = 4;   // 'size' of 3D homogeneous vector or transform matrix

  TYPE
    EVectorError = CLASS(Exception);
    EMatrixError = CLASS(Exception);

    TAxis       = (axisX, axisY, axisZ);
    TCoordinate = (coordCartesian, coordSpherical, coordCylindrical);
    TDimension  = (dimen2D, dimen3D);  // two- or three-dimensional TYPE
    TIndex      = 1..4;      // index of 'TMatrix' and 'TVector' TYPEs

    TMatrix     =            // transformation 'matrix'
      RECORD
        size:    TIndex;
        matrix:  ARRAY[TIndex,TIndex] OF DOUBLE
      END;

    Trotation   = (rotateClockwise, rotateCounterClockwise);

    // Normally the TVector TYPE is used to define 2D/3D homogenous
    // cartesian coordinates for graphics, i.e., (x,y,1) for 2D and
    // (x,y,z,1) for 3D.
    //
    // Cartesian coordinates can be converted to spherical (r, theta, phi),
    // or cylindrical coordinates (r,theta, z).  Spherical or cylindrical
    // coordinates can be converted back to cartesian coordinates.

    TVector     =
      RECORD
        size: TIndex;
        CASE INTEGER OF
          0:  (vector:  ARRAY[TIndex] OF DOUBLE);
          1:  (x:  DOUBLE;
               y:  DOUBLE;
               z:  DOUBLE;   // contains 'h' for 2D cartesian vector
               h:  DOUBLE)
        END;

                                                 // Vector Operations

  FUNCTION Vector2D  (CONST xValue, yValue:          DOUBLE):  TVector;
  FUNCTION Vector3D  (CONST xValue, yValue, zValue:  DOUBLE):  TVector;
  FUNCTION AddVectors (CONST u,v:  TVector):  TVector;
  FUNCTION Transform (CONST u:  TVector; CONST a:  TMatrix):  TVector;

  FUNCTION DotProduct  (CONST u,v:  TVector):  DOUBLE;
  FUNCTION CrossProduct(CONST u,v:  TVector):  TVector;

                                                 // Basic Matrix Operations

  FUNCTION Matrix2D (CONST m11,m12,m13,          // 2D "graphics" matrix
                           m21,m22,m23,
                           m31,m32,m33:  DOUBLE):  TMatrix;

  FUNCTION Matrix3D (CONST m11,m12,m13,m14,      // 3D "graphics" matrix
                           m21,m22,m23,m24,
                           m31,m32,m33,m34,
                           m41,m42,m43,m44:  DOUBLE):  TMatrix;

  FUNCTION MultiplyMatrices (CONST a,b:  TMatrix):  TMatrix;

  FUNCTION InvertMatrix  (CONST a,b:  TMatrix; VAR determinant:  DOUBLE):  TMatrix;


                                                 // Transformation Matrices

  FUNCTION RotateMatrix     (CONST dimension:  TDimension;
                             CONST xyz      :  TAxis;
                             CONST angle    :  DOUBLE;
                             CONST rotation :  Trotation):  TMatrix;

  FUNCTION ScaleMatrix      (CONST s:  TVector):  TMatrix;

  FUNCTION TranslateMatrix  (CONST t:  TVector):  TMatrix;

  FUNCTION ViewTransformMatrix (CONST coordinate:  TCoordinate;
       CONST azimuth {or x}, elevation {or y}, distance {or z}:  DOUBLE;
       CONST ScreenX, ScreenY, ScreenDistance:  DOUBLE):  TMatrix;


                                                 // conversions

  FUNCTION FromCartesian (CONST ToCoordinate:  TCoordinate; CONST u:  TVector):  TVector;
  FUNCTION ToCartesian   (CONST FromCoordinate:  TCoordinate; CONST u:  TVector):  TVector;

  FUNCTION ToDegrees(CONST angle {radians}:  DOUBLE):  DOUBLE {degrees};
  FUNCTION ToRadians(CONST angle {degrees}:  DOUBLE):  DOUBLE {radians};


                                                 // miscellaneous

  FUNCTION  Defuzz(CONST x:  DOUBLE):  DOUBLE;
  FUNCTION  GetFuzz:  DOUBLE;
  PROCEDURE SetFuzz(CONST x:  DOUBLE);


IMPLEMENTATION

  VAR
    fuzz : DOUBLE;


// *************************  Vector Operations *************************

  // This procedure defines two-dimensional homogeneous coordinates (x,y,1)
  // as a single 'vector' data element 'u'.  The 'size' of a two-dimensional
  // homogenous vector is 3.
  FUNCTION Vector2D  (CONST xValue, yValue:  DOUBLE):  TVector;

  BEGIN
    WITH RESULT DO
    BEGIN
      x    := xValue;
      y    := yValue;
      z    := 1.0;      // should be 'h' but let's use the 'z' slot to keep
      size := size2D    // subscripting possible
    END
  END {Vector2D};


  // This procedure defines three-dimensional homogeneous coordinates
  // (x,y,z,1) as a single 'vector' data element 'u'.  The 'size' of a
  // three-dimensional homogenous vector is 4.
  FUNCTION Vector3D  (CONST xValue, yValue, zValue:  DOUBLE):  TVector;
  BEGIN
    WITH RESULT DO
    BEGIN
      x    := xValue;
      y    := yValue;
      z    := zValue;
      h    := 1.0;       // homogeneous coordinate
      size := size3D
    END
  END {Vector3D};


  // AddVectors adds two vectors defined with homogeneous coordinates.
  FUNCTION AddVectors (CONST u,v:  TVector):  TVector;
    VAR
      i: TIndex;
  BEGIN
    IF  (u.size IN [size2D..size3D])  AND
        (v.size IN [size2D..size3D])  AND
        (u.size = v.size)
    THEN BEGIN
      RESULT.size := u.size;
      FOR i := 1 TO u.size-1 DO     {2D + 2D = 2D  or  3D + 3D = 3D}
      BEGIN
        RESULT.vector[i] := u.vector[i] + v.vector[i]
      END;
      RESULT.vector[u.size] := 1.0   {homogeneous coordinate}
    END
    ELSE raise EVectorError.Create('Vector Addition Mismatch')
  END {AddVectors};


  // 'Transform' multiplies a row 'vector' by a transformation 'matrix'
  // resulting in a new row 'vector'.  The 'size' of the 'vector' and 'matrix'
  // must agree.  To save execution time, the vectors are assumed to contain
  // a homogeneous coordinate.
  FUNCTION Transform (CONST u:  TVector; CONST a:  TMatrix):  TVector;
    VAR
      i,k :  TIndex;
      temp:  DOUBLE;
  BEGIN
    RESULT.size := a.size;
    IF  a.size = u.size
    THEN BEGIN
      FOR i := 1 TO a.size-1 DO
      BEGIN
        temp := 0.0;
        FOR k := 1 TO a.size DO
        BEGIN
          temp := temp + u.vector[k]*a.matrix[k,i];
        END;
        RESULT.vector[i] := Defuzz(temp)
      END;
      RESULT.vector[a.size] := 1.0 {assume homogeneous coordinate}
    END
    ELSE raise EMatrixError.Create('Transform multiply error')
  END {Transform};


  // Assume vector contains 'extra' homogeneous coordinate -- ignore it.
  FUNCTION DotProduct  (CONST u,v:  TVector):  DOUBLE;
    VAR
      i:  INTEGER;
  BEGIN
    IF  (u.size = v.size)
    THEN BEGIN
      RESULT := 0.0;
      FOR i := 1 TO u.size-1 DO
      BEGIN
        RESULT := RESULT + u.vector[i] * v.vector[i];
      END;
    END
    ELSE RAISE EMatrixError.Create('Vector dot product error')
  END {DotProduct};


  // Assume vector contains 'extra' homogeneous coordinate -- ignore it.
  FUNCTION CrossProduct(CONST u,v:  TVector):  TVector;
  BEGIN
    IF  (u.size = v.size) AND (u.size = size3D)
    THEN BEGIN
      RESULT := Vector3D( u.y*v.z - v.y*u.z,
                         -u.x*v.z + v.x*u.z,
                          u.x*v.y - v.x*u.y)
    END
    ELSE RAISE EMatrixError.Create('Vector cross product error')
  END {CrossProduct};


// *********************** Basic Matrix Operations **********************

  FUNCTION Matrix2D (CONST m11,m12,m13, m21,m22,m23, m31,m32,m33:  DOUBLE):
                     TMatrix;
  BEGIN
    WITH RESULT DO
    BEGIN
      matrix[1,1] := m11; matrix[1,2] := m12; matrix[1,3] := m13;
      matrix[2,1] := m21; matrix[2,2] := m22; matrix[2,3] := m23;
      matrix[3,1] := m31; matrix[3,2] := m32; matrix[3,3] := m33;
      size := size2D
    END
  END {Matrix2D};


  FUNCTION Matrix3D (CONST m11,m12,m13,m14, m21,m22,m23,m24,
                           m31,m32,m33,m34, m41,m42,m43,m44:  DOUBLE):  TMatrix;
  BEGIN
    WITH RESULT DO
    BEGIN
      matrix[1,1] := m11; matrix[1,2] := m12;
      matrix[1,3] := m13; matrix[1,4] := m14;

      matrix[2,1] := m21; matrix[2,2] := m22;
      matrix[2,3] := m23; matrix[2,4] := m24;

      matrix[3,1] := m31; matrix[3,2] := m32;
      matrix[3,3] := m33; matrix[3,4] := m34;

      matrix[4,1] := m41; matrix[4,2] := m42;
      matrix[4,3] := m43; matrix[4,4] := m44;
      size := size3D
    END
  END {Matrix3D};


  // Compound geometric transformation matrices can be formed by multiplying
  // simple transformation matrices.  This procedure only multiplies together
  // matrices for two- or three-dimensional transformations, i.e., 3x3 or 4x4
  // matrices.  The multiplier and multiplicand must be of the same dimension.
  FUNCTION MultiplyMatrices (CONST a,b:  TMatrix):  TMatrix;
    VAR
      i,j,k:  TIndex;
      temp :  DOUBLE;
  BEGIN
    RESULT.size := a.size;
    IF  a.size = b.size
    THEN

      FOR i := 1 TO a.size DO
      BEGIN
        FOR j := 1 TO a.size DO
        BEGIN

          temp := 0.0;
          FOR k := 1 TO a.size DO
          BEGIN
            temp := temp + a.matrix[i,k]*b.matrix[k,j];
          END;
          RESULT.matrix[i,j] := Defuzz(temp)

        END
      END

    ELSE EMatrixError.Create('MultiplyMatrices error')
  END {MultiplyMatrices};


  // This procedure inverts a general transformation matrix.  The user need
  // not form an inverse geometric transformation by keeping a product of
  // the inverses of simple geometric transformations:  translations, rotations
  // and scaling.  A determinant of zero indicates no inverse is possible for
  // a singular matrix.
  FUNCTION InvertMatrix  (CONST a,b:  TMatrix; VAR determinant:  DOUBLE):  TMatrix;
    VAR
      c        :  TMatrix;
      i,i_pivot:  TIndex;
      i_flag   :  ARRAY[TIndex] OF BOOLEAN;
      j,j_pivot:  TIndex;
      j_flag   :  ARRAY[TIndex] OF BOOLEAN;
      modulus  :  DOUBLE;
      n        :  TIndex;
      pivot    :  DOUBLE;
      pivot_col:  ARRAY[TIndex] OF TIndex;
      pivot_row:  ARRAY[TIndex] OF TIndex;
      temporary:  DOUBLE;
  BEGIN
    c := a;                         // The matrix inversion algorithm used here
    WITH c DO                       // is similar to the "maximum pivot strategy"
    BEGIN                           // described in "Applied Numerical Methods"
      FOR i := 1 TO size DO         // by Carnahan, Luther and Wilkes,
      BEGIN                         // pp. 282-284.
        i_flag[i] := TRUE;
        j_flag[i] := TRUE
      END;
      modulus := 1.0;
      i_pivot := 1;  // avoid initialization warning
      j_pivot := 1;  // avoid initialization warning

      FOR n := 1 TO size DO
      BEGIN
        pivot := 0.0;
        IF   ABS(modulus) > 0.0
        THEN BEGIN
          FOR i := 1 TO size DO
            IF  i_flag[i]
            THEN

              FOR j := 1 TO size DO
                IF   j_flag[j]
                THEN
                  IF   ABS(matrix[i,j]) > ABS(pivot)
                  THEN BEGIN
                    pivot := matrix[i,j];   // largest value on which to pivot
                    i_pivot := i;           // indices of pivot element
                    j_pivot := j
                  END;

          IF   Defuzz(pivot) = 0    // If pivot is too small, consider
          THEN modulus := 0         // the matrix to be singular
          ELSE BEGIN
            pivot_row[n] := i_pivot;
            pivot_col[n] := j_pivot;
            i_flag[i_pivot] := FALSE;
            j_flag[j_pivot] := FALSE;
            FOR i := 1 TO size DO
              IF   i <> i_pivot
              THEN
                FOR j := 1 TO size DO  // pivot column unchanged for elements
                  IF   j <> j_pivot    // not in pivot row or column ...
                  THEN matrix[i,j] := (matrix[i,j]*matrix[i_pivot,j_pivot] -
                                    matrix[i_pivot,j]*matrix[i,j_pivot])
                                    / modulus;  // 2x2 minor / modulus
            FOR j := 1 TO size DO
              IF   j <> j_pivot        // change signs of elements in pivot row
              THEN matrix[i_pivot,j] := -matrix[i_pivot,j];
            temporary := modulus;      // exchange pivot element and modulus
            modulus := matrix[i_pivot,j_pivot];
            matrix[i_pivot,j_pivot] := temporary
          END
        END
      END {FOR n}
    END {WITH};

    determinant := Defuzz(modulus);
    IF  determinant <> 0
    THEN BEGIN
      RESULT.size := c.size;       // The matrix inverse must be unscrambled
      FOR i := 1 TO c.size DO      // if pivoting was not along main diagonal.
        FOR j := 1 TO c.size DO
          RESULT.matrix[pivot_row[i],pivot_col[j]] := Defuzz(c.matrix[i,j]/determinant)
    END
    ELSE EMatrixError.Create('InvertMatrix error')

  END {InvertMatrix};


// ***********************  Transformation Matrices  ********************


  // This procedure defines a matrix for a two- or three-dimensional rotation.
  // To avoid possible confusion in the sense of the rotation, 'rotateClockwise'
  // or 'roCounterlcockwise' must always be specified along with the axis
  // of rotation. Two-dimensional rotations are assumed to be about the z-axis
  // in the x-y plane.
  //
  // A rotation about an arbitrary axis can be performed with the following
  // steps:
  //   (1) Translate the object into a new coordinate system where (x,y,z)
  //       maps into the origin (0,0,0).
  //   (2) Perform appropriate rotations about the x and y axes of the
  //       coordinate system so that the unit vector (a,b,c) is mapped into
  //       the unit vector along the z axis.
  //   (3) Perform the desired rotation about the z-axis of the new
  //       coordinate system.
  //   (4) Apply the inverse of step (2).
  //   (5) Apply the inverse of step (1).
  FUNCTION RotateMatrix     (CONST dimension:  TDimension;
                             CONST xyz      :  TAxis;
                             CONST angle    :  DOUBLE;
                             CONST rotation :  Trotation):  TMatrix;
    VAR
      cosx     :  DOUBLE;
      sinx     :  DOUBLE;
      TempAngle:  DOUBLE;

  BEGIN
    TempAngle := angle;  // Use TempAngle since "angle" is CONST parameter

    IF  rotation = rotateCounterClockwise
    THEN TempAngle := -TempAngle;

    cosx := Defuzz( COS(TempAngle) );
    sinx := Defuzz( SIN(TempAngle) );

    CASE dimension OF
      dimen2D:
        CASE xyz OF
          axisX,axisY:  EMatrixError.Create('Invalid 2D rotation matrix.  Specify axisZ');

          axisZ:  RESULT := Matrix2D ( cosx, -sinx,     0,
                                       sinx,  cosx,     0,
                                          0,     0,     1)
        END;

      dimen3D:
        CASE xyz OF
          axisX:  RESULT := Matrix3D (    1,     0,     0, 0,
                                          0,  cosx, -sinx, 0,
                                          0,  sinx,  cosx, 0,
                                          0,     0,     0, 1);

          axisY:  RESULT := Matrix3D ( cosx,     0,  sinx, 0,
                                          0,     1,     0, 0,
                                      -sinx,     0,  cosx, 0,
                                          0,     0,     0, 1);

          axisZ:  RESULT := Matrix3D ( cosx, -sinx,     0, 0,
                                       sinx,  cosx,     0, 0,
                                          0,     0,     1, 0,
                                          0,     0,     0, 1);
        END
    END
  END {RotateMatrix};


  // 'ScaleMatrix' accepts a 'vector' containing the scaling factors for
  //  each of the dimensions and creates a scaling matrix.  The size
  //  of the vector dictates the size of the resulting matrix.
  FUNCTION ScaleMatrix      (CONST s:  TVector):  TMatrix;
  BEGIN
    CASE s.size OF
      size2D: RESULT := Matrix2D (s.x,   0,   0,
                                    0, s.y,   0,
                                    0,   0,   1);

      size3D: RESULT := Matrix3D (s.x,   0,   0,  0,
                                    0, s.y,   0,  0,
                                    0,   0, s.z,  0,
                                    0,   0,   0,  1)
    END
  END {ScaleMatrix};


  // 'TranslateMatrix' defines a translation transformation matrix.  The
  // components of the vector 't' determine the translation components.
  // (Note:  'Translate' here is from kinematics in physics.)
  FUNCTION TranslateMatrix  (CONST t:  TVector):  TMatrix;
  BEGIN
    CASE t.size OF
      size2D: RESULT := Matrix2D (  1,   0, 0,
                                    0,   1, 0,
                                  t.x, t.y, 1);

      size3D: RESULT := Matrix3D (  1,   0,   0,  0,
                                    0,   1,   0,  0,
                                    0,   0,   1,  0,
                                  t.x, t.y, t.z,  1)
    END
  END {TranslateMatrix};


  // 'ViewTransformMatrix' creates a transformation matrix for changing
  // from world coordinates to eye coordinates. The location of the 'eye'
  // from the 'object' is given in spherical (azimuth,elevation,distance)
  // coordinates or Cartesian (x,y,z) coordinates.  The size of the screen
  // is 'ScreenX' units horizontally and 'ScreenY' units vertically.  The
  // eye is 'ScreenDistance' units from the viewing screen.  A large ratio
  // 'ScreenDistance/ScreenX (or ScreenY)' specifies a narrow aperature
  // -- a telephoto view.  Conversely, a small ratio specifies a large
  // aperature -- a wide-angle view.  This view transform matrix is very
  // useful as the default three-dimensional transformation matrix.  Once
  // set, all points are automatically transformed.
  FUNCTION ViewTransformMatrix (CONST coordinate:  TCoordinate;
       CONST azimuth {or x}, elevation {or y}, distance {or z}:  DOUBLE;
       CONST ScreenX, ScreenY, ScreenDistance:  DOUBLE):  TMatrix;

    CONST
      HalfPI   =  PI / 2.0;

    VAR
      a         :  TMatrix;
      b         :  TMatrix;
      cosm      :  DOUBLE;        // COS(-angle)
      hypotenuse:  DOUBLE;
      sinm      :  DOUBLE;        // SIN(-angle)
      temporary :  DOUBLE;
      u         :  TVector;
      x         :  DOUBLE  ABSOLUTE azimuth;     // x and azimuth are synonyms
      y         :  DOUBLE  ABSOLUTE elevation;   // synonyms
      z         :  DOUBLE  ABSOLUTE distance;    // synonyms

  BEGIN
    CASE coordinate OF
      coordCartesian:  u := Vector3D (-x, -y, -z);

      coordSpherical:
        BEGIN
          temporary := -distance * COS(elevation);
          u := Vector3D (temporary * COS(azimuth - HalfPI),
                         temporary * SIN(azimuth - HalfPI),
                        -distance  * SIN(elevation));
        END
    END;

    a := TranslateMatrix(u);      // translate origin to 'eye'
    b := RotateMatrix (dimen3D, axisX, HalfPI, rotateClockwise);
    a := MultiplyMatrices(a,b);

    CASE coordinate OF
      coordCartesian:
        BEGIN
          temporary := SQR(x) + SQR(y);
          hypotenuse := SQRT(temporary);
          cosm := -y/hypotenuse;
          sinm :=  x/hypotenuse;

          b := Matrix3D ( cosm, 0, sinm, 0,
                             0, 1,    0, 0,
                         -sinm, 0, cosm, 0,
                             0, 0,    0, 1);

          a := MultiplyMatrices (a,b);
          cosm := hypotenuse;
          hypotenuse := SQRT(temporary + SQR(z));
          cosm := cosm/hypotenuse;
          sinm := -z/hypotenuse;

          b := Matrix3D (    1,    0,     0,  0,
                             0, cosm, -sinm,  0,
                             0, sinm,  cosm,  0,
                             0,    0,     0,  1)
        END;
      coordSpherical:
        BEGIN
          b := RotateMatrix (dimen3D,axisY,-azimuth,rotateCounterClockwise);
          a := MultiplyMatrices(a,b);
          b := RotateMatrix (dimen3D,axisX,elevation,rotateCounterClockwise);
        END
    END {CASE};

    a := MultiplyMatrices (a,b);
    u := Vector3D (ScreenDistance/(0.5*ScreenX),
              ScreenDistance/(0.5*ScreenY),-1.0);
    b := ScaleMatrix (u);  // reverse sense of z-axis; screen transformation

    RESULT := MultiplyMatrices (a,b)
  END {ViewTransformMatrix};


// ***************************   Conversions   **************************

  // This function converts the vector parameter from Cartesian
  // coordinates to the specified type of coordinates.
  FUNCTION FromCartesian (CONST ToCoordinate:  TCoordinate; CONST u:  TVector):  TVector;
    VAR
      phi  :  DOUBLE;
      r    :  DOUBLE;
      temp :  DOUBLE;
      theta:  DOUBLE;

  BEGIN
    IF  ToCoordinate = coordCartesian
    THEN RESULT := u
    ELSE BEGIN
      RESULT.size := u.size;

      IF   (u.size = size3D) AND
           (ToCoordinate = coordSpherical)
      THEN BEGIN                    // spherical 3D
        temp := SQR(u.x)+SQR(u.y);  // (x,y,z) -> (r,theta,phi)
        r := SQRT(temp+SQR(u.z));
        IF   Defuzz(u.x) = 0.0
        THEN theta := PI/4
        ELSE theta := ARCTAN(u.y/u.x);
        IF   Defuzz(u.z) = 0.0
        THEN phi := PI/4
        ELSE phi := ARCTAN(SQRT(temp)/u.z);
        RESULT.x := r;
        RESULT.y := theta;
        RESULT.z := phi
      END
      ELSE BEGIN              // cylindrical 2D/3D or spherical 2D
                              // (x,y) -> (r,theta)  or  (x,y,z) -> (r,theta,z)
        r := SQRT( SQR(u.x) + SQR(u.y) );
        IF   Defuzz(u.x) = 0.0
        THEN theta := PI/4
        ELSE theta := ARCTAN(u.y/u.x);
        RESULT.x := r;
        RESULT.y := theta
      END

    END
  END {FromCartesian};


  // This function converts the vector parameter from specified coordinates
  // into Cartesian coordinates.
  FUNCTION ToCartesian   (CONST FromCoordinate:  TCoordinate; CONST u:  TVector):  TVector;
    VAR
      phi   :  DOUBLE;
      r     :  DOUBLE;
      sinphi:  DOUBLE;
      theta :  DOUBLE;

  BEGIN
    RESULT := u;

    IF  FromCoordinate = coordCartesian
    THEN RESULT := u
    ELSE BEGIN
      RESULT.size := u.size;

      IF   (u.size = size3D) AND
           (FromCoordinate = coordSpherical)
      THEN BEGIN       // spherical 3D
        r :=  u.x;     //  (r,theta,phi) -> (x,y,z)
        theta := u.y;
        phi := u.z;
        sinphi := SIN(phi);
        RESULT.x := r * COS(theta) * sinphi;
        RESULT.y := r * SIN(theta) * sinphi;
        RESULT.z := r * COS(phi)
      END
      ELSE BEGIN       // cylindrical 2D/3D or spherical 2D
        r :=  u.x;     // (r,theta) -> (x,y)  or  (r,theta,z) -> (x,y,z)
        theta := u.y;
        RESULT.x := r * COS(theta);
        RESULT.y := r * SIN(theta)
      END
    END
  END {ToCartesian};


  // Convert angle in radians to degrees.
  FUNCTION ToDegrees(CONST angle {radians}:  DOUBLE):  DOUBLE {degrees};
  BEGIN
    RESULT := 180.0/PI * angle
  END {ToRadians};


  // Convert angle in degrees to radians.
  FUNCTION ToRadians (CONST angle:  DOUBLE):  DOUBLE;
  BEGIN
    RESULT := PI/180.0 * angle
  END {ToRadians};


// ***************************  Miscellaneous  **************************

  // 'Defuzz' is used for comparisons and to avoid propagation of 'fuzzy',
  //  nearly-zero values.  DOUBLE calculations often result in 'fuzzy' values.
  //  The term 'fuzz' was adapted from the APL language.
  FUNCTION  Defuzz(CONST x:  DOUBLE):  DOUBLE;
  BEGIN
    IF  ABS(x) < fuzz
    THEN RESULT := 0.0
    ELSE RESULT := x
  END {Defuzz};


  FUNCTION  GetFuzz:   DOUBLE;
  BEGIN
    RESULT := fuzz
  END {GetFuzz};


  PROCEDURE SetFuzz(CONST x:  DOUBLE);
  BEGIN
    fuzz := x
  END {SetFuzz};


INITIALIZATION
  SetFuzz(1.0E-6)
END {GraphicsMath UNIT}.

I have written a very simple method for using this matrix library.
Delphi code:

unit Unit1;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, StdCtrls;

type
  TForm1 = class(TForm)
    Button1: TButton;
    procedure Button1Click(Sender: TObject);
  private
    { Private declarations }
  public
    { Public declarations }
  end;

var
  Form1: TForm1;

implementation

uses GraphicsMathLibrary;

{$R *.dfm}

procedure TForm1.Button1Click(Sender: TObject);
var
  mMatrix: TMatrix;
  vData: TVector;
  vVector: TVector;
begin

  // Change X,Y,Z point to a vector unit
  vData := Vector3D(1, 0, 2.415);

  // Setup Matrix for rotation on X Axis at 10 degrees counterclockwise
  mMatrix := RotateMatrix(dimen3D, axisX, ToRadians(10), rotateCounterClockwise);

  // Run the Vector unit through the X Rotate matrix
  vData := Transform(vData, mMatrix);

  // Show the vData output
  showmessage(floattostr(vData.x) + ',' + floattostr(vData.y) + ',' + floattostr(vData.z));

  // Setup Matrix for rotation on Y Axis at 20 degrees counterclockwise
  mMatrix := RotateMatrix(dimen3D, axisY, ToRadians(20), rotateCounterClockwise);

  // Use X Rotated Data
  vData := Vector3D(vData.x, vData.y, vData.z);

  // Run the Vector unit through the Y Rotate matrix
  vData := Transform(vData, mMatrix);

  // Show the vData output
  showmessage(floattostr(vData.x) + ',' + floattostr(vData.y) + ',' + floattostr(vData.z));

  // Setup Matrix for rotation on Z Axis at 30 degrees counterclockwise
  mMatrix := RotateMatrix(dimen3D, axisZ, ToRadians(30), rotateCounterClockwise);

  // Use Y Rotated Data
  vData := Vector3D(vData.x, vData.y, vData.z);

  // Run the Vector unit through the Z Rotate matrix
  vData := Transform(vData, mMatrix);

  // Show the vData output
  showmessage(floattostr(vData.x) + ',' + floattostr(vData.y) + ',' + floattostr(vData.z));

  // Setup a vVector for Translation move Y of 2
  vVector := Vector3D(0, 2, 0);

  // Setup Tramslation Matrix with vVector unit
  mMatrix := TranslateMatrix(vVector);

  // Run the vData unit through the Translation Matrix
  vData := Transform(vData, mMatrix);

  // Show the VData output
  showmessage(floattostr(vData.x) + ',' + floattostr(vData.y) + ',' + floattostr(vData.z));

  // Setup Scale Matrix with vVector unit to 2 times scale
  vVector := Vector3D(2, 2, 2);

  // Setup Scale Matrix with vVector unit
  mMatrix := ScaleMatrix(vVector);

  // Run the Vector unit through the Z Rotate matrix
  vData := Transform(vData, mMatrix);

  // Show the VData output
  showmessage(floattostr(vData.x) + ',' + floattostr(vData.y) + ',' + floattostr(vData.z));

end;

end.