* Add system.math.vectors for Delphi compatibility.

packages/rtl-objpas/fpmake.pp

@@ -71,6 +71,8 @@ begin
     
     T:=P.Targets.AddUnit('system.actions.pp',UItypesOSes);
       T.Dependencies.AddUnit('system.uitypes');
+    T:=P.Targets.AddUnit('system.math.vectors.pp',UItypesOSes);
+      T.Dependencies.AddUnit('system.uitypes');
 
     T:=P.Targets.AddUnit('strutils.pp',StrUtilsOses);
       T.ResourceStrings:=true;

packages/rtl-objpas/src/inc/system.math.vectors.pp

@@ -0,0 +1,2407 @@
+{
+    This file is part of the Free Pascal run time library.
+    Copyright (c) 2023 by Michael Van Canneyt
+    member of the Free Pascal development team
+
+    See the file COPYING.FPC, included in this distribution,
+    for details about the copyright.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+ **********************************************************************}
+{
+  This unit contains modified version of some algorithms from the GLScene
+  GLS.VectorGeometry.pas.
+  The GLScene unit is covered under the MPL 1.0. license.
+}
+
+unit System.Math.Vectors;
+
+{$mode objfpc}
+{$modeswitch advancedrecords}
+
+interface
+
+uses
+{$IFDEF FPC_DOTTEDUNITS}
+  System.Types, System.Math;
+{$ELSE}
+  Types, Math;
+{$ENDIF}
+
+
+type
+  TEpsilon = record
+  const
+    Matrix = 1E-5;
+    Vector = 1E-4;
+    Scale = 1E-4;
+    FontSize = 1E-2;
+    Position = 1E-3;
+    Angle = 1E-4;
+  end;
+
+  TVector3DType = array [0..3] of Single;
+
+  TVectorArray = array [0..2] of Single;
+
+  TPolygon = array of TPointF;
+
+  TCubicBezier = array [0..3] of TPointF;
+
+  tagVECTOR = record
+    case Integer of
+      0: (V: TVectorArray;);
+      1: (X: Single;
+          Y: Single;
+          W: Single;);
+  end;
+
+  PVector = ^TVector;
+
+  { TVector }
+
+  TVector = record
+    class function Create(const aX, aY: Single; const aW: Single = 1.0): TVector; overload; static; inline;
+    class function Create(const aPoint: TPointF): TVector; overload; static; inline;
+    class function Zero: TVector; inline; static;
+
+    class operator +(const aVector1, aVector2: TVector): TVector;
+    class operator /(const aVector: TVector; const aFactor: Single): TVector;
+    class operator =(const aVector1, aVector2: TVector): Boolean; inline;
+    class operator Explicit(const aPoint: TVector): TPointF;
+    class operator :=(const aPoint: TPointF): TVector; inline;
+    class operator :=(const aPoint: TVector): TPointF; inline;
+    class operator :=(const aSize: TSizeF): TVector; inline;
+    class operator *(const aFactor: Single; const aVector: TVector): TVector; inline;
+    class operator *(const aVector: TVector; const aFactor: Single): TVector;
+    class operator <>(const aVector1, aVector2: TVector): Boolean; inline;
+    class operator -(const aVector1, aVector2: TVector): TVector;
+
+
+    function Length: Single;
+    function Normalize: TVector;
+    function CrossProduct(const aVector: TVector): TVector;
+    function DotProduct(const aVector: TVector): Single;
+    function MidVector(const aVector: TVector): TVector;
+
+    function ToPointF: TPointF; inline;
+    Function ToString: String;
+
+    case Integer of
+      0: (V: TVectorArray;);
+      1: (X: Single;
+          Y: Single;
+          W: Single;);
+  end;
+
+  { TPoint3D }
+
+  TPoint3D = record
+  Public
+    type
+      TPoint3DArray = array [0..2] of Single;
+    class function Create(const aX, aY, aZ: Single): TPoint3D; overload; static; inline;
+    class function Create(const P: TVector3DType): TPoint3D; overload; static; inline;
+    class function Create(const aPoint: TPointF; const aZ: Single = 0.0): TPoint3D; overload; static; inline;
+    class function Zero: TPoint3D; inline; static;
+
+    class operator +(const aPoint1, aPoint2: TPoint3D): TPoint3D;
+    class operator /(const aPoint: TPoint3D; const aFactor: Single): TPoint3D;
+    class operator =(const aPoint1, aPoint2: TPoint3D): Boolean; inline;
+    class operator *(const aFactor: Single; const aPoint: TPoint3D): TPoint3D; inline;
+    class operator *(const aPoint: TPoint3D; const aFactor: Single): TPoint3D; inline;
+    class operator *(const aPoint1, aPoint2: TPoint3D): TPoint3D;
+    class operator -(const aPoint: TPoint3D): TPoint3D;
+    class operator <>(const aPoint1, aPoint2: TPoint3D): Boolean; inline;
+    class operator -(const aPoint1, aPoint2: TPoint3D): TPoint3D;
+
+    function AngleCosine(const aPoint: TPoint3D): Single;
+    function CrossProduct(const aPoint: TPoint3D): TPoint3D;
+    function Distance(const aPoint: TPoint3D): Single;
+    function DotProduct(const aPoint: TPoint3D): Single;
+    function EqualsTo(const aPoint: TPoint3D; const aEpsilon: Single = 0): Boolean; inline;
+    function Length: Single; inline;
+    function MidPoint(const aPoint: TPoint3D): TPoint3D; inline;
+    function Normalize: TPoint3D;
+    procedure Offset(const aDelta: TPoint3D); overload; inline;
+    procedure Offset(const aDeltaX, aDeltaY, aDeltaZ: Single); overload; inline;
+    function Reflect(const aPoint: TPoint3D): TPoint3D; inline;
+    function Rotate(const aAxis: TPoint3D; const aAngle: Single): TPoint3D; inline;
+    Function ToString: String;
+  Public
+    case Integer of
+      0: (V: TPoint3DArray;);
+      1: (X: Single;
+          Y: Single;
+          Z: Single;);
+  end;
+  PPoint3D = ^TPoint3D;
+
+  tagVECTOR3D = record
+    case Integer of
+      0: (V: TVector3DType;);
+      1: (X: Single;
+          Y: Single;
+          Z: Single;
+          W: Single;);
+  end;
+
+
+  TMatrixArray = array [0..2] of TVector;
+  TMaxtrixArrayBase = array[0..2] of tagVECTOR;
+
+  { TMatrix }
+
+  TMatrix = record
+  public
+    class function CreateRotation(const aAngle: Single): TMatrix; static;
+    class function CreateScaling(const aScaleX, aScaleY: Single): TMatrix; static;
+    class function CreateTranslation(const aDeltaX, aDeltaY: Single): TMatrix; static;
+
+    class operator =(const RightMatrix, LeftMatrix: TMatrix): Boolean;
+    class operator *(const aMatrix1, aMatrix2: TMatrix): TMatrix;
+    class operator *(const aPoint: TPointF; const aMatrix: TMatrix): TPointF;
+    class operator *(const aVector: TVector; const aMatrix: TMatrix): TVector;
+    class operator *(const aVector: TPoint3D; const aMatrix: TMatrix): TPoint3D;
+    class operator *(const aFactor: Single; const aMatrix: TMatrix): TMatrix;
+    class operator *(const aMatrix: TMatrix; const aFactor: Single): TMatrix;
+    class operator /(const aMatrix: TMatrix; const aFactor: Single): TMatrix;
+
+    function Adjoint: TMatrix;
+    function Determinant: Single;
+    function EqualsTo(const aMatrix: TMatrix; const Epsilon: Single = TEpsilon.Matrix): Boolean;
+    function ExtractScale: TPointF;
+    function Inverse: TMatrix;
+    function Scale(const aFactor: Single): TMatrix;
+    Function ToString(Multiline : Boolean = true) : String;
+  Public
+    case Integer of
+      0: (M: TMatrixArray;);
+      1: (m11, m12, m13: Single;
+          m21, m22, m23: Single;
+          m31, m32, m33: Single);
+  end;
+
+  TMatrixConstants = record helper for TMatrix
+    const Identity: TMatrix = (m11: 1; m12: 0; m13: 0;
+                               m21: 0; m22: 1; m23: 0;
+                               m31: 0; m32: 0; m33: 1);
+  end;
+
+  PVector3D = ^TVector3D;
+
+  { TVector3D }
+
+  TVector3D = record
+  Public
+    class function Create(const aX, aY, aZ: Single; const aW: Single = 1.0): TVector3D; overload; static; inline;
+    class function Create(const aPoint: TPoint3D; const aW: Single = 1.0): TVector3D; overload; static; inline;
+    class function Zero: TVector3D; inline; static;
+
+    class operator +(const aVector1, aVector2: TVector3D): TVector3D;
+    class operator /(const aVector: TVector3D; const aFactor: Single): TVector3D;
+    class operator =(const aVector1, aVector2: TVector3D): Boolean;
+    class operator explicit(const aVector: TVector3D): TPoint3D;
+    class operator :=(const aPoint: TPoint3D): TVector3D;
+    class operator :=(const aVector: TVector3D): TPoint3D; inline;
+    class operator *(const aFactor: Single; const aVector: TVector3D): TVector3D; inline;
+    class operator *(const aVector: TVector3D; const aFactor: Single): TVector3D; inline;
+    class operator *(const aVector1, aVector2: TVector3D): TVector3D;
+    class operator -(const aVector: TVector3D): TVector3D;
+    class operator <>(const aVector1, aVector2: TVector3D): Boolean;
+    class operator -(const aVector1, aVector2: TVector3D): TVector3D;
+
+
+    function AngleCosine(const aVector: TVector3D): Single;
+    function CrossProduct(const aVector: TVector3D): TVector3D;
+    function Distance(const aVector: TVector3D): Single;
+    function DotProduct(const aVector: TVector3D): Single;
+    function EqualsTo(const aVector: TVector3D; const Epsilon: Single = 0): Boolean; inline;
+    function Length: Single;
+    function MidVector(const aVector: TVector3D): TVector3D;
+    function Normalize: TVector3D;
+    procedure Offset(const aDelta: TPoint3D); overload; inline;
+    procedure Offset(const aDeltaX, aDeltaY, aDeltaZ: Single); overload; inline;
+    function Reflect(const aVector: TVector3D): TVector3D; inline;
+    function Rotate(const aAxis: TPoint3D; const aAngle: Single): TVector3D; inline;
+    function ToPoint3D(const aTransform: Boolean = False): TPoint3D;
+  Public
+    case Integer of
+      0: (V: TVector3DType;);
+      1: (X: Single;
+          Y: Single;
+          Z: Single;
+          W: Single;);
+  end;
+
+  TVector3DArray = array [0..2] of TVector3D;
+  TVector3DArrayBase = array[0..2] of tagVECTOR3D;
+
+  TMatrix3DType = array [0..3] of TVector3D;
+  TMatrix3DTypeBase = array[0..3] of tagVECTOR3D;
+
+  { TMatrix3D }
+
+  TMatrix3D = record
+  public
+    constructor Create(const aM11, aM12, aM13, aM14,
+                             aM21, aM22, aM23, aM24,
+                             aM31, aM32, aM33, aM34,
+                             aM41, aM42, aM43, aM44: Single); overload;
+     { Transposed:
+       0 4 8  12
+       1 5 9  13
+       2 6 10 14
+       3 7 11 15 }
+    constructor Create(const aArray: TSingleDynArray); overload;
+
+    class Function Zero : TMatrix3D; static;
+
+    class function CreateLookAtDirLH(const aSource, aDirection, aCeiling: TPoint3D): TMatrix3D; static;
+    class function CreateLookAtDirRH(const aSource, aDirection, aCeiling: TPoint3D): TMatrix3D; static;
+    class function CreateLookAtLH(const aSource, aTarget, aCeiling: TPoint3D): TMatrix3D; static;
+    class function CreateLookAtRH(const aSource, aTarget, aCeiling: TPoint3D): TMatrix3D; static;
+
+    class function CreateOrthoLH(const aWidth, aHeight, aZNear, aZFar: Single): TMatrix3D; static;
+    class function CreateOrthoOffCenterLH(const aLeft, aTop, aRight, aBottom, aZNear, aZFar: Single): TMatrix3D; static;
+    class function CreateOrthoOffCenterRH(const aLeft, aTop, aRight, aBottom, aZNear, aZFar: Single): TMatrix3D; static;
+    class function CreateOrthoRH(const aWidth, aHeight, aZNear, aZFar: Single): TMatrix3D; static;
+
+    class function CreatePerspectiveFovLH(const aFOV, aAspect, aZNear, aZFar: Single; const aHorizontalFOV: Boolean = False): TMatrix3D; static;
+    class function CreatePerspectiveFovRH(const aFOV, aAspect, aZNear, aZFar: Single; const aHorizontalFOV: Boolean = False): TMatrix3D; static;
+
+    class function CreateRotation(const aAxis: TPoint3D; const aAngle: Single): TMatrix3D; static;
+    class function CreateRotationHeadingPitchBank(const aHeading, aPitch, aBank: Single): TMatrix3D; static;
+    class function CreateRotationX(const aAngle: Single): TMatrix3D; static;
+    class function CreateRotationY(const aAngle: Single): TMatrix3D; static;
+    class function CreateRotationYawPitchRoll(const aYaw, aPitch, aRoll: Single): TMatrix3D; static;
+    class function CreateRotationZ(const aAngle: Single): TMatrix3D; static;
+    class function CreateScaling(const aScale: TPoint3D): TMatrix3D; static;
+    class function CreateTranslation(const aTranslation: TPoint3D): TMatrix3D; static;
+
+    class operator *(const aPoint: TPoint3D; const aMatrix: TMatrix3D): TPoint3D;
+    class operator *(const aMatrix1, aMatrix2: TMatrix3D): TMatrix3D;
+    class operator *(const aVector: TVector3D; const aMatrix: TMatrix3D): TVector3D;
+    class operator *(const aFactor: single; const aMatrix: TMatrix3D) : TMatrix3D;
+    class operator *(const aMatrix: TMatrix3D; const aFactor: single) : TMatrix3D;
+    class operator /(const aMatrix: TMatrix3D; const aFactor: single) : TMatrix3D;
+
+    function Adjoint: TMatrix3D;
+    function Determinant: Single;
+    function EyePosition: TPoint3D;
+    function Inverse: TMatrix3D;
+    function Scale(const aFactor : Single) : TMatrix3D;
+    function ToMatrix: TMatrix;
+    function Transpose: TMatrix3D;
+    Function ToString(Multiline : Boolean = true) : String;
+
+    case Integer of
+      0: (M: TMatrix3DType;);
+      // m[row,col]
+      1: (m11, m12, m13, m14: Single;
+          m21, m22, m23, m24: Single;
+          m31, m32, m33, m34: Single;
+          m41, m42, m43, m44: Single);
+  end;
+
+  TMatrix3DConstants = record helper for TMatrix3D
+    const Identity: TMatrix3D = (m11: 1; m12: 0; m13: 0; m14: 0;
+                                 m21: 0; m22: 1; m23: 0; m24: 0;
+                                 m31: 0; m32: 0; m33: 1; m34: 0;
+                                 m41: 0; m42: 0; m43: 0; m44: 1;);
+  end;
+
+  { TQuaternion3D }
+
+  TQuaternion3D = record
+  Public
+    constructor Create(const aAxis: TPoint3D; const aAngle: Single); overload;
+    constructor Create(const P1,P2: TPoint3D); overload;
+    constructor Create(const aYaw, aPitch, aRoll: Single); overload;
+    constructor Create(const aMatrix: TMatrix3D); overload;
+
+    class operator :=(const aQuaternion: TQuaternion3D): TMatrix3D;
+    class operator *(const aQuaternion1, aQuaternion2: TQuaternion3D): TQuaternion3D;
+
+    function Length: Single;
+    function Normalize: TQuaternion3D;
+    function ToString : string;
+  Public
+    case Integer of
+      0: (V: TVector3DType;);
+      1: (ImagPart: TPoint3D;
+          RealPart: Single;);
+  end;
+
+  TQuaternion3DConstants = record helper for TQuaternion3D
+    const Identity: TQuaternion3D = (ImagPart: (X: 0; Y: 0; Z: 0);
+                                     RealPart: 1);
+  end;
+
+const
+  DefaultVectorWidth = Single(1.0);
+  NullVector3D: TVector3D = (X: 0; Y: 0; Z: 0; W: DefaultVectorWidth);
+  NullPoint3D: TPoint3D = (X: 0; Y: 0; Z: 0);
+
+function Vector(const X, Y: Single; const W: Single = DefaultVectorWidth): TVector; overload;
+function Vector(const P: TPointF; const W: Single = DefaultVectorWidth): TVector; overload;
+function Vector3D(const X, Y, Z: Single; const W: Single = DefaultVectorWidth): TVector3D; overload;
+function Vector3D(const P: TPoint3D; const W: Single = DefaultVectorWidth): TVector3D; overload;
+function Point3D(const X, Y, Z: Single): TPoint3D; overload;
+function Point3D(const aVector3D: TVector3D; const aTransform: Boolean = False): TPoint3D; overload;
+function PointF(const V: TVector): TPointF; inline; overload;
+
+implementation
+
+Function Display(D : Single) : String;
+
+begin
+  WriteStr(Result,D:7:4);
+end;
+
+{ ---------------------------------------------------------------------
+  TVector
+  ---------------------------------------------------------------------}
+
+class function TVector.Create(const aX, aY: Single; const aW: Single): TVector;
+
+begin
+  Result.X:=aX;
+  Result.Y:=aY;
+  Result.W:=aW;
+end;
+
+
+class function TVector.Create(const aPoint: TPointF): TVector;
+
+begin
+  Result.X:=aPoint.X;
+  Result.Y:=aPoint.Y;
+  Result.W:=DefaultVectorWidth;
+end;
+
+
+class function TVector.Zero: TVector;
+
+begin
+  Result.X:=0;
+  Result.Y:=0;
+  Result.W:=0;
+end;
+
+
+class operator TVector.+(const aVector1, aVector2: TVector): TVector;
+
+begin
+  Result.X:=aVector1.X+aVector2.X;
+  Result.Y:=aVector1.Y+aVector2.Y;
+  Result.W:=aVector1.W+aVector2.W;
+end;
+
+
+class operator TVector./(const aVector: TVector; const aFactor: Single): TVector;
+
+begin
+  Result.X:=aVector.X/aFactor;
+  Result.Y:=aVector.Y/aFactor;
+  Result.W:=aVector.W/aFactor;
+end;
+
+
+class operator TVector.=(const aVector1, aVector2: TVector): Boolean;
+
+begin
+  Result:=SameValue(aVector1.X,aVector2.X, TEpsilon.Vector)
+          and SameValue(aVector1.Y,aVector2.Y, TEpsilon.Vector)
+          and SameValue(aVector1.W,aVector2.W, TEpsilon.Vector);
+end;
+
+
+class operator TVector.Explicit(const aPoint: TVector): TPointF;
+
+var
+  S : Single;
+
+begin
+  S:=1;
+  if not SameValue(APoint.W,0,TEpsilon.Vector) then
+    S:=1/aPoint.W;
+  Result.X := APoint.X * S;
+  Result.Y := APoint.Y * S;
+end;
+
+class operator TVector.:=(const aPoint: TPointF): TVector;
+
+begin
+  Result:=TVector.Create(aPoint);
+end;
+
+
+class operator TVector.:=(const aPoint: TVector): TPointF;
+
+begin
+  Result:=TPointF(aPoint);
+end;
+
+class operator TVector.:=(const aSize: TSizeF): TVector;
+
+begin
+  Result:=TVector.Create(TPointF(aSize));
+end;
+
+
+class operator TVector.*(const aFactor: Single; const aVector: TVector): TVector;
+
+begin
+  Result.X:=aFactor*aVector.X;
+  Result.Y:=aFactor*aVector.Y;
+  Result.W:=aFactor*aVector.W;
+end;
+
+
+class operator TVector.*(const aVector: TVector; const aFactor: Single): TVector;
+begin
+  Result.X:=aVector.X*aFactor;
+  Result.Y:=aVector.Y*aFactor;
+  Result.W:=aVector.W*aFactor;
+end;
+
+
+class operator TVector.<>(const aVector1, aVector2: TVector): Boolean;
+
+begin
+  Result:=Not (aVector1=aVector2);
+end;
+
+
+class operator TVector.-(const aVector1, aVector2: TVector): TVector;
+
+begin
+  Result.X:=aVector1.X-aVector2.X;
+  Result.Y:=aVector1.Y-aVector2.Y;
+  Result.W:=aVector1.W-aVector2.W; // [Spock:] Highly illogical...
+end;
+
+
+function TVector.Length: Single;
+begin
+  Result:=Sqrt(DotProduct(Self));
+end;
+
+
+function TVector.Normalize: TVector;
+
+var
+  S: Single;
+
+begin
+  S:=Length;
+  if SameValue(S,0,TEpsilon.Vector) then
+    S:=1
+  else
+    S:=1/S;
+  Result:=S*Self;
+end;
+
+
+function TVector.CrossProduct(const aVector: TVector): TVector;
+
+begin
+  Result.X:=(Y*aVector.W)-(W*aVector.Y);
+  Result.Y:=(W*aVector.X)-(X*aVector.W);
+  Result.W:=(X*aVector.Y)-(Y*aVector.X);
+end;
+
+
+function TVector.DotProduct(const aVector: TVector): Single;
+
+begin
+  Result:=(X*aVector.X)+(Y*aVector.Y)+(W*aVector.W);
+end;
+
+
+function TVector.MidVector(const aVector: TVector): TVector;
+
+begin
+  Result:=(Self+aVector)/2;
+end;
+
+
+function TVector.ToPointF: TPointF;
+
+begin
+  Result:=TPointF(Self);
+end;
+
+
+Function TVector.ToString: String;
+
+begin
+  WriteStr(Result,'(',X:5:2,',',Y:5:2,' W:',W:5:2,')');
+end;
+
+
+{ ---------------------------------------------------------------------
+  TPoint3D
+  ---------------------------------------------------------------------}
+
+
+class function TPoint3D.Create(const aX, aY, aZ: Single): TPoint3D;
+
+begin
+  Result.X:=aX;
+  Result.Y:=aY;
+  Result.Z:=aZ;
+end;
+
+
+class function TPoint3D.Create(const P: TVector3DType): TPoint3D;
+
+begin
+  Result.X:=P[0];
+  Result.Y:=P[1];
+  Result.Z:=P[2];
+end;
+
+
+class function TPoint3D.Create(const aPoint: TPointF; const aZ: Single): TPoint3D;
+
+begin
+  Result.X:=aPoint.X;
+  Result.Y:=aPoint.Y;
+  Result.Z:=aZ;
+end;
+
+
+class function TPoint3D.Zero: TPoint3D;
+
+begin
+  Result.X:=0;
+  Result.Y:=0;
+  Result.Z:=0;
+end;
+
+
+class operator TPoint3D.+(const aPoint1, aPoint2: TPoint3D): TPoint3D;
+
+begin
+  Result.X:=aPoint1.X+aPoint2.X;
+  Result.Y:=aPoint1.Y+aPoint2.Y;
+  Result.Z:=aPoint1.Z+aPoint2.Z;
+end;
+
+
+class operator TPoint3D./(const aPoint: TPoint3D; const aFactor: Single): TPoint3D;
+
+begin
+  Result.X:=aPoint.X/aFactor;
+  Result.Y:=aPoint.Y/aFactor;
+  Result.Z:=aPoint.Z/aFactor;
+end;
+
+
+class operator TPoint3D.=(const aPoint1, aPoint2: TPoint3D): Boolean;
+
+begin
+  Result:=SameValue(aPoint1.X,aPoint2.X, TEpsilon.Vector)
+          and SameValue(aPoint1.Y,aPoint2.Y, TEpsilon.Vector)
+          and SameValue(aPoint1.Z,aPoint2.Z, TEpsilon.Vector);
+end;
+
+
+class operator TPoint3D.*(const aFactor: Single; const aPoint: TPoint3D): TPoint3D;
+
+begin
+  Result.X:=aFactor*aPoint.X;
+  Result.Y:=aFactor*aPoint.Y;
+  Result.Z:=aFactor*aPoint.Z;
+end;
+
+
+class operator TPoint3D.*(const aPoint: TPoint3D; const aFactor: Single): TPoint3D;
+
+begin
+  Result.X:=aPoint.X*aFactor;
+  Result.Y:=aPoint.Y*aFactor;
+  Result.Z:=aPoint.Z*aFactor;
+end;
+
+
+class operator TPoint3D.*(const aPoint1, aPoint2: TPoint3D): TPoint3D;
+
+begin
+  Result.X:=aPoint1.X*aPoint2.X;
+  Result.Y:=aPoint1.Y*aPoint2.Y;
+  Result.Z:=aPoint1.Z*aPoint2.Z;
+end;
+
+
+class operator TPoint3D.-(const aPoint: TPoint3D): TPoint3D;
+
+begin
+  Result.X:=-aPoint.X;
+  Result.Y:=-aPoint.Y;
+  Result.Z:=-aPoint.Z;
+end;
+
+
+class operator TPoint3D.<>(const aPoint1, aPoint2: TPoint3D): Boolean;
+
+begin
+  Result:=Not(aPoint1=aPoint2);
+end;
+
+
+class operator TPoint3D.-(const aPoint1, aPoint2: TPoint3D): TPoint3D;
+
+begin
+  Result.X:=aPoint1.X-aPoint2.X;
+  Result.Y:=aPoint1.Y-aPoint2.Y;
+  Result.Z:=aPoint1.Z-aPoint2.Z;
+end;
+
+
+function TPoint3D.DotProduct(const aPoint: TPoint3D): Single;
+
+begin
+  Result:=(X*aPoint.X)+(Y*aPoint.Y)+(Z*aPoint.Z);
+//  Writeln('Dot(',X,',',Y,',',Z,': ',Result);
+end;
+
+
+function TPoint3D.Length: Single;
+
+begin
+  Result:=Sqrt(DotProduct(Self));
+//  Writeln('Len:',Result);
+end;
+
+
+function TPoint3D.AngleCosine(const aPoint: TPoint3D): Single;
+
+var
+  Fact : Single;
+
+begin
+  Fact:=Length*APoint.Length;
+  if Abs(Fact)<=Epsilon then
+    Fact:=Epsilon;
+  Result:=DotProduct(APoint)/Fact;
+  Result:=Max(Result,-1);
+  Result:=Min(Result,1);
+end;
+
+
+function TPoint3D.CrossProduct(const aPoint: TPoint3D): TPoint3D;
+
+// https://en.wikipedia.org/wiki/Cross_product
+begin
+  Result.X:=(Y*aPoint.Z)-(Z*aPoint.Y);
+  Result.Y:=(Z*aPoint.X)-(X*aPoint.Z);
+  Result.Z:=(X*aPoint.Y)-(Y*aPoint.X);
+end;
+
+
+function TPoint3D.Distance(const aPoint: TPoint3D): Single;
+
+begin
+   Result:=(aPoint-Self).Length;
+end;
+
+
+function TPoint3D.EqualsTo(const aPoint: TPoint3D; const aEpsilon: Single): Boolean;
+
+begin
+  Result:=SameValue(X,aPoint.X,aEpsilon)
+          and SameValue(Y,aPoint.Y,aEpsilon)
+          and SameValue(Z,aPoint.Z,aEpsilon)
+end;
+
+
+function TPoint3D.MidPoint(const aPoint: TPoint3D): TPoint3D;
+
+begin
+  Result:=(aPoint+Self)/2;
+end;
+
+
+function TPoint3D.Normalize: TPoint3D;
+
+var
+  S : Single;
+
+begin
+  S:=Length;
+  if SameValue(S,0,TEpsilon.Vector) then
+    S:=1
+  else
+    S:=1/S;
+  Result:=Self*S;
+end;
+
+
+procedure TPoint3D.Offset(const aDelta: TPoint3D);
+
+begin
+  X:=X+aDelta.X;
+  Y:=Y+aDelta.Y;
+  Z:=Z+aDelta.Z;
+end;
+
+
+procedure TPoint3D.Offset(const aDeltaX, aDeltaY, aDeltaZ: Single);
+
+begin
+  X:=X+aDeltaX;
+  Y:=Y+aDeltaY;
+  Z:=Z+aDeltaZ;
+end;
+
+
+function TPoint3D.Reflect(const aPoint: TPoint3D): TPoint3D;
+
+begin
+  // Seems to be this, with aPoint the normal line:
+  // https://math.stackexchange.com/questions/13261/how-to-get-a-reflection-vector
+  Result:=Self-(aPoint*Self.DotProduct(APoint)*2);
+end;
+
+
+function TPoint3D.Rotate(const aAxis: TPoint3D; const aAngle: Single): TPoint3D;
+
+begin
+  Result:=Self*TMatrix3D.CreateRotation(aAxis,aAngle);
+end;
+
+
+function TPoint3D.ToString: String;
+
+begin
+  WriteStr(Result,'(',X:5:2,',',Y:5:2,',',Z:5:2,')');
+end;
+
+
+{ ---------------------------------------------------------------------
+  TMatrix
+  ---------------------------------------------------------------------}
+
+
+class function TMatrix.CreateRotation(const aAngle: Single): TMatrix;
+
+var
+  cA,sA : Single;
+
+begin
+  SinCos(aAngle,sA,cA);
+  Result:=TMatrix.Identity;
+  Result.m11 := cA;
+  Result.m12 := sA;
+  Result.m21 := -sA;
+  Result.m22 := cA;
+end;
+
+
+class function TMatrix.CreateScaling(const aScaleX, aScaleY: Single): TMatrix;
+
+begin
+  Result:=TMatrix.Identity;
+  Result.m11:=aScaleX;
+  Result.m22:=aScaleY;
+end;
+
+
+class function TMatrix.CreateTranslation(const aDeltaX, aDeltaY: Single): TMatrix;
+
+begin
+  Result:=TMatrix.Identity;
+  Result.m31:=aDeltaX;
+  Result.m32:=aDeltaY;
+end;
+
+
+class operator TMatrix.=(const RightMatrix, LeftMatrix: TMatrix): Boolean;
+
+begin
+  Result:=RightMatrix.EqualsTo(LeftMatrix,TEpsilon.Matrix);
+end;
+
+
+class operator TMatrix.*(const aMatrix1, aMatrix2: TMatrix): TMatrix;
+
+var
+  M1 : TMatrix absolute aMatrix1;
+  M2 : TMatrix absolute aMatrix2;
+  R : TMatrix;
+
+begin
+  With R do
+    begin
+    m11:=M1.m11*M2.m11 + M1.m12*M2.m21 + M1.m13*M2.m31;
+    m12:=M1.m11*M2.m12 + M1.m12*M2.m22 + M1.m13*M2.m32;
+    m13:=M1.m11*M2.m13 + M1.m12*M2.m23 + M1.m13*M2.m33;
+
+    m21:=M1.m21*M2.m11 + M1.m22*M2.m21 + M1.m23*M2.m31;
+    m22:=M1.m21*M2.m12 + M1.m22*M2.m22 + M1.m23*M2.m32;
+    m23:=M1.m21*M2.m13 + M1.m22*M2.m23 + M1.m23*M2.m33;
+
+    m31:=M1.m31*M2.m11 + M1.m32*M2.m21 + M1.m33*M2.m31;
+    m32:=M1.m31*M2.m12 + M1.m32*M2.m22 + M1.m33*M2.m32;
+    m33:=M1.m31*M2.m13 + M1.m32*M2.m23 + M1.m33*M2.m33;
+    end;
+  Result:=R;
+end;
+
+
+class operator TMatrix.*(const aPoint: TPointF; const aMatrix: TMatrix): TPointF;
+
+var
+  P : TPointF absolute aPoint;
+  N : TMatrix absolute aMatrix;
+
+begin
+  Result.X:=P.X*N.M11 + P.Y*N.M21 + N.M31;
+  Result.Y:=P.X*N.M12 + P.Y*N.M22 + N.M32
+end;
+
+
+class operator TMatrix.*(const aVector: TVector; const aMatrix: TMatrix): TVector;
+
+var
+  V : TVector absolute aVector;
+  N : TMatrix absolute aMatrix;
+
+begin
+  Result.X:=V.X*N.M11 + V.Y*N.M21 + V.W*N.M31;
+  Result.Y:=V.X*N.M12 + V.Y*N.M22 + V.W*N.M32;
+  Result.W:=V.X*N.M13 + V.Y*N.M23 + V.W*N.M33;
+end;
+
+
+class operator TMatrix.*(const aVector: TPoint3D; const aMatrix: TMatrix): TPoint3D;
+var
+  V : TPoint3D absolute aVector;
+  N : TMatrix absolute aMatrix;
+
+begin
+  Result.X:=V.X*N.M11 + V.Y*N.M21 + V.Z*N.M31;
+  Result.Y:=V.X*N.M12 + V.Y*N.M22 + V.Z*N.M32;
+  Result.Z:=V.X*N.M13 + V.Y*N.M23 + V.Z*N.M33;
+end;
+
+class operator TMatrix.*(const aFactor: Single; const aMatrix: TMatrix): TMatrix;
+
+begin
+  Result:=aMatrix.Scale(aFactor);
+end;
+
+
+class operator TMatrix.*(const aMatrix: TMatrix; const aFactor: Single): TMatrix;
+
+begin
+  Result:=aMatrix.Scale(aFactor);
+end;
+
+
+class operator TMatrix./(const aMatrix: TMatrix; const aFactor: Single): TMatrix;
+
+begin
+  Result:=aMatrix.Scale(1/aFactor);
+end;
+
+
+function TMatrix.Adjoint: TMatrix;
+// https://en.wikipedia.org/wiki/Adjugate_matrix
+var
+  a1, a2, a3, b1, b2, b3, c1, c2, c3: Single;
+  R : TMatrix;
+
+begin
+  a1:=Self.M11;
+  a2:=Self.M12;
+  a3:=Self.M13;
+  b1:=Self.M21;
+  b2:=Self.M22;
+  b3:=Self.M23;
+  c1:=Self.M31;
+  c2:=Self.M32;
+  c3:=Self.M33;
+  With R do
+    begin
+    M11:= (b2 * c3 - c2 * b3);
+    M12:=-(a2 * c3 - c2 * a3);
+    M13:= (a2 * b3 - b2 * a3);
+
+    M21:=-(b1 * c3 - c1 * b3);
+    M22:= (a1 * c3 - c1 * a3);
+    M23:=-(a1 * b3 - b1 * a3);
+
+    M31:= (b1 * c2 - c1 * b2);
+    M32:=-(a1 * c2 - c1 * a2);
+    M33:= (a1 * b2 - b1 * a2);
+    end;
+  Result:=R;
+end;
+
+
+// Determinant is calculated recursively.
+Function Det2x2(const a1, a2, b1, b2 : single) : Single; inline;
+
+begin
+  Result:=a1*b2-b1*a2;
+end;
+
+
+function Det3x3(const a1, a2, a3, b1, b2, b3, c1, c2, c3: Single): Single;
+
+begin
+  Result:=a1*Det2x2(b2,b3,c2,c3)-b1*Det2x2(a2,a3,c2,c3)+c1*Det2x2(a2,a3,b2,b3);
+end;
+
+
+function TMatrix.Determinant: Single;
+
+begin
+  Result:=Det3x3(m11,m12,m13,m21,m22,m23,m31,m32,m33);
+end;
+
+
+function TMatrix.EqualsTo(const aMatrix: TMatrix; const Epsilon: Single
+  ): Boolean;
+
+  Function Eq(s1,S2 : single) : Boolean; inline;
+  begin
+    Result:=SameValue(S1,S2,Epsilon);
+  end;
+
+Var
+  N : TMatrix absolute aMatrix;
+
+begin
+  Result:=Eq(m11,N.m11) and Eq(m12,N.m12) and Eq(m13,N.m13)
+          and Eq(m21,N.m21) and Eq(m22,N.m22) and Eq(m23,N.m23)
+          and Eq(m31,N.m31) and Eq(m32,N.m32) and Eq(m33,N.m33);
+end;
+
+function TMatrix.ExtractScale: TPointF;
+
+const
+  Origin : TPointF = (X:0;Y:0);
+
+begin
+  Result.X:=(PointF(1,0)*Self).Distance(Origin*Self);
+  Result.Y:=(PointF(0,1)*Self).Distance(Origin*Self);
+end;
+
+
+function TMatrix.Inverse: TMatrix;
+
+var
+  D: Single;
+
+begin
+  D:=Self.Determinant;
+  if SameValue(D,0,Epsilon) then
+    Result:=TMatrix.Identity
+  else
+    Result:=(Self.Adjoint/D);
+end;
+
+
+function TMatrix.Scale(const aFactor: Single): TMatrix;
+
+begin
+  Result.M11:=M11*aFactor;
+  Result.M12:=M12*aFactor;
+  Result.M13:=M13*aFactor;
+
+  Result.M21:=M21*aFactor;
+  Result.M22:=M22*aFactor;
+  Result.M23:=M23*aFactor;
+
+  Result.M31:=M31*aFactor;
+  Result.M32:=M32*aFactor;
+  Result.M33:=M33*aFactor;
+end;
+
+
+function TMatrix.ToString(Multiline: Boolean): String;
+
+var
+  S,Sep : String;
+
+begin
+  Sep:='';
+  if MultiLine then
+    Sep:=sLineBreak;
+  WriteStr(S,'[',Display(M11),',',Display(M12),',',Display(M13));
+  Result:=S+','+Sep;
+  WriteStr(S,Display(M21),',',Display(M22),',',Display(M23));
+  Result:=Result+S+Sep;
+  WriteStr(S,Display(M31),',',Display(M32),',',Display(M33));
+  Result:=Result+S+']';
+end;
+
+
+{ ---------------------------------------------------------------------
+  TVector3D
+  ---------------------------------------------------------------------}
+
+class function TVector3D.Create(const aX, aY, aZ: Single; const aW: Single): TVector3D;
+
+begin
+  Result.X:=aX;
+  Result.Y:=aY;
+  Result.Z:=aZ;
+  Result.W:=aW;
+end;
+
+
+class function TVector3D.Create(const aPoint: TPoint3D; const aW: Single): TVector3D;
+
+begin
+  Result.X:=aPoint.X;
+  Result.Y:=aPoint.Y;
+  Result.Z:=aPoint.Z;
+  Result.W:=aW;
+end;
+
+
+class function TVector3D.Zero: TVector3D;
+
+begin
+  Result.X:=0;
+  Result.Y:=0;
+  Result.Z:=0;
+  Result.W:=0;
+end;
+
+
+function TVector3D.EqualsTo(const aVector: TVector3D; const Epsilon: Single): Boolean;
+
+begin
+  Result:=SameValue(X,aVector.X,Epsilon)
+          and SameValue(Y,aVector.Y,Epsilon)
+          and SameValue(Z,aVector.Z,Epsilon)
+          and SameValue(W,aVector.W,Epsilon);
+end;
+
+
+class operator TVector3D.+(const aVector1, aVector2: TVector3D): TVector3D;
+
+var
+  V1 : TVector3D absolute aVector1;
+  V2 : TVector3D absolute aVector2;
+
+begin
+  Result.X:=V1.X+V2.X;
+  Result.Y:=V1.Y+V2.Y;
+  Result.Z:=V1.Z+V2.Z;
+  Result.W:=V1.W+V2.W;
+end;
+
+
+class operator TVector3D./(const aVector: TVector3D; const aFactor: Single): TVector3D;
+
+begin
+  With aVector do
+    begin
+    Result.X:=X/aFactor;
+    Result.Y:=Y/aFactor;
+    Result.Z:=Z/aFactor;
+    Result.W:=W/aFactor;
+    end;
+end;
+
+
+class operator TVector3D.=(const aVector1, aVector2: TVector3D): Boolean;
+
+begin
+  Result:=aVector1.EqualsTo(aVector2,TEpsilon.Vector);
+end;
+
+
+class operator TVector3D.explicit(const aVector: TVector3D): TPoint3D;
+
+begin
+  Result:=aVector.ToPoint3D;
+end;
+
+
+class operator TVector3D.:=(const aPoint: TPoint3D): TVector3D;
+
+begin
+  Result:=TVector3D.Create(APoint);
+end;
+
+
+class operator TVector3D.:=(const aVector: TVector3D): TPoint3D;
+
+begin
+  Result:=TPoint3D(AVector);
+end;
+
+
+class operator TVector3D.*(const aVector: TVector3D; const aFactor: Single): TVector3D;
+begin
+  With aVector do
+    begin
+    Result.X:=X*aFactor;
+    Result.Y:=Y*aFactor;
+    Result.Z:=Z*aFactor;
+    Result.W:=W*aFactor;
+    end;
+end;
+
+
+class operator TVector3D.*(const aFactor: Single; const aVector: TVector3D): TVector3D;
+
+begin
+  Result:=aVector*aFactor;
+end;
+
+
+class operator TVector3D.*(const aVector1, aVector2: TVector3D): TVector3D;
+
+var
+  V1 : TVector3D absolute aVector1;
+  V2 : TVector3D absolute aVector2;
+
+begin
+  Result.X:=V1.X*V2.X;
+  Result.Y:=V1.Y*V2.Y;
+  Result.Z:=V1.Z*V2.Z;
+  Result.W:=V1.W*V2.W;
+end;
+
+
+class operator TVector3D.-(const aVector: TVector3D): TVector3D;
+
+var
+  aV : TVector3D absolute aVector;
+
+begin
+  With Result do
+    begin
+    X:=-aV.X;
+    Y:=-aV.Y;
+    Z:=-aV.Z;
+    W:=-aV.W;
+    end;
+end;
+
+
+class operator TVector3D.<>(const aVector1, aVector2: TVector3D): Boolean;
+
+begin
+  Result:=Not (aVector1=aVector2)
+end;
+
+
+class operator TVector3D.-(const aVector1, aVector2: TVector3D): TVector3D;
+
+var
+  V1 : TVector3D absolute aVector1;
+  V2 : TVector3D absolute aVector2;
+
+begin
+  Result.X:=V1.X-V2.X;
+  Result.Y:=V1.Y-V2.Y;
+  Result.Z:=V1.Z-V2.Z;
+  Result.W:=V1.W-V2.W;
+end;
+
+
+function TVector3D.AngleCosine(const aVector: TVector3D): Single;
+
+var
+  P1,P2 : TPoint3D;
+
+begin
+  P1:=Self;
+  P2:=aVector;
+  Result:=P1.AngleCosine(P2);
+end;
+
+
+function TVector3D.CrossProduct(const aVector: TVector3D): TVector3D;
+
+var
+  aV : TVector3D absolute aVector;
+
+begin
+  Result.X:=Y*aV.Z-Z*aV.Y;
+  Result.Y:=Z*aV.X-X*aV.Z;
+  Result.Z:=X*aV.Y-Y*aV.X;
+  Result.W:=1;
+end;
+
+
+function TVector3D.Distance(const aVector: TVector3D): Single;
+
+begin
+  Result:=(AVector-Self).Length;
+end;
+
+
+function TVector3D.DotProduct(const aVector: TVector3D): Single;
+
+begin
+  Result:=X*aVector.X
+          + Y*aVector.Y
+          + Z*aVector.Z;
+end;
+
+
+function TVector3D.Length: Single;
+begin
+  Result:=Sqrt(Sqr(X)+Sqr(Y)+Sqr(Z)+Sqr(W));
+end;
+
+
+function TVector3D.MidVector(const aVector: TVector3D): TVector3D;
+begin
+  Result:=(Self+aVector)/2;
+end;
+
+
+function TVector3D.Normalize: TVector3D;
+
+var
+  Len: Single;
+begin
+  Len:=Self.Length;
+  if SameValue(0,Len,TEpsilon.Vector)  then
+    Result:=Self
+  else
+    Result:=Result/Len;
+end;
+
+
+function TVector3D.Reflect(const aVector: TVector3D): TVector3D;
+
+var
+  P1,P2 : TPoint3D;
+
+begin
+  P1:=Self;
+  P2:=aVector;
+  Result:=P1.Reflect(P2);
+end;
+
+
+function TVector3D.Rotate(const aAxis: TPoint3D; const aAngle: Single
+  ): TVector3D;
+begin
+  Result:=Self*TMatrix3D.CreateRotation(aAxis,aAngle);
+end;
+
+
+function TVector3D.ToPoint3D(const aTransform: Boolean): TPoint3D;
+
+var
+  Scale: Single;
+
+begin
+  Scale:=1;
+  if ATransform and not SameValue(W,0,TEpsilon.Vector) then
+    Scale:=1/W;
+  Result.X:=X*Scale;
+  Result.Y:=Y*Scale;
+  Result.Z:=Z*Scale;
+end;
+
+
+procedure TVector3D.Offset(const aDelta: TPoint3D);
+begin
+  X:=X+aDelta.X;
+  Y:=Y+aDelta.Y;
+  Z:=Z+aDelta.Z;
+end;
+
+
+procedure TVector3D.Offset(const aDeltaX, aDeltaY, aDeltaZ: Single);
+
+begin
+  X:=X+aDeltaX;
+  Y:=Y+aDeltaY;
+  Z:=Z+aDeltaZ;
+end;
+
+
+{ ---------------------------------------------------------------------
+  TMatrix3D
+  ---------------------------------------------------------------------}
+
+constructor TMatrix3D.Create(const aM11, aM12, aM13, aM14,
+                                   aM21, aM22, aM23, aM24,
+                                   aM31, aM32, aM33, aM34,
+                                   aM41, aM42, aM43, aM44: Single);
+
+begin
+  m11:=aM11;
+  m12:=aM12;
+  m13:=aM13;
+  m14:=aM14;
+  m21:=aM21;
+  m22:=aM22;
+  m23:=aM23;
+  m24:=aM24;
+  m31:=aM31;
+  m32:=aM32;
+  m33:=aM33;
+  m34:=aM34;
+  m41:=aM41;
+  m42:=aM42;
+  m43:=aM43;
+  m44:=aM44;
+end;
+
+
+constructor TMatrix3D.Create(const aArray: TSingleDynArray);
+begin
+ Self:=TMatrix3D.Create(aArray[0],aArray[4],aArray[8], aArray[12],
+                        aArray[1],aArray[5],aArray[9], aArray[13],
+                        aArray[2],aArray[6],aArray[10],aArray[14],
+                        aArray[3],aArray[7],aArray[11],aArray[15]);
+end;
+
+
+class function TMatrix3D.CreateLookAtDirLH(const aSource, aDirection, aCeiling: TPoint3D): TMatrix3D;
+
+// See e.g.
+// https://www.gamedev.net/tutorials/programming/graphics/perspective-projections-in-lh-and-rh-systems-r3598/
+// or See glscene https://github.com/glscene
+// GLS.VectorGeometry.pas : CreateLookAtMatrix
+
+var
+  Z,X,Y : TPoint3D;
+  R : TMatrix3D;
+
+begin
+  // Create reference axes X,Y,Z
+  Z:=-aDirection.Normalize;
+  X:=aCeiling.CrossProduct(Z).Normalize;
+  Y:=Z.CrossProduct(X); // Is normalized because X,Z normalized
+  R:=TMatrix3D.Identity;
+
+  // Translate
+  R.m41:= X.DotProduct(aSource);
+  R.m42:=-Y.DotProduct(aSource);
+  R.m43:=-Z.DotProduct(aSource);
+
+  // Rotate
+  R.m11 := X.X;
+  R.m12 := Y.X;
+  R.m13 := Z.X;
+
+  R.m21 := X.Y;
+  R.m22 := Y.Y;
+  R.m23 := Z.Y;
+
+  R.m31 := X.Z;
+  R.m32 := Y.Z;
+  R.m33 := Z.Z;
+
+  Result:=R;
+end;
+
+
+class function TMatrix3D.CreateLookAtDirRH(const aSource, aDirection,
+  aCeiling: TPoint3D): TMatrix3D;
+// See e.g.
+// https://www.gamedev.net/tutorials/programming/graphics/perspective-projections-in-lh-and-rh-systems-r3598/
+// or See glscene https://github.com/glscene
+// GLS.VectorGeometry.pas : CreateLookAtMatrix
+
+var
+  Z,X,Y : TPoint3D;
+  R : TMatrix3D;
+
+begin
+  // Create reference axes X,Y,Z
+  Z:=aDirection.Normalize;
+  X:=aCeiling.CrossProduct(Z).Normalize;
+  Y:=Z.CrossProduct(X); // Is normalized because X,Z normalized
+  R:=TMatrix3D.Identity;
+
+  // Translate
+  R.m41:=-X.DotProduct(aSource);
+  R.m42:=-Y.DotProduct(aSource);
+  R.m43:=-Z.DotProduct(aSource);
+
+  // Rotate
+  R.m11 := X.X;
+  R.m12 := Y.X;
+  R.m13 := Z.X;
+
+  R.m21 := X.Y;
+  R.m22 := Y.Y;
+  R.m23 := Z.Y;
+
+  R.m31 := X.Z;
+  R.m32 := Y.Z;
+  R.m33 := Z.Z;
+
+  Result:=R;
+end;
+
+
+class function TMatrix3D.CreateLookAtLH(const aSource, aTarget, aCeiling: TPoint3D): TMatrix3D;
+
+// See e.g.
+// https://www.gamedev.net/tutorials/programming/graphics/perspective-projections-in-lh-and-rh-systems-r3598/
+// or See glscene https://github.com/glscene
+// GLS.VectorGeometry.pas : CreateLookAtMatrix
+
+var
+  Z,X,Y : TPoint3D;
+  R : TMatrix3D;
+
+begin
+  // Create reference axes X,Y,Z
+  Z:=(aTarget-aSource).Normalize;
+  X:=aCeiling.CrossProduct(Z).Normalize;
+  Y:=Z.CrossProduct(X); // Is normalized because X,Z normalized
+  R:=TMatrix3D.Identity;
+
+  // Translate
+  R.m41:=-X.DotProduct(aSource);
+  R.m42:=-Y.DotProduct(aSource);
+  R.m43:=-Z.DotProduct(aSource);
+
+  // Rotate
+  // row 1
+  R.m11 := X.X;
+  R.m12 := Y.X;
+  R.m13 := Z.X;
+
+  // row 2
+  R.m21 := X.Y;
+  R.m22 := Y.Y;
+  R.m23 := Z.Y;
+
+  // row 3
+  R.m31 := X.Z;
+  R.m32 := Y.Z;
+  R.m33 := Z.Z;
+
+  Result:=R;
+end;
+
+
+class function TMatrix3D.CreateLookAtRH(const aSource, aTarget, aCeiling: TPoint3D): TMatrix3D;
+
+var
+  Z,X,Y : TPoint3D;
+  R : TMatrix3D;
+
+begin
+  // Create reference axes X,Y,Z
+  Z:=(aSource-aTarget).Normalize;
+  X:=aCeiling.CrossProduct(Z).Normalize;
+  Y:=Z.CrossProduct(X); // Is normalized because X,Z normalized
+  R:=TMatrix3D.Identity;
+
+  // Translate
+  R.m41:=-X.DotProduct(aSource);
+  R.m42:=-Y.DotProduct(aSource);
+  R.m43:=-Z.DotProduct(aSource);
+
+  // Rotate
+  // Row 1
+  R.m11 := X.X;
+  R.m12 := Y.X;
+  R.m13 := Z.X;
+
+  // Row 2
+  R.m21 := X.Y;
+  R.m22 := Y.Y;
+  R.m23 := Z.Y;
+
+  // Row 3
+  R.m31 := X.Z;
+  R.m32 := Y.Z;
+  R.m33 := Z.Z;
+
+  Result:=R;
+end;
+
+
+class function TMatrix3D.CreateOrthoLH(const aWidth, aHeight, aZNear, aZFar: Single): TMatrix3D;
+
+// See glscene https://github.com/glscene
+// GLS.VectorGeometry.pas, CreateMatrixFromFrustum and friends.
+
+var
+  R : TMatrix3D;
+  aLen : Single;
+
+begin
+  R:=TMatrix3D.Identity;
+  aLen:=(aZFar-aZNear);
+  R.m11:=2/aWidth;
+  R.m22:=2/aHeight;
+  R.m33:=1/aLen;
+  R.m42:=aZNear/(-alen);
+  Result:=R;
+end;
+
+
+class function TMatrix3D.CreateOrthoOffCenterLH(const aLeft, aTop, aRight, aBottom, aZNear, aZFar: Single): TMatrix3D;
+// See glscene https://github.com/glscene
+// GLS.VectorGeometry.pas, CreateMatrixFromFrustum and friends.
+
+var
+  R : TMatrix3D;
+  aWidth,aHeight,aLen : Single;
+
+begin
+  // Same as TMatrix3D.CreateOrthoLH, but offset
+  R:=TMatrix3D.Identity;
+  aWidth:=aRight-aleft;
+  aHeight:=aTop-aBottom;
+  aLen:=(aZFar-aZNear);
+  R.m11:=2/aWidth;
+  R.m22:=2/aHeight;
+  R.m33:=1/aLen;
+  R.m41:=(aLeft+aRight)/(-aWidth);
+  R.m42:=(aTop+aBottom)/(-aHeight);
+  R.m43:=aZNear/(-alen);
+  Result:=R;
+end;
+
+
+class function TMatrix3D.CreateOrthoOffCenterRH(const aLeft, aTop, aRight, aBottom, aZNear, aZFar: Single): TMatrix3D;
+// See glscene https://github.com/glscene
+// GLS.VectorGeometry.pas, CreateMatrixFromFrustum and friends.
+
+var
+  R : TMatrix3D;
+  aWidth,aHeight,aLen : Single;
+
+begin
+  // Same as TMatrix3D.CreateOrthoRH, but offset
+  R:=TMatrix3D.Identity;
+  aWidth:=aRight-aleft;
+  aHeight:=aTop-aBottom;
+  aLen:=(aZFar-aZNear);
+  R.m11:=2/aWidth;
+  R.m22:=2/aHeight;
+  R.m33:=1/-aLen;
+  R.m41:=(aLeft+aRight)/(-aWidth);
+  R.m42:=(aTop+aBottom)/(-aHeight);
+  R.m43:=aZNear/(-alen);
+  Result:=R;
+end;
+
+
+class function TMatrix3D.CreateOrthoRH(const aWidth, aHeight, aZNear, aZFar: Single): TMatrix3D;
+
+// See glscene https://github.com/glscene
+// GLS.VectorGeometry.pas, CreateMatrixFromFrustum and friends.
+
+var
+  R : TMatrix3D;
+
+begin
+  R:=TMatrix3D.Identity;
+  R.m11:=2/aWidth;
+  R.m22:=2/aHeight;
+  R.m33:=1/(aZNear-aZFar);
+  R.m42:=AZNear/(AZNear-AZFar);
+  Result:=R;
+end;
+
+
+class function TMatrix3D.CreatePerspectiveFovLH(const aFOV, aAspect, aZNear,
+  aZFar: Single; const aHorizontalFOV: Boolean): TMatrix3D;
+
+// see https//www.gamedev.net/tutorials/programming/graphics/perspective-projections-in-lh-and-rh-systems-r3598/
+// Perspective projection for A,B,C,D,E, DirectX calculation.
+
+var
+  aDist,sX,SY: Single;
+  R: TMatrix3D;
+
+begin
+  if aHorizontalFOV then
+    begin
+    sX:=1/Tan(aFOV/2);
+    sY:=sX/aAspect;
+    end
+  else
+    begin
+    sY:=1/Tan(aFOV/2);
+    sX:=sY/aAspect;
+    end;
+  aDist:=aZFar-aZNear;
+
+  R:=TMatrix3D.Identity;
+  R.m11:=sY; // A
+  R.m22:=sY; // B
+  R.m33:=aZFar/aDist; // C
+  R.m34:=1; // D
+  R.m43:=-aZNear*aZFar/aDist; // E
+  R.m44:=0;
+  Result:=R;
+end;
+
+
+class function TMatrix3D.CreatePerspectiveFovRH(const aFOV, aAspect, aZNear,
+  aZFar: Single; const aHorizontalFOV: Boolean): TMatrix3D;
+
+// see https//www.gamedev.net/tutorials/programming/graphics/perspective-projections-in-lh-and-rh-systems-r3598/
+// Perspective projection for A,B,C,D,E, DirectX calculation.
+
+var
+  aDist,sX,SY: Single;
+  R: TMatrix3D;
+
+begin
+  if aHorizontalFOV then
+    begin
+    sX:=1/Tan(aFOV/2);
+    sY:=sX/aAspect;
+    end
+  else
+    begin
+    sY:=1/Tan(aFOV/2);
+    sX:=sY/aAspect;
+    end;
+  aDist:=aZNear-aZFar;
+
+  R:=TMatrix3D.Identity;
+  R.m11:=sY; // A
+  R.m22:=sY; // B
+  R.m33:=aZFar/aDist; // C
+  R.m34:=-1; // D
+  R.m43:=aZNear*aZFar/aDist; // E
+  R.m44:=0;
+  Result:=R;
+end;
+
+
+class function TMatrix3D.CreateRotation(const aAxis: TPoint3D; const aAngle: Single): TMatrix3D;
+
+// See 'Rotation matrix from axis and angle' on
+// https://en.wikipedia.org/wiki/Rotation_matrix
+// Y-axis reversed
+
+var
+  U: TPoint3D;
+  C,S,Cos1: Single;
+  R : TMatrix3D;
+
+begin
+  SinCos(aAngle,S,C);
+  Cos1:=1-C;
+  U:=aAxis.Normalize;
+  R:=TMatrix3D.Identity;
+  With U do
+    begin
+    // Row 1
+    R.m11:=(X*X*Cos1) + C;
+    R.m12:=(X*Y*Cos1) + (Z*S);
+    R.m13:=(Z*X*Cos1) - (Y*S);
+
+    // Row 2
+    R.m21:=(X*Y*Cos1) - (Z*S);
+    R.m22:=(Y*Y*Cos1) + C;
+    R.m23:=(Y*Z*Cos1) + (X*S);
+
+    // Row 3
+    R.m31:=(Z*X*Cos1) + (Y*S);
+    R.m32:=(Y*Z*Cos1) - (X*S);
+    R.m33:=(Z*Z*Cos1) + C;
+    end;
+  Result:=R;
+end;
+
+
+class function TMatrix3D.CreateRotationHeadingPitchBank(const aHeading, aPitch,
+  aBank: Single): TMatrix3D;
+
+// Header-Pitch-Bank uses Tait-Brian angles
+// See https://en.wikipedia.org/wiki/Conversion_between_quaternions_and_Euler_angles
+// Y-axis reversed
+
+var
+  cHeading, sHeading, cPitch, sPitch, cBank, sBank : Single;
+  R : TMatrix3D;
+
+begin
+  SinCos(aHeading,sHeading,cHeading);  // Z axis (psi)
+  SinCos(aPitch,sPitch,cPitch);        // Y axis (Theta)
+  SinCos(aBank,sBank,cBank);           // X axis (phi)
+  R:=TMatrix3D.Identity;
+  // Row 1
+  R.m11:= cHeading*cBank + sHeading*sPitch*sBank;
+  R.m12:=-cHeading*sBank + sBank*sHeading*cBank;
+  R.m13:= sHeading*cPitch;
+
+  // Row 2
+  R.m21:=sBank*cPitch;
+  R.m22:=cBank*cPitch;
+  R.m23:=-sPitch;
+
+  // Row 3
+  R.m31:=-sHeading*cBank + cHeading*sPitch*sBank;
+  R.m32:=(sBank*sHeading) + cHeading*sPitch*cBank;
+  R.m33:=cHeading*cPitch;
+
+  Result:=R;
+end;
+
+
+class function TMatrix3D.CreateRotationX(const aAngle: Single): TMatrix3D;
+
+var
+  cA,sA: Single;
+  R: TMatrix3D;
+
+begin
+  SinCos(aAngle,sA,cA);
+  R:=TMatrix3D.Identity;
+  (*
+  [1 0   0 ]
+  [0 ca sa]
+  [0 -sa  ca]
+  *)
+  R.m11:=1;
+  R.m22:=cA;
+  R.m23:=sA;
+  R.m32:=-sA;
+  R.m33:=cA;
+  Result:=R;
+end;
+
+
+class function TMatrix3D.CreateRotationY(const aAngle: Single): TMatrix3D;
+var
+  cA,sA: Single;
+  R: TMatrix3D;
+
+begin
+  SinCos(aAngle,sA,cA);
+  R:=TMatrix3D.Identity;
+  (*
+  [ca 0  -sa]
+  [0  1   0 ]
+  [sa 0   ca]
+  *)
+  R.m11:=cA;
+  R.m13:=-sA;
+  R.m31:=sA;
+  R.m33:=cA;
+  Result:=R;
+end;
+
+
+class function TMatrix3D.CreateRotationYawPitchRoll(const aYaw, aPitch,
+  aRoll: Single): TMatrix3D;
+
+// See general 3D rotations in
+// https://en.wikipedia.org/wiki/Rotation_matrix
+// aYaw = Alpha = A, aPitch = Beta = B, aRoll = Gamma = G
+
+var
+  sYaw, cYaw, sPitch, cPitch, sRoll, cRoll: Single;
+  R : TMatrix3D;
+
+begin
+  SinCos(aYaw,SYaw,CYaw);
+  SinCos(aPitch,sPitch,CPitch);
+  SinCos(aRoll,SRoll,CRoll);
+  R:=TMatrix3D.Identity;
+  // Row 1
+  R.m11:=cRoll*cYaw        + sPitch*sRoll*sYaw;
+  R.m12:=cYaw*sPitch*sRoll - cRoll*sYaw;
+  R.m13:= -cPitch*sRoll;
+  // Row 2
+  R.m21:=cPitch*sYaw;
+  R.m22:=cPitch*cYaw;
+  R.m23:=sPitch;
+  // Row 3
+  R.m31:=cYaw*sRoll         - cRoll*sPitch*sYaw;
+  R.m32:=-cRoll*cYaw*sPitch - sRoll*sYaw;
+  R.m33:=cPitch*cRoll;
+
+  Result:=R;
+end;
+
+
+class function TMatrix3D.CreateRotationZ(const aAngle: Single): TMatrix3D;
+
+var
+  cA,sA: Single;
+  R : TMatrix3D;
+
+begin
+  SinCos(aAngle,sA,cA);
+  R:=TMatrix3D.Identity;
+  (*
+  [ca sa 0]
+  [-sa  ca 0]
+  [0   0  1]
+  *)
+
+  R.m11:=cA;
+  R.m12:=sA;
+  R.m21:=-sA;
+  R.m22:=cA;
+  Result:=R;
+end;
+
+
+class function TMatrix3D.CreateScaling(const aScale: TPoint3D): TMatrix3D;
+
+var
+  R : TMatrix3D;
+
+begin
+  R:=Zero;
+  R.m11:=aScale.X;
+  R.m22:=aScale.Y;
+  R.m33:=aScale.Z;
+  R.m44:=1;
+  Result:=R;
+end;
+
+
+class function TMatrix3D.CreateTranslation(const aTranslation: TPoint3D): TMatrix3D;
+
+var
+  R : TMatrix3D;
+
+begin
+  R:=TMatrix3D.Identity;
+  With aTranslation do
+    begin
+    R.m41:=X;
+    R.m42:=Y;
+    R.m43:=Z;
+    R.m44:=1;
+    end;
+  Result:=R;
+end;
+
+
+class operator TMatrix3D.*(const aPoint: TPoint3D; const aMatrix: TMatrix3D
+  ): TPoint3D;
+var
+  P : TPoint3D absolute aPoint;
+  N : TMatrix3D absolute aMatrix;
+begin
+  Result.X:=(P.X*N.m11)+(P.Y*N.m21)+(P.Z*N.m31)+N.m41;
+  Result.Y:=(P.X*N.m12)+(P.Y*N.m22)+(P.Z*N.m32)+N.m42;
+  Result.Z:=(P.X*N.m13)+(P.Y*N.m23)+(P.Z*N.m33)+N.m43;
+end;
+
+
+class operator TMatrix3D.*(const aMatrix1, aMatrix2: TMatrix3D): TMatrix3D;
+
+var
+  R,C,K : Integer;
+  S : Single;
+
+begin
+  For R:=0 to 3 do
+    For C:=0 to 3 do
+      begin
+      S:=0;
+      For K:=0 to 3 do
+        S:=S+(aMatrix1.M[R].V[K]*aMatrix2.M[K].V[C]);
+      Result.M[R].V[C]:=S;
+      end;
+end;
+
+
+class operator TMatrix3D.*(const aVector: TVector3D; const aMatrix: TMatrix3D
+  ): TVector3D;
+
+var
+  V : TVector3D absolute aVector;
+  N : TMatrix3D absolute aMatrix;
+
+begin
+  Result.X:=(V.X*N.m11)+(V.Y*N.m21)+(V.Z*N.m31)+(V.W*N.m41);
+  Result.Y:=(V.X*N.m12)+(V.Y*N.m22)+(V.Z*N.m32)+(V.W*N.m42);
+  Result.Z:=(V.X*N.m13)+(V.Y*N.m23)+(V.Z*N.m33)+(V.W*N.m43);
+  Result.W:=(V.X*N.m14)+(V.Y*N.m24)+(V.Z*N.m34)+(V.W*N.m44);
+end;
+
+
+class operator TMatrix3D.*(const aFactor: single; const aMatrix: TMatrix3D) : TMatrix3D;
+
+begin
+  Result:=aMatrix*aFactor;
+end;
+
+
+class operator TMatrix3D.*(const aMatrix: TMatrix3D; const aFactor: single) : TMatrix3D;
+
+var
+  R : TMatrix3D;
+  A : TMatrix3D absolute aMatrix;
+
+begin
+  R.M11:=A.M11*aFactor;
+  R.M12:=A.M12*aFactor;
+  R.M13:=A.M13*aFactor;
+  R.M14:=A.M14*aFactor;
+
+  R.M21:=A.M21*aFactor;
+  R.M22:=A.M22*aFactor;
+  R.M23:=A.M23*aFactor;
+  R.M24:=A.M24*aFactor;
+
+  R.M31:=A.M31*aFactor;
+  R.M32:=A.M32*aFactor;
+  R.M33:=A.M33*aFactor;
+  R.M34:=A.M34*aFactor;
+
+  R.M41:=A.M41*aFactor;
+  R.M42:=A.M42*aFactor;
+  R.M43:=A.M43*aFactor;
+  R.M44:=A.M44*aFactor;
+
+  Result:=R;
+end;
+
+
+class operator TMatrix3D./(const aMatrix: TMatrix3D; const aFactor: single
+  ): TMatrix3D;
+var
+  R : TMatrix3D;
+  A : TMatrix3D absolute aMatrix;
+
+begin
+  R.M11:=A.M11/aFactor;
+  R.M12:=A.M12/aFactor;
+  R.M13:=A.M13/aFactor;
+  R.M14:=A.M14/aFactor;
+
+  R.M21:=A.M21/aFactor;
+  R.M22:=A.M22/aFactor;
+  R.M23:=A.M23/aFactor;
+  R.M24:=A.M24/aFactor;
+
+  R.M31:=A.M31/aFactor;
+  R.M32:=A.M32/aFactor;
+  R.M33:=A.M33/aFactor;
+  R.M34:=A.M34/aFactor;
+
+  R.M41:=A.M41/aFactor;
+  R.M42:=A.M42/aFactor;
+  R.M43:=A.M43/aFactor;
+  R.M44:=A.M44/aFactor;
+
+  Result:=R;
+end;
+
+
+function TMatrix3D.Adjoint: TMatrix3D;
+
+var
+  a1, a2, a3, a4,
+  b1, b2, b3, b4,
+  c1, c2, c3, c4,
+  d1, d2, d3, d4: Single;
+  R : TMatrix3D;
+
+begin
+  a1 := M11;
+  b1 := M12;
+  c1 := M13;
+  d1 := M14;
+
+  a2 := M21;
+  b2 := M22;
+  c2 := M23;
+  d2 := M24;
+
+  a3 := M31;
+  b3 := M32;
+  c3 := M33;
+  d3 := M34;
+
+  a4 := M41;
+  b4 := M42;
+  c4 := M43;
+  d4 := M44;
+
+  With R do
+    begin
+    M11:= Det3x3(b2,b3,b4, c2,c3,c4, d2,d3,d4);
+    M21:=-Det3x3(a2,a3,a4, c2,c3,c4, d2,d3,d4);
+    M31:= Det3x3(a2,a3,a4, b2,b3,b4, d2,d3,d4);
+    M41:=-Det3x3(a2,a3,a4, b2,b3,b4, c2,c3,c4);
+
+    M12:=-Det3x3(b1,b3,b4, c1,c3,c4, d1,d3,d4);
+    M22:= Det3x3(a1,a3,a4, c1,c3,c4, d1,d3,d4);
+    M32:=-Det3x3(a1,a3,a4, b1,b3,b4, d1,d3,d4);
+    M42:= Det3x3(a1,a3,a4, b1,b3,b4, c1,c3,c4);
+
+    M13:= Det3x3(b1,b2,b4, c1,c2,c4, d1,d2,d4);
+    M23:=-Det3x3(a1,a2,a4, c1,c2,c4, d1,d2,d4);
+    M33:= Det3x3(a1,a2,a4, b1,b2,b4, d1,d2,d4);
+    M43:=-Det3x3(a1,a2,a4, b1,b2,b4, c1,c2,c4);
+
+    M14:=-Det3x3(b1,b2,b3, c1,c2,c3, d1,d2,d3);
+    M24:= Det3x3(a1,a2,a3, c1,c2,c3, d1,d2,d3);
+    M34:=-Det3x3(a1,a2,a3, b1,b2,b3, d1,d2,d3);
+    M44:= Det3x3(a1,a2,a3, b1,b2,b3, c1,c2,c3);
+    end;
+  Result:=R;
+end;
+
+
+function TMatrix3D.Determinant: Single;
+
+begin
+  Result:= (M11*Det3x3(m22,m23,m24, m32,m33,m34, m42,m43,m44))
+          -(M12*Det3x3(m21,m23,m24, m31,m33,m34, m41,m43,m44))
+          +(M13*Det3x3(m21,m22,m24, m31,m32,m34, m41,m42,m44))
+          -(M14*Det3x3(m21,m22,m23, m31,m32,m33, m41,m42,m43))
+end;
+
+
+function TMatrix3D.EyePosition: TPoint3D;
+
+begin
+  Result.X:=-(M11*M41 + M12*M42 + M13*M43);
+  Result.Y:=-(M21*M41 + M22*M42 + M23*M43);
+  Result.Z:=-(M31*M41 + M32*M42 + M33*M43);
+end;
+
+
+function TMatrix3D.Inverse: TMatrix3D;
+
+var
+  D: Single;
+
+begin
+  D:=Self.Determinant;
+  if SameValue(D,0,Epsilon) then
+    Result:=TMatrix3D.Identity
+  else
+    Result:=(Self.Adjoint/D);
+end;
+
+function TMatrix3D.Scale(const aFactor: Single): TMatrix3D;
+begin
+  Result:=Self*aFactor;
+end;
+
+
+function TMatrix3D.ToMatrix: TMatrix;
+
+begin
+  Result.m11:=m11;
+  Result.m12:=m12;
+  Result.m13:=m13;
+  Result.m21:=m21;
+  Result.m22:=m22;
+  Result.m23:=m23;
+  Result.m31:=m31;
+  Result.m32:=m32;
+  Result.m33:=m33;
+end;
+
+
+function TMatrix3D.Transpose: TMatrix3D;
+
+begin
+  Result.M11:=M11;
+  Result.M21:=M12;
+  Result.M31:=M13;
+  Result.M41:=M14;
+
+  Result.M12:=M21;
+  Result.M22:=M22;
+  Result.M32:=M23;
+  Result.M42:=M24;
+
+  Result.M13:=M31;
+  Result.M23:=M32;
+  Result.M33:=M33;
+  Result.M43:=M34;
+
+  Result.M14:=M41;
+  Result.M24:=M42;
+  Result.M34:=M43;
+  Result.M44:=M44;
+end;
+
+function TMatrix3D.ToString(Multiline: Boolean): String;
+
+var
+  S,Sep : String;
+
+begin
+  Sep:='';
+  if MultiLine then
+    Sep:=sLineBreak;
+  WriteStr(S,Display(m11),',',Display(M12),',',Display(M13),',',Display(M14));
+  Result:=S+','+Sep;
+  WriteStr(S,Display(M21),',',Display(M22),',',Display(M23),',',Display(M24));
+  Result:=Result+S+','+Sep;
+  WriteStr(S,Display(M31),',',Display(M32),',',Display(M33),',',Display(M34));
+  Result:=Result+S+','+Sep;
+  WriteStr(S,Display(M41),',',Display(M42),',',Display(M43),',',Display(M44));
+  Result:='['+Result+S+']';
+end;
+
+
+class function TMatrix3D.Zero: TMatrix3D; static;
+
+begin
+  Result:=Default(TMatrix3D);
+end;
+
+
+{ ---------------------------------------------------------------------
+  TQuaternion3D
+  ---------------------------------------------------------------------}
+
+constructor TQuaternion3D.Create(const aAxis: TPoint3D; const aAngle: Single);
+
+var
+  L,sA,cA : Single;
+
+begin
+  L:=aAxis.Length;
+  if SameValue(L,0,Epsilon2) then
+    Self:=TQuaternion3D.Identity
+  else
+    begin
+    SinCos(aAngle/2,sA,cA);
+    Self.RealPart:=cA;
+    Self.ImagPart:=aAxis*(sA/L);
+    end;
+end;
+
+constructor TQuaternion3D.Create(const P1, P2: TPoint3D);
+begin
+  ImagPart:=P1.CrossProduct(P2);
+  RealPart:=Sqrt((P1.DotProduct(P2)+1)/2)
+end;
+
+
+constructor TQuaternion3D.Create(const aYaw, aPitch, aRoll: Single);
+
+begin
+  Self:=TQuaternion3D.Create(Point3D(0,1,0),AYaw) // rotation around Y
+        *TQuaternion3D.Create(Point3D(1,0,0),APitch) // rotation around X
+        *TQuaternion3D.Create(Point3D(0,0,1),ARoll); // rotation around Z
+end;
+
+constructor TQuaternion3D.Create(const aMatrix: TMatrix3D);
+
+// See glscene  QuaternionFromMatrix
+
+var
+  Q : TQuaternion3D;
+  M : TMatrix3D absolute aMatrix;
+  Diag, S : Double;
+
+begin
+  Diag:=1+m.m11+m.m22+m.m33;
+  if Diag>Epsilon then
+    begin
+    // All large
+    S:=2*Sqrt(Diag);
+    Q.RealPart:=S/4;
+    S:=1/S;
+    Q.ImagPart.X:=(M.m23 - M.m32)*S;
+    Q.ImagPart.Y:=(M.m31 - M.m13)*S;
+    Q.ImagPart.Z:=(M.m12 - M.m21)*S;
+    end
+  else if (M.m11>M.m22) and (M.m11>M.m33) then
+    begin
+    // XX is largest
+    S:=2*Sqrt(Max(Epsilon, (1 + M.m11-M.m22-M.m33)));
+    Q.ImagPart.X:=(S/4);
+    S:=1/S;
+    Q.ImagPart.Y:=(M.m12 + M.m21)*S;
+    Q.ImagPart.Z:=(M.m31 + M.m13)*S;
+    Q.RealPart  :=(M.m23 - M.m32)*S;
+    end
+  else if (M.m22 > M.m33) then
+    begin
+    // YY is largest
+    S:=2*Sqrt(Max(Epsilon, (1 + M.m22-M.m11-M.m33)));
+    Q.ImagPart.Y:=S/4;
+    S:=1/S;
+    Q.ImagPart.X:=(M.m12 + M.m21)*S;
+    Q.ImagPart.Z:=(M.m23 + M.m32)*S;
+    Q.RealPart  :=(M.m31 - M.m13)*S;
+    end
+  else
+    begin
+    // ZZ is largest
+    S:=2*Sqrt(Max(Epsilon,(1+M.m33-M.m11-M.m22)));
+    Q.ImagPart.Z:=S/4;
+    S:=1/S;
+    Q.ImagPart.X:=(M.m31 + M.m13)*S;
+    Q.ImagPart.Y:=(M.m23 + M.m32)*S;
+    Q.RealPart  :=(M.m12 - M.m21)*S;
+    end;
+  Self:=Q.Normalize;
+end;
+
+
+class operator TQuaternion3D.:=(const aQuaternion: TQuaternion3D): TMatrix3D;
+
+var
+  N: TQuaternion3D;
+  xx, xy, xz, xw, yy, yz, yw, zz, zw: Single;
+  R : TMatrix3D;
+
+begin
+  N:=AQuaternion.Normalize;
+  With N do
+    begin
+    xx:=Sqr(ImagPart.X);
+    xy:=ImagPart.X*ImagPart.Y;
+    xz:=ImagPart.X*ImagPart.Z;
+    xw:=ImagPart.X*RealPart;
+
+    yy:=Sqr(ImagPart.Y);
+    yz:=ImagPart.Y*ImagPart.Z;
+    yw:=ImagPart.Y*RealPart;
+
+    zz:=Sqr(ImagPart.Z);
+    zw:=ImagPart.Z*RealPart;
+    end;
+  R:=TMatrix3D.Zero;
+  With R do
+    begin
+    M11:=1-2*(yy+zz);
+    M21:=2*(xy-zw);
+    M31:=2*(xz+yw);
+
+    M12:=2*(xy+zw);
+    M22:=1-2*(xx+zz);
+    M32:=2*(yz-xw);
+
+    M13:=2*(xz-yw);
+    M23:=2*(yz+xw);
+    M33:=1-2*(xx+yy);
+
+    M44:=1;
+    end;
+  Result:=R;
+end;
+
+
+class operator TQuaternion3D.*(const aQuaternion1, aQuaternion2: TQuaternion3D): TQuaternion3D;
+
+var
+  Q1 : TQuaternion3D absolute aQuaternion1;
+  Q2 : TQuaternion3D absolute aQuaternion2;
+
+begin
+  Result.RealPart:= Q1.RealPart*Q2.RealPart - Q1.ImagPart.DotProduct(Q2.ImagPart);
+  Result.ImagPart.X:=Q1.RealPart*Q2.ImagPart.X
+                     + Q2.RealPart*Q1.ImagPart.X
+                     + Q1.ImagPart.Y*Q2.ImagPart.Z
+                     - Q1.ImagPart.Z*Q2.ImagPart.Y;
+  Result.ImagPart.Y:=Q1.RealPart*Q2.ImagPart.Y
+                     + Q2.RealPart*Q1.ImagPart.Y
+                     + Q1.ImagPart.Z*Q2.ImagPart.X
+                     - Q1.ImagPart.X*Q2.ImagPart.Z;
+  Result.ImagPart.Z:=Q1.RealPart*Q2.ImagPart.Z
+                     + Q2.RealPart*Q1.ImagPart.Z
+                     + Q1.ImagPart.X*Q2.ImagPart.Y
+                     - Q1.ImagPart.Y*Q2.ImagPart.X;
+end;
+
+
+function TQuaternion3D.Length: Single;
+
+begin
+  Result:=Sqrt(Sqr(ImagPart.X)+Sqr(ImagPart.Y)+Sqr(ImagPart.Z)+Sqr(RealPart));
+end;
+
+
+function TQuaternion3D.Normalize: TQuaternion3D;
+
+var
+  Len : Single;
+
+begin
+  Len:=Length;
+  If SameValue(Len,0,Epsilon) then
+    Result:=TQuaternion3D.Identity
+  else
+    begin
+    Len:=1/Len;
+    Result.ImagPart:=ImagPart*Len;
+    Result.RealPart:=RealPart*Len;
+    end;
+end;
+
+function TQuaternion3D.ToString: string;
+begin
+  WriteStr(Result,'(',Display(RealPart),',i ',Display(ImagPart.X),',j ',Display(ImagPart.Y),',k ',Display(ImagPart.Z),')');
+end;
+
+{ ---------------------------------------------------------------------
+  Auxiliary functions
+  ---------------------------------------------------------------------}
+
+
+function Vector(const X, Y: Single; const W: Single = 1.0): TVector;
+begin
+  Result.X:=X;
+  Result.Y:=Y;
+  Result.W:=W;
+end;
+
+
+function Vector(const P: TPointF; const W: Single = 1.0): TVector;
+begin
+  Result.X:=P.X;
+  Result.Y:=P.Y;
+  Result.W:=W;
+end;
+
+function Vector3D(const X, Y, Z: Single; const W: Single = 1.0): TVector3D;
+
+begin
+  Result:=TVector3D.Create(X,Y,Z,W);
+end;
+
+
+function Vector3D(const P: TPoint3D; const W: Single = 1.0): TVector3D;
+
+begin
+  Result:=TVector3D.Create(P,W);
+end;
+
+
+function Point3D(const X, Y, Z: Single): TPoint3D;
+
+begin
+  Result:=TPoint3D.Create(X,Y,Z);
+end;
+
+
+function Point3D(const aVector3D: TVector3D; const aTransform: Boolean): TPoint3D;
+
+begin
+  Result:=AVector3D.ToPoint3D(aTransform);
+end;
+
+
+function PointF(const V: TVector): TPointF;
+
+begin
+  Result:=TPointF(V);
+end;
+
+end.

packages/rtl-objpas/tests/testrunner.rtlobjpas.pp

@@ -33,7 +33,14 @@ uses
 {$ifdef testimpl}
   tests.rtti.impl,
 {$endif}
-  tests.rtti, tests.rtti.value, tests.rtti.types;
+  tests.rtti,
+  tests.rtti.value,
+  tests.rtti.types,
+  utmathvectorbase,
+  utcmatrix,
+  utcpoint,
+  utcvector,
+  utcquaternion;
 
 var
   Application: TTestRunner;

packages/rtl-objpas/tests/tv.dpr

@@ -0,0 +1,562 @@
+{
+    This file is part of the Free Pascal run time library.
+    Copyright (c) 2023 by Michael Van Canneyt
+    member of the Free Pascal development team
+
+    See the file COPYING.FPC, included in this distribution,
+    for details about the copyright.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+ **********************************************************************}
+
+{ 
+  Tests were made to determine where delphi implementation differs 
+  from standards and glscene,  and adapt the implementation
+  The numbers of tests done here can also be found in the unit tests.
+  
+  run this application with number of test to run only that test.
+}
+
+{$IFDEF FPC}
+{$mode delphi}
+uses SysUtils, System.Math.Vectors;
+{$ELSE}
+uses System.SysUtils, System.Math.Vectors;
+{$ENDIF}
+
+{$IFNDEF FPC}
+Function MIdentity : TMatrix;
+
+begin
+  Result:=TMatrix.Identity;
+end;
+
+Function M3DIdentity : TMatrix3D;
+
+begin
+  Result:=TMatrix3D.Identity;
+end;
+
+// Delphi does not have ToString, use helper to implement it.
+// We cannot inherit record helpers in Delphi, 
+// so we define this after the above functions to have Identity.  
+Type
+  TVHelper = record Helper for TVector 
+    Function ToString: String;
+  end;
+
+  TPHelper = record Helper for TPoint3D 
+    Function ToString: String;
+  end;
+
+  TMHelper = record Helper {(TMatrixConstants)} for TMatrix 
+    function ToString(Multiline: Boolean = True): String; 
+    class function Identity : TMatrix; static;
+  end;
+
+  TM3DHelper = record Helper {(TMatrixConstants)} for TMatrix3D 
+    function ToString(Multiline: Boolean = True): String; 
+    class function Identity : TMatrix3D; static;
+  end;
+
+  TQ3DHelper = record Helper {(TMatrixConstants)} for TQuaternion3D 
+    function ToString(Multiline: Boolean = True): String; 
+  end;
+
+Function TVHelper.ToString: String;
+begin
+  Result:=Format('(%7.4f,%7.4f W:%7.4f)',[X,Y,W]);
+end;
+
+Function TPHelper.ToString: String;
+begin
+  Result:=Format('(%7.4f,%7.4f,%7.4f)',[X,Y,Z]);
+end;
+
+
+class function TMHelper.Identity : TMatrix;
+begin
+  Result:=Midentity;
+end;
+
+function TMHelper.ToString(Multiline: Boolean): String;
+var
+  S,Sep : String;
+
+begin
+  Sep:='';
+  if MultiLine then
+    Sep:=sLineBreak;
+  S:='['+Format('%7.4f,%7.4f,%7.4f',[m11,m12,m13]);
+  Result:=S+','+Sep;
+  S:=Format('%7.4f,%7.4f,%7.4f',[m21,m22,m23]);
+  Result:=Result+S+Sep;
+  S:=Format('%7.4f,%7.4f,%7.4f',[m31,m32,m33]);
+  Result:=Result+S+']';
+end; 
+
+class function TM3DHelper.Identity : TMatrix3D;
+begin
+  Result:=M3Didentity;
+end;
+
+function TM3DHelper.ToString(Multiline: Boolean): String;
+var
+  S,Sep : String;
+
+begin
+  Sep:='';
+  if MultiLine then
+    Sep:=sLineBreak;
+  S:='['+Format('%7.4f,%7.4f,%7.4f,%7.4f',[m11,m12,m13,m14]);
+  Result:=S+','+Sep;
+  S:=Format('%7.4f,%7.4f,%7.4f,%7.4f',[m21,m22,m23,m24]);
+  Result:=Result+S+','+Sep;
+  S:=Format('%7.4f,%7.4f,%7.4f,%7.4f',[m31,m32,m33,m34]);
+  Result:=Result+S+','+Sep;
+  S:=Format('%7.4f,%7.4f,%7.4f,%7.4f',[m41,m42,m43,m44]);
+  Result:=Result+S+']';
+end; 
+
+ 
+ function TQ3DHelper.ToString(Multiline: Boolean = True): String; 
+ 
+ begin
+   Result:=Format('(%7.4f,i: %7.4f,j: %7.4f, k: %7.4f)',[Realpart,ImagPart.X,ImagPart.Y,ImagPart.Z]);
+ end;
+
+{$ENDIF}
+
+Function DoTest(aId : integer; const aCaption : String) : boolean;
+
+var
+  cID : Integer;
+
+begin
+  cID:=StrToIntDef(Paramstr(1),-1);
+  Result:=(aId=cID) or (Cid=-1);
+  if Result then
+    Writeln('[',aID,'] ',aCaption);
+end;
+
+Procedure TV;
+
+var
+  V1,V2,V3 : TVector;
+
+begin
+  if DoTest(1,'Vector CrossProduct') then
+    begin
+    V1:=TVector.Create(1,0,0);
+    V2:=TVector.Create(0,1,0);
+    V3:=V2.CrossProduct(V1);
+    With V3 do
+      Writeln('1: V3 ',V3.ToString);
+    end;
+end;
+
+Procedure TP;
+
+var
+  P1,P2,P3 : TPoint3D;
+
+begin
+  P1:=TPoint3D.Create(1,0,0);
+  P2:=TPoint3D.Create(0,0,1); // Z-axis
+  if DoTest(101,'Rotate Z axis') then
+    begin
+    P3:=P1.Rotate(P2,Pi/2);
+    Writeln('P3 = ',P3.ToString);
+    end;
+  if DoTest(102,'Rotate Y axis') then
+    begin
+    P2:=TPoint3D.Create(0,1,0); // Y-axis
+    P3:=P1.Rotate(P2,Pi/2);
+    Writeln('P3 = ',P3.ToString);
+    end;
+end;
+
+Procedure TM;
+
+var
+  M,M2 : TMatrix;
+
+begin
+  if DoTest(11,'Identity') then
+    begin
+    M:=TMatrix.Identity;
+    Writeln(M.ToString);
+    end;
+  If DoTest(12,'Rotation (Pi/2)') then
+    begin
+    M:=TMatrix.CreateRotation(Pi/2);
+    Writeln(M.ToString);
+    end;
+  if DoTest(13,'Rotation Pi') then
+    begin
+    M:=TMatrix.CreateRotation(Pi);
+    Writeln(M.ToString);
+    end;
+  If DoTest(14,'Scaling 2') then
+    begin
+    M:=TMatrix.CreateScaling(2,3);
+    Writeln(M.ToString);
+    end;
+  If DoTest(15,'Translation (3,4)') then
+    begin
+    M:=TMatrix.CreateTranslation(3,4);
+    Writeln(M.ToString);
+    end;
+  If DoTest(16,'Adjoint') then
+    begin
+    M:=TMatrix.CreateRotation(Pi/6);
+    M2:=M.Adjoint;
+    Writeln(M2.ToString);
+    end;
+end;
+
+Procedure Params(const ASource, ATarget, ACeiling: TPoint3D);
+var
+  ZAxis, XAxis, YAxis: TPoint3D;
+begin
+  Writeln('aSource: ',aSource.ToString);
+  Writeln('aTarget: ',aTarget.ToString);
+  Writeln('aCeiling: ',aCeiling.ToString);
+ 
+  ZAxis := (ASource - ATarget).Normalize;
+  Writeln('Zaxis: ',ZAxis.ToString);
+  XAxis := ACeiling.CrossProduct(ZAxis).Normalize;
+  Writeln('Xaxis: ',XAxis.ToString);
+  YAxis := ZAxis.CrossProduct(XAxis);
+  Writeln('Yaxis: ',YAxis.ToString);
+end;
+
+Procedure TM3D;
+
+var
+  M,M2,M3 : TMatrix3D;
+
+begin
+  if DoTest(31,'M3D.Identity') then
+    begin
+    M:=TMatrix3D.Identity;
+    Writeln(M.ToString);
+    end;
+  if DoTest(32,'M3D.Rotation ((1,0,0),Pi/2)') then
+    begin
+    M:=TMatrix3D.CreateRotation(Point3D(1,0,0),Pi/2);
+    Writeln(M.ToString);
+    end;
+  if DoTest(33,'M3D.Rotation ((0,1,0),Pi/2)') then
+    begin
+    M:=TMatrix3D.CreateRotation(Point3D(0,1,0),Pi/2);
+    Writeln(M.ToString);
+    end;
+  if DoTest(34,'M3D.Rotation ((0,0,1),Pi/2)') then
+    begin
+    M:=TMatrix3D.CreateRotation(Point3D(0,0,1),pi/2);
+    Writeln(M.ToString);
+    end;
+  if DoTest(35,'M3D.Scaling (1,2,3)') then
+    begin
+    M:=TMatrix3D.CreateScaling(Point3D(1,2,3));
+    Writeln(M.ToString);
+    end;
+  if DoTest(36,'MD3D.Translation (3,4,5)') then
+    begin
+    M:=TMatrix3d.CreateTranslation(Point3D(3,4,5));
+    Writeln(M.ToString);
+    end;
+  if DoTest(37,'M3D.HeadingPitchBank(Pi/2,0,0)') then
+    begin
+    M:=TMatrix3D.CreateRotationHeadingPitchBank(Pi/2,0,0);
+    Writeln(M.ToString);
+    end;
+  if DoTest(38,'M3D.HeadingPitchBank(0,Pi/2,0)') then
+    begin
+    M:=TMatrix3D.CreateRotationHeadingPitchBank(0,Pi/2,0);
+    Writeln(M.ToString);
+    end;
+  if DoTest(39,'M3D.HeadingPitchBank(0,0,Pi/2)') then
+    begin
+    M:=TMatrix3D.CreateRotationHeadingPitchBank(0,0,Pi/2);
+    Writeln(M.ToString);
+    end;
+  if DoTest(40,'M3D.RotationX(Pi/2)') then
+    begin
+    M:=TMatrix3D.CreateRotationX(Pi/2);
+    Writeln(M.ToString);
+    end;
+  if DoTest(41,'M3D.RotationY(Pi/2)') then
+    begin
+    M:=TMatrix3D.CreateRotationY(Pi/2);
+    Writeln(M.ToString);
+    end;
+  if DoTest(42,'M3D.RotationZ(Pi/2)') then
+    begin
+    M:=TMatrix3D.CreateRotationZ(Pi/2);
+    Writeln(M.ToString);
+    end;
+  if DoTest(43,'M3D.YawPitchRoll(Pi/2,0,0)') then
+    begin
+    M:=TMatrix3D.CreateRotationYawPitchRoll(Pi/2,0,0);
+    Writeln(M.ToString);
+    end;
+  if DoTest(44,'M3D.YawPitchRoll(0,Pi/2,0)') then
+    begin
+    M:=TMatrix3D.CreateRotationYawPitchRoll(0,Pi/2,0);
+    Writeln(M.ToString);
+    end;
+  if DoTest(45,'M3D.YawPitchRoll(0,0,Pi/2)') then
+    begin
+    M:=TMatrix3D.CreateRotationYawPitchRoll(0,0,Pi/2);
+    Writeln(M.ToString);
+    end;
+  if DoTest(46,'M3D.LookAtDirLH((0,0,0),(0,0,1),(0,0,1)') then
+    begin
+    M:=TMatrix3D.CreateLookAtDirLH(Point3D(0,0,0),Point3D(0,0,1),Point3D(0,0,1));
+    Writeln(M.ToString);
+    end;
+  if DoTest(47,'M3D.LookAtDirLH((0,0,0),(0,1,0),(0,0,1)') then
+    begin
+    M:=TMatrix3D.CreateLookAtDirLH(Point3D(0,0,0),Point3D(0,1,0),Point3D(0,0,1));
+    Writeln(M.ToString);
+    end;
+  if DoTest(48,'M3D.LookAtDirLH((0,0,0),(1,0,0),(0,0,1)') then
+    begin
+    M:=TMatrix3D.CreateLookAtDirLH(Point3D(0,0,0),Point3D(1,0,0),Point3D(0,0,1));
+    Writeln(M.ToString);
+    end;
+  if DoTest(49,'M3D.LookAtDirRH((0,0,0),(0,0,1),(0,0,1)') then
+    begin
+    M:=TMatrix3D.CreateLookAtDirRH(Point3D(0,0,0),Point3D(0,0,1),Point3D(0,0,1));
+    Writeln(M.ToString);
+    end;
+  if DoTest(50,'M3D.LookAtDirRH((0,0,0),(0,1,0),(0,0,1)') then
+    begin
+    M:=TMatrix3D.CreateLookAtDirRH(Point3D(0,0,0),Point3D(0,1,0),Point3D(0,0,1));
+    Writeln(M.ToString);
+    end;
+  if DoTest(51,'M3D.LookAtDirRH((0,0,0),(1,0,0),(0,0,1)') then
+    begin
+    M:=TMatrix3D.CreateLookAtDirRH(Point3D(0,0,0),Point3D(1,0,0),Point3D(0,0,1));
+    Writeln(M.ToString);
+    end;
+
+  if DoTest(52,'M3D.LookAtLH((0,0,0),(0,0,1),(0,0,1)') then
+    begin
+    M:=TMatrix3D.CreateLookAtLH(Point3D(0,0,0),Point3D(0,0,1),Point3D(0,0,1));
+    Writeln(M.ToString);
+    end;
+  if DoTest(53,'M3D.LookAtLH((0,0,0),(0,1,0),(0,0,1)') then
+    begin
+    M:=TMatrix3D.CreateLookAtLH(Point3D(0,0,0),Point3D(0,1,0),Point3D(0,0,1));
+    Writeln(M.ToString);
+    end;
+  if DoTest(54,'M3D.LookAtLH((0,0,0),(1,0,0),(0,0,1)') then
+    begin
+    M:=TMatrix3D.CreateLookAtLH(Point3D(0,0,0),Point3D(1,0,0),Point3D(0,0,1));
+    Writeln(M.ToString);
+    end;
+  if DoTest(55,'M3D.LookAtRH((0,0,0),(0,0,1),(0,0,1)') then
+    begin
+    Params(Point3D(0,0,0),Point3D(0,0,1),Point3D(0,0,1));
+    M:=TMatrix3D.CreateLookAtRH(Point3D(0,0,0),Point3D(0,0,1),Point3D(0,0,1));
+    Writeln(M.ToString);
+    end;
+  if DoTest(56,'M3D.LookAtRH((0,0,0),(0,1,0),(0,0,1)') then
+    begin
+    M:=TMatrix3D.CreateLookAtRH(Point3D(0,0,0),Point3D(0,1,0),Point3D(0,0,1));
+    Writeln(M.ToString);
+    end;
+  if DoTest(57,'M3D.LookAtRH((0,0,0),(1,0,0),(0,0,1)') then
+    begin
+    M:=TMatrix3D.CreateLookAtRH(Point3D(0,0,0),Point3D(1,0,0),Point3D(0,0,1));
+    Writeln(M.ToString);
+    end;
+  if DoTest(58,'M3D.CreateOrthoLH(1,1,0,1)') then
+    begin
+    M:=TMatrix3D.CreateOrthoLH(1,1,0,1);
+    Writeln(M.ToString);
+    end;
+  if DoTest(59,'M3D.CreateOrthoLH(1,1,0,2)') then
+    begin
+    M:=TMatrix3D.CreateOrthoLH(1,1,0,2);
+    Writeln(M.ToString);
+    end;
+  if DoTest(60,'M3D.CreateOrthoRH(1,1,0,1)') then
+    begin
+    M:=TMatrix3D.CreateOrthoRH(1,1,0,1);
+    Writeln(M.ToString);
+    end;
+  if DoTest(61,'M3D.CreateOrthoRH(1,1,0,2)') then
+    begin
+    M:=TMatrix3D.CreateOrthoRH(1,1,0,2);
+    Writeln(M.ToString);
+    end;
+  if DoTest(62,'M3D.CreateOrthoOffCenterLH(0,0,1,1,0,1)') then
+    begin
+    M:=TMatrix3D.CreateOrthoOffCenterLH(0,0,1,1,0,1);
+    Writeln(M.ToString);
+    end;
+  if DoTest(63,'M3D.CreateOrthoOffCenterLH(0,0,1,1,0,2)') then
+    begin
+    M:=TMatrix3D.CreateOrthoOffCenterLH(0,0,1,1,0,2);
+    Writeln(M.ToString);
+    end;
+  if DoTest(64,'M3D.CreateOrthoOffCenterRH(0,0,1,1,0,1)') then
+    begin
+    M:=TMatrix3D.CreateOrthoOffCenterRH(0,0,1,1,0,1);
+    Writeln(M.ToString);
+    end;
+  if DoTest(65,'M3D.CreateOrthoOffCenterRH(0,0,1,1,0,2)') then
+    begin
+    M:=TMatrix3D.CreateOrthoOffCenterRH(0,0,1,1,0,2);
+    Writeln(M.ToString);
+    end;
+  if DoTest(66,'M3D.CreatePerspectiveFovLH(pi/2,1,0,1)') then
+    begin
+    M:=TMatrix3D.CreatePerspectiveFovLH(pi/2,1,0,1);
+    Writeln(M.ToString);
+    end;
+  if DoTest(67,'M3D.CreatePerspectiveFovLH(pi/21,0,1,True)') then
+    begin
+    M:=TMatrix3D.CreatePerspectiveFovLH(pi/2,1,0,1,True);
+    Writeln(M.ToString);
+    end;
+  if DoTest(68,'M3D.CreatePerspectiveFovRH(pi/2,1,0,1)') then
+    begin
+    M:=TMatrix3D.CreatePerspectiveFovRH(pi/2,1,0,1);
+    Writeln(M.ToString);
+    end;
+  if DoTest(69,'M3D.CreatePerspectiveFovRH(pi/2,1,0,1,True)') then
+    begin
+    M:=TMatrix3D.CreatePerspectiveFovRH(pi/2,1,0,1,True);
+    Writeln(M.ToString);
+    end;
+  if DoTest(70,'M3D.Create(array)') then
+    begin
+    M:=TMatrix3D.Create([1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16]);
+    Writeln(M.ToString);
+    end;
+  if DoTest(71,'M3D.CreateScaling(Point3D(1,2,3))') then
+    begin
+    M:=TMatrix3D.CreateScaling(Point3D(1,2,3));
+    Writeln(M.ToString);
+    end;
+  if DoTest(72,'M3D.CreateTranslation(Point3D(1,2,3))') then
+    begin
+    M:=TMatrix3D.CreateTranslation(Point3D(1,2,3));
+    Writeln(M.ToString);
+    end;
+  if DoTest(73,'M3D*M3D') then
+    begin
+    M:=TMatrix3D.Create( [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16]);
+    M2:=TMatrix3D.Create([16,15,14,13,12,11,10,9,8,7,6,5,4,3,2,1]);
+    M3:=M*M2;
+    Writeln(M3.ToString);
+    end;
+  if DoTest(74,'M3D.Adjoint') then
+    begin
+    M:=TMatrix3D.CreateRotationZ(Pi/3);
+    M2:=M.Adjoint;
+    Writeln(M2.ToString);
+    end;
+  if DoTest(75,'M3D.Determinant') then
+    begin
+    M:=TMatrix3D.CreateRotationZ(Pi/3);
+    Writeln(M.Determinant);
+    end;
+  if DoTest(76,'M3D.Determinant') then
+    begin
+    M:=TMatrix3D.Create(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16);
+//    M:=TMatrix3D.Create(1,0,0,0, 0,1,0,0, 0,0,1,0, 0,0,0,1);
+    Writeln(M.Determinant);
+    end;
+  if DoTest(77,'M3D.Rotation.EyePosition') then
+    begin
+    M:=TMatrix3D.CreateRotationZ(Pi/2);
+    Writeln(M.EyePosition.ToString);
+    end;
+  if DoTest(78,'M3D.Translation.EyePosition') then
+    begin
+    M:=TMatrix3D.CreateRotationZ(Pi/3)*TMatrix3D.CreateTranslation(Point3D(1,0,0));
+    Writeln(M.EyePosition.ToString);
+    end;
+  if DoTest(79,'M3D.Inverse') then
+    begin
+    M:=TMatrix3D.CreateRotationZ(Pi/3);
+    Writeln(M.Inverse.ToString);
+    end;
+      
+    
+end;
+
+Procedure TQ;
+
+var
+  Q1,Q2 : TQuaternion3D;
+
+begin
+{$IFDEF FPC}
+  If DoTest(201,'TQuaternion3D.Create(Point3D(1,0,0),Point3D(0,1,0))') then
+    begin
+    Q1:=TQuaternion3D.Create(Point3D(1,0,0),Point3D(0,1,0));
+    Writeln('Q = ',Q1.ToString);
+    end;
+{$ENDIF}    
+  If DoTest(202,'TQuaternion3D.Create(Point3D(1,0,0),Pi/2)') then
+    begin
+    Q1:=TQuaternion3D.Create(Point3D(1,0,0),Pi/2);
+    Writeln('Q = ',Q1.ToString);
+    end;
+  If DoTest(203,'TQuaternion3D.Create(Point3D(0,1,0),Pi/2)') then
+    begin
+    Q1:=TQuaternion3D.Create(Point3D(0,1,0),Pi/2);
+    Writeln('Q = ',Q1.ToString);
+    end;
+  If DoTest(204,'TQuaternion3D.Create(Point3D(0,0,1),Pi/2)') then
+    begin
+    Q1:=TQuaternion3D.Create(Point3D(0,0,1),Pi/2);
+    Writeln('Q = ',Q1.ToString);
+    end;
+  If DoTest(205,'TQuaternion3D.Create(Pi/2,0,0)') then
+    begin
+    Q1:=TQuaternion3D.Create(Pi/2,0,0);
+    Writeln('Q = ',Q1.ToString);
+    end;
+  If DoTest(206,'TQuaternion3D.Create(0,Pi/2,0)') then
+    begin
+    Q1:=TQuaternion3D.Create(0,Pi/2,0);
+    Writeln('Q = ',Q1.ToString);
+    end;
+  If DoTest(207,'TQuaternion3D.Create(0,0,Pi/2)') then
+    begin
+    Q1:=TQuaternion3D.Create(0,0,Pi/2);
+    Writeln('Q = ',Q1.ToString);
+    end;
+  If DoTest(208,'TQuaternion3D.Multiply (yawpitchroll)') then
+    begin
+    Q1:=TQuaternion3D.Create(0,0,Pi/2);
+    Q2:=TQuaternion3D.Create(0,Pi/2,0);
+    Writeln('Q1*Q2 = ',(Q1*Q2).ToString);
+    end;
+  If DoTest(209,'TQuaternion3D.Multiply (angle/vector)') then
+    begin
+    Q1:=TQuaternion3D.Create(Point3D(1,0,0),Pi/2);
+    Q2:=TQuaternion3D.Create(Point3D(0,1,0),Pi/2);
+    Writeln('Q1*Q2 = ',(Q1*Q2).ToString);
+    end;
+
+end;
+
+begin
+  TV;
+  TP;
+  TM;
+  TM3D;
+  TQ;
+end.

packages/rtl-objpas/tests/utcmatrix.pp

@@ -0,0 +1,328 @@
+{
+    This file is part of the Free Pascal run time library.
+    Copyright (c) 2023 by Michael Van Canneyt
+    member of the Free Pascal development team
+
+    See the file COPYING.FPC, included in this distribution,
+    for details about the copyright.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+ **********************************************************************}
+unit utcmatrix;
+
+{$mode ObjFPC}{$H+}
+
+interface
+
+uses
+  Classes, SysUtils, fpcunit, testregistry, utmathvectorbase, types, system.math.vectors;
+
+Type
+
+  { TTestMatrix }
+
+  TTestMatrix = class(TCMathVectorsBase)
+  Private
+    FM : Array[1..3] of TMatrix;
+    procedure ClearMatrices;
+    function GetM(AIndex: Integer): TMatrix;
+    procedure SetM(AIndex: Integer; AValue: TMatrix);
+  Protected
+    procedure SetUp; override;
+    procedure TearDown; override;
+    Property M1 : TMatrix Index 1 Read GetM Write SetM;
+    Property M2 : TMatrix Index 2 Read GetM Write SetM;
+    Property M3 : TMatrix Index 3 Read GetM Write SetM;
+  Published
+    procedure TestHookUp;
+    Procedure TestCreateRotation;
+    Procedure TestCreateScaling;
+    Procedure TestCreateTranslation;
+    Procedure TTestEqual;
+    Procedure TTestMultiply;
+    Procedure TTestMultiplyPointf;
+    Procedure TTestMultiplyVector;
+    Procedure TTestMultiplyPoint3D;
+    Procedure TTestMultiplyFactor1;
+    Procedure TTestMultiplyFactor2;
+    Procedure TTestDiv;
+    Procedure TTestAdjoint;
+    Procedure TTestDeterminant;
+    Procedure TTestEqualsTo;
+    Procedure TTestExtractScale;
+    Procedure TTestInverse;
+    Procedure TTestScale;
+  end;
+
+
+implementation
+
+{ TTestMatrix }
+
+procedure TTestMatrix.ClearMatrices;
+
+var
+  I : integer;
+
+begin
+  For I:=1 to 3 do
+    FM[I]:=Default(TMatrix);
+end;
+
+function TTestMatrix.GetM(AIndex: Integer): TMatrix;
+begin
+  Result:=FM[aIndex];
+end;
+
+procedure TTestMatrix.SetM(AIndex: Integer; AValue: TMatrix);
+begin
+  FM[aIndex]:=aValue;
+end;
+
+procedure TTestMatrix.SetUp;
+begin
+  inherited SetUp;
+  ClearMatrices;
+end;
+
+procedure TTestMatrix.TearDown;
+begin
+  inherited TearDown;
+  ClearMatrices;
+end;
+
+procedure TTestMatrix.TestHookUp;
+
+var
+  I,C,R : Integer;
+begin
+  For I:=1 to 3 do
+    For R:=0 to 2 do
+      For C:=0 to 2 do
+       AssertEquals(Format('M%d[%d,%d]',[I,C,R]),0.0,FM[I].M[R].V[C]);
+end;
+
+procedure TTestMatrix.TestCreateRotation;
+
+begin
+  M1:=TMatrix.CreateRotation(0);
+  AssertMatrix('Create 1',[1,0,0,
+                           0,1,0,
+                           0,0,1],M1);
+  M1:=TMatrix.CreateRotation(Pi/2);
+  AssertMatrix('Create 2',[0,1,0,
+                           -1,0,0,
+                           0,0,1],M1);
+  M1:=TMatrix.CreateRotation(Pi);
+  AssertMatrix('Create 3',[-1,0,0,
+                           0,-1,0,
+                           0,0,1],M1);
+end;
+
+procedure TTestMatrix.TestCreateScaling;
+begin
+  M1:=TMatrix.CreateScaling(2,2);
+  AssertMatrix('Create 1',[2,0,0,
+                           0,2,0,
+                           0,0,1],M1);
+  M1:=TMatrix.CreateScaling(2,3);
+  AssertMatrix('Create 2',[2,0,0,
+                           0,3,0,
+                           0,0,1],M1);
+end;
+
+procedure TTestMatrix.TestCreateTranslation;
+begin
+  M1:=TMatrix.CreateTranslation(2,3);
+  AssertMatrix('Create 1',[1,0,0,
+                           0,1,0,
+                           2,3,1],M1);
+  M1:=TMatrix.CreateTranslation(0,0);
+  AssertMatrix('Create 2',[1,0,0,
+                           0,1,0,
+                           0,0,1],M1);
+end;
+
+procedure TTestMatrix.TTestEqual;
+begin
+  M1:=TMatrix.CreateTranslation(2,3);
+  M2:=TMatrix.CreateTranslation(2,3);
+  AssertTrue('Equal translation',M1=M2);
+  M2:=TMatrix.CreateTranslation(0,0);
+  AssertFalse('NEqual translation',M1=M2);
+  M1:=TMatrix.CreateScaling(2,3);
+  M2:=TMatrix.CreateScaling(2,3);
+  AssertTrue('Equal Scaling',M1=M2);
+  M2:=TMatrix.CreateScaling(0,0);
+  AssertFalse('NEqual Scaling',M1=M2);
+end;
+
+procedure TTestMatrix.TTestMultiply;
+begin
+  M1:=TMatrix.CreateRotation(Pi/4);
+  M2:=TMatrix.CreateRotation(Pi/4);
+  M3:=M1*M2;
+  AssertMatrix('Multiply',[0,1,0,
+                           -1,0,0,
+                           0,0,1],M3);
+  M1:=TMatrix.CreateRotation(Pi/3);
+  M2:=TMatrix.CreateRotation(Pi/6);
+  M3:=M1*M2;
+  AssertMatrix('Multiply',[0,1,0,
+                           -1,0,0,
+                           0,0,1],M3);
+end;
+
+procedure TTestMatrix.TTestMultiplyPointf;
+
+Var
+  P1,P2 : TPointF;
+
+begin
+  P1:=PointF(1,0);
+  M1:=TMatrix.CreateRotation(Pi/2);
+  P2:=P1*M1;
+  AssertPointF('Multiplied',0,1,P2);
+  M1:=TMatrix.CreateRotation(Pi);
+  P2:=P1*M1;
+  AssertPointF('Multiplied 2',-1,0,P2);
+  M1:=TMatrix.CreateTranslation(3,4);
+  P2:=P1*M1;
+  AssertPointF('Multiplied 2',4,4,P2);
+end;
+
+procedure TTestMatrix.TTestMultiplyVector;
+
+Var
+  V1,V2 : TVector;
+
+begin
+  V1:=Vector(1,0);
+  M1:=TMatrix.CreateRotation(Pi/2);
+  V2:=V1*M1;
+  AssertVector('Multiplied',0,1,1,V2);
+  M1:=TMatrix.CreateRotation(Pi);
+  V2:=V1*M1;
+  AssertVector('Multiplied 2',-1,0,1,V2);
+end;
+
+procedure TTestMatrix.TTestMultiplyPoint3D;
+Var
+  P1,P2 : TPoint3D;
+
+begin
+  P1:=Point3D(1,0,0);
+  M1:=TMatrix.CreateRotation(Pi/2);
+  P2:=P1*M1;
+  AssertPoint3D('Multiplied',0,1,0,P2);
+  M1:=TMatrix.CreateRotation(Pi);
+  P2:=P1*M1;
+  AssertPoint3D('Multiplied 2',-1,0,0,P2);
+end;
+
+procedure TTestMatrix.TTestMultiplyFactor1;
+begin
+  M1:=TMatrix.CreateRotation(Pi/2);
+  M2:=M1*2;
+  AssertMatrix('Multiply',[0,2,0,
+                           -2,0,0,
+                           0,0,2],M2);
+end;
+
+procedure TTestMatrix.TTestMultiplyFactor2;
+begin
+  M1:=TMatrix.CreateRotation(Pi/2);
+  M2:=2*M1;
+  AssertMatrix('Multiply',[0,2,0,
+                           -2,0,0,
+                           0,0,2],M2);
+end;
+
+procedure TTestMatrix.TTestDiv;
+begin
+  M1:=TMatrix.CreateRotation(Pi/2);
+  M2:=M1/2;
+  AssertMatrix('Divide',[ 0,  1/2,0,
+                         -1/2, 0, 0,
+                          0,   0, 1/2],M2);
+end;
+
+procedure TTestMatrix.TTestAdjoint;
+begin
+  M1:=TMatrix.CreateRotation(Pi/2);
+  M2:=M1.Adjoint;
+  AssertMatrix('Adjoint 1',[0,-1,0,
+                            1,0,0,
+                            0,0,1],M2);
+  M1:=TMatrix.CreateRotation(Pi/6);
+  M2:=M1.Adjoint;
+  AssertMatrix('Adjoint 2',[ 0.866, -0.5, 0.00,
+                             0.50, 0.866, 0.00,
+                            -0.00, 0.00, 1.00],M2);
+end;
+
+procedure TTestMatrix.TTestDeterminant;
+begin
+  M1:=TMatrix.CreateRotation(Pi/2);
+  AssertEquals('Det 1',1,M1.Determinant);
+  M1:=TMatrix.CreateRotation(Pi/3);
+  AssertEquals('Det 1',1,M1.Determinant);
+end;
+
+procedure TTestMatrix.TTestEqualsTo;
+begin
+  M1:=TMatrix.CreateTranslation(2,3);
+  M2:=TMatrix.CreateTranslation(2,3);
+  AssertTrue('Equal translation',M1.EqualsTo(M2));
+  M2:=TMatrix.CreateTranslation(0,0);
+  AssertFalse('NEqual translation',M1.EqualsTo(M2));
+  M1:=TMatrix.CreateScaling(2,3);
+  M2:=TMatrix.CreateScaling(2,3);
+  AssertTrue('Equal Scaling',M1.EqualsTo(M2));
+  M2:=TMatrix.CreateScaling(0,0);
+  AssertFalse('NEqual Scaling',M1.EqualsTo(M2));
+
+end;
+
+procedure TTestMatrix.TTestExtractScale;
+begin
+  M1:=TMatrix.CreateRotation(Pi/2);
+  AssertPointF('Scale 1',1,1,M1.ExtractScale);
+  M1:=M1.Scale(2);
+  AssertPointF('Scale 2',2,2,M1.ExtractScale);
+end;
+
+procedure TTestMatrix.TTestInverse;
+begin
+  M1:=TMatrix.CreateRotation(Pi/2);
+  M2:=M1.Inverse;
+  M3:=TMatrix.CreateRotation(-Pi/2);
+  AssertMatrix('Inverse 1',M3,M2);
+  M1:=TMatrix.CreateRotation(Pi/6);
+  M2:=M1.Inverse;
+  M3:=TMatrix.CreateRotation(-Pi/6);
+  AssertMatrix('Inverse 2',M3,M2);
+end;
+
+procedure TTestMatrix.TTestScale;
+
+
+begin
+  M1:=TMatrix.CreateRotation(Pi/2);
+  AssertMatrix('Create',[0,1,0,
+                         -1,0,0,
+                         0,0,1],M1);
+  M2:=M1.Scale(2);
+  AssertMatrix('Scaled',[0,2,0,
+                         -2,0,0,
+                         0,0,2],M2);
+end;
+
+
+initialization
+  RegisterTest(TTestMatrix);
+end.
+

packages/rtl-objpas/tests/utcmatrix3d.pp

@@ -0,0 +1,636 @@
+{
+    This file is part of the Free Pascal run time library.
+    Copyright (c) 2023 by Michael Van Canneyt
+    member of the Free Pascal development team
+
+    See the file COPYING.FPC, included in this distribution,
+    for details about the copyright.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+ **********************************************************************}
+unit utcmatrix3d;
+
+{$mode ObjFPC}{$H+}
+
+interface
+
+uses
+  Classes, SysUtils, fpcunit, testregistry, utmathvectorbase, types, system.math.vectors;
+
+Type
+  { TTestMatrix3D }
+
+  TTestMatrix3D = class(TCMathVectorsBase)
+  Private
+    FM : Array[1..3] of TMatrix3D;
+    procedure ClearMatrices;
+    function GetM(AIndex: Integer): TMatrix3D;
+    procedure SetM(AIndex: Integer; AValue: TMatrix3D);
+  Protected
+    procedure SetUp; override;
+    procedure TearDown; override;
+    Property M1 : TMatrix3D Index 1 Read GetM Write SetM;
+    Property M2 : TMatrix3D Index 2 Read GetM Write SetM;
+    Property M3 : TMatrix3D Index 3 Read GetM Write SetM;
+  Published
+    procedure TestHookUp;
+    constructor TestCreate;
+    constructor TestCreateArray;
+
+    Procedure TestZero;
+
+    Procedure TestCreateLookAtDirLH;
+    Procedure TestCreateLookAtDirRH;
+    Procedure TestCreateLookAtLH;
+    Procedure TestCreateLookAtRH;
+
+    Procedure TestCreateOrthoLH;
+    Procedure TestCreateOrthoOffCenterLH;
+    Procedure TestCreateOrthoOffCenterRH;
+    Procedure TestCreateOrthoRH;
+
+    Procedure TestCreatePerspectiveFovLH;
+    Procedure TestCreatePerspectiveFovRH;
+
+    Procedure TestCreateRotation;
+    Procedure TestCreateRotationX;
+    Procedure TestCreateRotationY;
+    Procedure TestCreateRotationZ;
+    Procedure TestCreateRotationHeadingPitchBank;
+    Procedure TestCreateRotationYawPitchRoll;
+    Procedure TestCreateScaling;
+    Procedure TestCreateTranslation;
+
+    Procedure TestMultiplyPoint;
+    Procedure TestMultiplyMatrix;
+    Procedure TestMultiplyVector3D;
+    Procedure TestMultiplyFactor;
+    Procedure TestMultiplyFactor2;
+    Procedure TestDiv;
+
+    Procedure TestAdjoint;
+    Procedure TestDeterminant;
+    Procedure TestEyePosition;
+    Procedure TestInverse;
+    Procedure TestScale;
+    Procedure TestToMatrix;
+    Procedure TestTranspose;
+  end;
+
+
+implementation
+
+procedure TTestMatrix3D.ClearMatrices;
+
+var
+  I : integer;
+
+begin
+  For I:=1 to 3 do
+    FM[I]:=Default(TMatrix3D);
+end;
+
+Function TTestMatrix3D.GetM(AIndex: Integer): TMatrix3D;
+begin
+  Result:=FM[aIndex];
+end;
+
+procedure TTestMatrix3D.SetM(AIndex: Integer; AValue: TMatrix3D);
+begin
+  FM[aIndex]:=aValue;
+end;
+
+procedure TTestMatrix3D.SetUp;
+begin
+  inherited SetUp;
+  ClearMatrices;
+end;
+
+procedure TTestMatrix3D.TearDown;
+begin
+  inherited TearDown;
+  ClearMatrices;
+end;
+
+procedure TTestMatrix3D.TestHookUp;
+
+var
+  I,C,R : Integer;
+
+begin
+  For I:=1 to 3 do
+    For R:=0 to 3 do
+      For C:=0 to 3 do
+        AssertEquals(Format('M%d[%d,%d]',[I,C,R]),0.0,FM[I].M[R].V[C]);
+end;
+
+constructor TTestMatrix3D.TestCreate;
+begin
+  M1:=TMatrix3D.Create(1,2,3,4,
+                       5,6,7,8,
+                       9,10,11,12,
+                       13,14,15,16);
+  // actual test
+  AssertMatrix3D('Create',1,2,3,4,
+                          5,6,7,8,
+                          9,10,11,12,
+                          13,14,15,16,M1);
+  // Test assert using array
+  AssertMatrix3D('Create',[1,2,3,4,
+                          5,6,7,8,
+                          9,10,11,12,
+                          13,14,15,16],M1);
+  // Test assert using second matrix
+  M2:=M1;
+  AssertMatrix3D('Create',M2,M1);
+end;
+
+constructor TTestMatrix3D.TestCreateArray;
+begin
+  M1:=TMatrix3D.Create([ 1,  2,  3,  4,
+                         5,  6,  7,  8,
+                         9, 10, 11, 12,
+                        13, 14, 15, 16]);
+
+  AssertMatrix3D('Create',1,  5,  9, 13,
+                          2,  6, 10, 14,
+                          3,  7, 11, 15,
+                          4,  8, 12, 16, M1);
+end;
+
+procedure TTestMatrix3D.TestZero;
+begin
+  M1:=TMatrix3D.Zero;
+  AssertMatrix3D('Create',0,0,0,0,
+                          0,0,0,0,
+                          0,0,0,0,
+                          0,0,0,0, M1);
+end;
+
+procedure TTestMatrix3D.TestCreateLookAtDirLH;
+begin
+  M1:=TMatrix3D.CreateLookAtDirLH(Point3D(0,0,0),Point3D(0,0,1),Point3D(0,0,1));
+  // test 46
+  AssertMatrix3D('Create 1',[ 0.0000, 0.0000, 0.0000, 0.0000,
+                            0.0000, 0.0000, 0.0000, 0.0000,
+                            0.0000, 0.0000,-1.0000, 0.0000,
+                            0.0000, 0.0000, 0.0000, 1.0000],M1);
+  // test 47
+  M1:=TMatrix3D.CreateLookAtDirLH(Point3D(0,0,0),Point3D(0,1,0),Point3D(0,0,1));
+  AssertMatrix3D('Create 2',[ 1.0000, 0.0000, 0.0000, 0.0000,
+                              0.0000, 0.0000,-1.0000, 0.0000,
+                              0.0000, 1.0000, 0.0000, 0.0000,
+                              0.0000, 0.0000, 0.0000, 1.0000],M1);
+  // test 48
+  M1:=TMatrix3D.CreateLookAtDirLH(Point3D(0,0,0),Point3D(1,0,0),Point3D(0,0,1));
+  AssertMatrix3D('Create 3',[ 0.0000, 0.0000,-1.0000, 0.0000,
+                             -1.0000, 0.0000, 0.0000, 0.0000,
+                             0.0000, 1.0000, 0.0000, 0.0000,
+                             0.0000, 0.0000, 0.0000, 1.0000],M1);
+
+end;
+
+procedure TTestMatrix3D.TestCreateLookAtDirRH;
+begin
+  // test 49
+  M1:=TMatrix3D.CreateLookAtDirRH(Point3D(0,0,0),Point3D(0,0,1),Point3D(0,0,1));
+  AssertMatrix3D('Create 1',[ 0.0000, 0.0000, 0.0000, 0.0000,
+                              0.0000, 0.0000, 0.0000, 0.0000,
+                              0.0000, 0.0000, 1.0000, 0.0000,
+                              0.0000, 0.0000, 0.0000, 1.0000],M1);
+  // test 50
+  M1:=TMatrix3D.CreateLookAtDirRH(Point3D(0,0,0),Point3D(0,1,0),Point3D(0,0,1));
+  AssertMatrix3D('Create 2',[ -1.0000, 0.0000, 0.0000, 0.0000,
+                               0.0000, 0.0000, 1.0000, 0.0000,
+                               0.0000, 1.0000, 0.0000, 0.0000,
+                               0.0000, 0.0000, 0.0000, 1.0000],M1);
+  // test 51
+  M1:=TMatrix3D.CreateLookAtDirRH(Point3D(0,0,0),Point3D(1,0,0),Point3D(0,0,1));
+  AssertMatrix3D('Create 3',[ 0.0000, 0.0000, 1.0000, 0.0000,
+                              1.0000, 0.0000, 0.0000, 0.0000,
+                              0.0000, 1.0000, 0.0000, 0.0000,
+                              0.0000, 0.0000, 0.0000, 1.0000],M1);
+end;
+
+procedure TTestMatrix3D.TestCreateLookAtLH;
+begin
+  M1:=TMatrix3D.CreateLookAtLH(Point3D(0,0,0),Point3D(0,0,1),Point3D(0,0,1));
+  // test 52
+  AssertMatrix3D('Create 1',[  0.0000, 0.0000, 0.0000, 0.0000,
+                               0.0000, 0.0000, 0.0000, 0.0000,
+                               0.0000, 0.0000, 1.0000, 0.0000,
+                               0.0000, 0.0000, 0.0000, 1.0000],M1);
+  // test 53
+  M1:=TMatrix3D.CreateLookAtLH(Point3D(0,0,0),Point3D(0,1,0),Point3D(0,0,1));
+  AssertMatrix3D('Create 2',[ -1.0000, 0.0000, 0.0000, 0.0000,
+                               0.0000, 0.0000, 1.0000, 0.0000,
+                               0.0000, 1.0000, 0.0000, 0.0000,
+                               0.0000, 0.0000, 0.0000, 1.0000],M1);
+  // test 54
+  M1:=TMatrix3D.CreateLookAtLH(Point3D(0,0,0),Point3D(1,0,0),Point3D(0,0,1));
+  AssertMatrix3D('Create 3',[ 0.0000, 0.0000, 1.0000, 0.0000,
+                              1.0000, 0.0000, 0.0000, 0.0000,
+                              0.0000, 1.0000, 0.0000, 0.0000,
+                              0.0000, 0.0000, 0.0000, 1.0000],M1);
+
+end;
+
+procedure TTestMatrix3D.TestCreateLookAtRH;
+begin
+  M1:=TMatrix3D.CreateLookAtRH(Point3D(0,0,0),Point3D(0,0,1),Point3D(0,0,1));
+  // test 55
+  AssertMatrix3D('Create 1',[  0.0000, 0.0000, 0.0000, 0.0000,
+                               0.0000, 0.0000, 0.0000, 0.0000,
+                               0.0000, 0.0000,-1.0000, 0.0000,
+                               0.0000, 0.0000, 0.0000, 1.0000],M1);
+  // test 56
+  M1:=TMatrix3D.CreateLookAtRH(Point3D(0,0,0),Point3D(0,1,0),Point3D(0,0,1));
+  AssertMatrix3D('Create 2',[  1.0000, 0.0000, 0.0000, 0.0000,
+                               0.0000, 0.0000,-1.0000, 0.0000,
+                               0.0000, 1.0000, 0.0000, 0.0000,
+                               0.0000, 0.0000, 0.0000, 1.0000],M1);
+  // test 57
+  M1:=TMatrix3D.CreateLookAtRH(Point3D(0,0,0),Point3D(1,0,0),Point3D(0,0,1));
+  AssertMatrix3D('Create 3',[  0.0000, 0.0000,-1.0000, 0.0000,
+                               -1.0000, 0.0000, 0.0000, 0.0000,
+                                0.0000, 1.0000, 0.0000, 0.0000,
+                                0.0000, 0.0000, 0.0000, 1.0000],M1);
+end;
+
+procedure TTestMatrix3D.TestCreateOrthoLH;
+begin
+  // test 58
+  M1:=TMatrix3D.CreateOrthoLH(1,1,0,1);
+  AssertMatrix3D('Create 1',[  2.0000, 0.0000, 0.0000, 0.0000,
+                               0.0000, 2.0000, 0.0000, 0.0000,
+                               0.0000, 0.0000, 1.0000, 0.0000,
+                               0.0000, 0.0000, 0.0000, 1.0000],M1);
+  // test 59
+  M1:=TMatrix3D.CreateOrthoLH(1,1,0,2);
+  AssertMatrix3D('Create 2',[2.0000, 0.0000, 0.0000, 0.0000,
+                             0.0000, 2.0000, 0.0000, 0.0000,
+                             0.0000, 0.0000, 0.5000, 0.0000,
+                             0.0000, 0.0000, 0.0000, 1.0000],M1);
+end;
+
+procedure TTestMatrix3D.TestCreateOrthoRH;
+begin
+  // test 60
+  M1:=TMatrix3D.CreateOrthoRH(1,1,0,1);
+  AssertMatrix3D('Create 1',[  2.0000, 0.0000, 0.0000, 0.0000,
+                               0.0000, 2.0000, 0.0000, 0.0000,
+                               0.0000, 0.0000,-1.0000, 0.0000,
+                               0.0000, 0.0000, 0.0000, 1.0000],M1);
+  // test 61
+  M1:=TMatrix3D.CreateOrthoRH(1,1,0,2);
+  AssertMatrix3D('Create 2',[  2.0000, 0.0000, 0.0000, 0.0000,
+                               0.0000, 2.0000, 0.0000, 0.0000,
+                               0.0000, 0.0000,-0.5000, 0.0000,
+                               0.0000, 0.0000, 0.0000, 1.0000],M1);
+end;
+
+procedure TTestMatrix3D.TestCreateOrthoOffCenterLH;
+begin
+  // test 62
+  M1:=TMatrix3D.CreateOrthoOffCenterLH(0,0,1,1,0,1);
+  AssertMatrix3D('Create 1',[2.0000, 0.0000, 0.0000, 0.0000,
+                             0.0000,-2.0000, 0.0000, 0.0000,
+                             0.0000, 0.0000, 1.0000, 0.0000,
+                             -1.0000, 1.0000, 0.0000, 1.0000],M1);
+  // test 63
+  M1:=TMatrix3D.CreateOrthoOffCenterLH(0,0,1,1,0,2);
+  AssertMatrix3D('Create 2',[2.0000, 0.0000, 0.0000, 0.0000,
+                             0.0000,-2.0000, 0.0000, 0.0000,
+                             0.0000, 0.0000, 0.5000, 0.0000,
+                             -1.0000, 1.0000, 0.0000, 1.0000],M1);
+end;
+
+procedure TTestMatrix3D.TestCreateOrthoOffCenterRH;
+begin
+  // test 64
+  M1:=TMatrix3D.CreateOrthoOffCenterRH(0,0,1,1,0,1);
+  AssertMatrix3D('Create 1',[2.0000, 0.0000, 0.0000, 0.0000,
+                             0.0000,-2.0000, 0.0000, 0.0000,
+                             0.0000, 0.0000,-1.0000, 0.0000,
+                             -1.0000, 1.0000, 0.0000, 1.0000],M1);
+  // test 65
+  M1:=TMatrix3D.CreateOrthoOffCenterRH(0,0,1,1,0,2);
+  AssertMatrix3D('Create 2',[ 2.0000, 0.0000, 0.0000, 0.0000,
+                              0.0000,-2.0000, 0.0000, 0.0000,
+                              0.0000, 0.0000,-0.5000, 0.0000,
+                              -1.0000, 1.0000, 0.0000, 1.0000],M1);
+end;
+
+
+procedure TTestMatrix3D.TestCreatePerspectiveFovLH;
+begin
+  // Test 66
+  M1:=TMatrix3D.CreatePerspectiveFovLH(pi/2,1,0,1);
+  AssertMatrix3D('Create 1',[ 1.0000, 0.0000, 0.0000, 0.0000,
+                              0.0000, 1.0000, 0.0000, 0.0000,
+                              0.0000, 0.0000, 1.0000, 1.0000,
+                              0.0000, 0.0000, 0.0000, 0.0000],M1);
+  // Test 67
+  M1:=TMatrix3D.CreatePerspectiveFovLH(pi/2,1,0,1,True);
+  AssertMatrix3D('Create 2',[ 1.0000, 0.0000, 0.0000, 0.0000,
+                              0.0000, 1.0000, 0.0000, 0.0000,
+                              0.0000, 0.0000, 1.0000, 1.0000,
+                              0.0000, 0.0000,-0.0000, 0.0000], M1);
+end;
+
+procedure TTestMatrix3D.TestCreatePerspectiveFovRH;
+begin
+  // Test 68
+  M1:=TMatrix3D.CreatePerspectiveFovRH(pi/2,1,0,1);
+  AssertMatrix3D('Create 1',[ 1.0000, 0.0000, 0.0000, 0.0000,
+                              0.0000, 1.0000, 0.0000, 0.0000,
+                              0.0000, 0.0000,-1.0000,-1.0000,
+                              0.0000, 0.0000, 0.0000, 0.0000],M1);
+  // Test 69
+  M1:=TMatrix3D.CreatePerspectiveFovRH(pi/2,1,0,1,True);
+  AssertMatrix3D('Create 2',[ 1.0000, 0.0000, 0.0000, 0.0000,
+                              0.0000, 1.0000, 0.0000, 0.0000,
+                              0.0000, 0.0000,-1.0000,-1.0000,
+                              0.0000, 0.0000, 0.0000, 0.0000 ],M1);
+end;
+
+procedure TTestMatrix3D.TestCreateRotation;
+begin
+  // test 32
+  M1:=TMatrix3D.CreateRotation(Point3D(1,0,0),Pi/2);
+  AssertMatrix3D('Create 1',[ 1.0000, 0.0000, 0.0000, 0.0000,
+                              0.0000, 0.0000, 1.0000, 0.0000,
+                              0.0000,-1.0000, 0.0000, 0.0000,
+                              0.0000, 0.0000, 0.0000, 1.0000],M1);
+  // test 33
+  M1:=TMatrix3D.CreateRotation(Point3D(0,1,0),Pi/2);
+  AssertMatrix3D('Create 2',[ 0.0000, 0.0000,-1.0000, 0.0000,
+                              0.0000, 1.0000, 0.0000, 0.0000,
+                              1.0000, 0.0000, 0.0000, 0.0000,
+                              0.0000, 0.0000, 0.0000, 1.0000],M1);
+  // test 34
+  M1:=TMatrix3D.CreateRotation(Point3D(0,0,1),pi/2);
+  AssertMatrix3D('Create 3',[ -0.0000, 1.0000, 0.0000, 0.0000,
+                              -1.0000,-0.0000, 0.0000, 0.0000,
+                              0.0000, 0.0000, 1.0000, 0.0000,
+                              0.0000, 0.0000, 0.0000, 1.0000],M1);
+end;
+
+procedure TTestMatrix3D.TestCreateRotationX;
+begin
+  // test 40
+  M1:=TMatrix3D.CreateRotationX(Pi/2);
+  AssertMatrix3D('Create 1',[ 1.0000, 0.0000, 0.0000, 0.0000,
+                              0.0000, 0.0000, 1.0000, 0.0000,
+                              0.0000,-1.0000, 0.0000, 0.0000,
+                              0.0000, 0.0000, 0.0000, 1.0000],M1);
+
+
+end;
+
+procedure TTestMatrix3D.TestCreateRotationY;
+begin
+  // test 41
+  M1:=TMatrix3D.CreateRotationY(Pi/2);
+  AssertMatrix3D('Create 1',[ 0.0000, 0.0000,-1.0000, 0.0000,
+                              0.0000, 1.0000, 0.0000, 0.0000,
+                              1.0000, 0.0000, 0.0000, 0.0000,
+                              0.0000, 0.0000, 0.0000, 1.0000],M1);
+
+end;
+
+procedure TTestMatrix3D.TestCreateRotationZ;
+begin
+  // test 42
+  M1:=TMatrix3D.CreateRotationZ(Pi/2);
+  AssertMatrix3D('Create 1',[ -0.0000, 1.0000, 0.0000, 0.0000,
+                              -1.0000,-0.0000, 0.0000, 0.0000,
+                              0.0000, 0.0000, 1.0000, 0.0000,
+                              0.0000, 0.0000, 0.0000, 1.0000],M1);
+end;
+
+procedure TTestMatrix3D.TestCreateRotationHeadingPitchBank;
+begin
+  // Test 37
+  M1:=TMatrix3D.CreateRotationHeadingPitchBank(Pi/2,0,0);
+  AssertMatrix3D('Create 1',[ 0.0000, 0.0000, 1.0000, 0.0000,
+                              0.0000, 1.0000, 0.0000, 0.0000,
+                              -1.0000, 0.0000, 0.0000, 0.0000,
+                              0.0000, 0.0000, 0.0000, 1.0000],M1);
+  // Test 38
+  M1:=TMatrix3D.CreateRotationHeadingPitchBank(0,Pi/2,0);
+  AssertMatrix3D('Create 2',[ 1.0000, 0.0000, 0.0000, 0.0000,
+                              0.0000, 0.0000,-1.0000, 0.0000,
+                              0.0000, 1.0000, 0.0000, 0.0000,
+                              0.0000, 0.0000, 0.0000, 1.0000],M1);
+  // Test 39
+  M1:=TMatrix3D.CreateRotationHeadingPitchBank(0,0,Pi/2);
+  AssertMatrix3D('Create 3',[ 0.0000,-1.0000, 0.0000, 0.0000,
+                              1.0000, 0.0000, 0.0000, 0.0000,
+                              0.0000, 0.0000, 1.0000, 0.0000,
+                              0.0000, 0.0000, 0.0000, 1.0000 ],M1);
+end;
+
+procedure TTestMatrix3D.TestCreateRotationYawPitchRoll;
+begin
+  // test 43
+  M1:=TMatrix3D.CreateRotationYawPitchRoll(Pi/2,0,0);
+  AssertMatrix3D('Create 1',[ 0.0000,-1.0000, 0.0000, 0.0000,
+                              1.0000, 0.0000, 0.0000, 0.0000,
+                              0.0000, 0.0000, 1.0000, 0.0000,
+                              0.0000, 0.0000, 0.0000, 1.0000 ],M1);
+  // test 44
+  M1:=TMatrix3D.CreateRotationYawPitchRoll(0,Pi/2,0);
+  AssertMatrix3D('Create 2',[ 1.0000, 0.0000, 0.0000, 0.0000,
+                              0.0000, 0.0000, 1.0000, 0.0000,
+                              0.0000,-1.0000, 0.0000, 0.0000,
+                              0.0000, 0.0000, 0.0000, 1.0000],M1);
+  // test 45
+  M1:=TMatrix3D.CreateRotationYawPitchRoll(0,0,Pi/2);
+  AssertMatrix3D('Create 3',[ 0.0000, 0.0000,-1.0000, 0.0000,
+                              0.0000, 1.0000, 0.0000, 0.0000,
+                              1.0000, 0.0000, 0.0000, 0.0000,
+                              0.0000, 0.0000, 0.0000, 1.0000],M1);
+end;
+
+procedure TTestMatrix3D.TestCreateScaling;
+begin
+  // test 71
+  M1:=TMatrix3D.CreateScaling(Point3D(1,2,3));
+  AssertMatrix3D('Create 1',[ 1, 0, 0, 0,
+                              0, 2, 0, 0,
+                              0, 0, 3, 0,
+                              0, 0, 0, 1],M1);
+
+end;
+
+procedure TTestMatrix3D.TestCreateTranslation;
+begin
+  // test 72
+  M1:=TMatrix3D.CreateTranslation(Point3D(1,2,3));
+  AssertMatrix3D('Create 1',[ 1.0000, 0.0000, 0.0000, 0.0000,
+                              0.0000, 1.0000, 0.0000, 0.0000,
+                              0.0000, 0.0000, 1.0000, 0.0000,
+                              1.0000, 2.0000, 3.0000, 1.0000 ],M1);
+end;
+
+procedure TTestMatrix3D.TestMultiplyPoint;
+var
+  P1,P2 : TPoint3D;
+begin
+  P1:=Point3D(1,0,0);
+  M1:=TMatrix3D.CreateRotationZ(pi/2);
+  P2:=P1*M1;
+  AssertPoint3D('Multiply 1',0,1,0,P2);
+  M1:=TMatrix3D.CreateRotationY(pi/2);
+  P2:=P1*M1;
+  AssertPoint3D('Multiply 2',0,0,-1,P2);
+  P1:=Point3D(0,1,0);
+  M1:=TMatrix3D.CreateRotationX(pi/2);
+  P2:=P1*M1;
+  AssertPoint3D('Multiply 3',0,0,1,P2);
+end;
+
+procedure TTestMatrix3D.TestMultiplyMatrix;
+begin
+  // Test 73
+  M1:=TMatrix3D.Create( [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16]);
+  M2:=TMatrix3D.Create([16,15,14,13,12,11,10,9,8,7,6,5,4,3,2,1]);
+  M3:=M1*M2;
+  AssertMatrix3D('Multiplied',[ 386.0000,274.0000,162.0000,50.0000,
+                                444.0000,316.0000,188.0000,60.0000,
+                                502.0000,358.0000,214.0000,70.0000,
+                                560.0000,400.0000,240.0000,80.0000],m3);
+end;
+
+procedure TTestMatrix3D.TestMultiplyVector3D;
+var
+  P1,P2 : TVector3D;
+begin
+  P1:=Vector3D(1,0,0);
+  M1:=TMatrix3D.CreateRotationZ(pi/2);
+  P2:=P1*M1;
+  AssertVector3D('Multiply 1',0,1,0,1,P2);
+  M1:=TMatrix3D.CreateRotationY(pi/2);
+  P2:=P1*M1;
+  AssertVector3D('Multiply 2',0,0,-1,1,P2);
+  P1:=Point3D(0,1,0);
+  M1:=TMatrix3D.CreateRotationX(pi/2);
+  P2:=P1*M1;
+  AssertVector3D('Multiply 3',0,0,1,1,P2);
+end;
+
+procedure TTestMatrix3D.TestMultiplyFactor;
+begin
+  M1:=TMatrix3D.Create( 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16);
+  M2:=2*M1;
+  AssertMatrix3D('Multiplied',[ 2, 4, 6, 8,
+                                10, 12, 14, 16,
+                                18, 20, 22, 24,
+                                26, 28, 30, 32],M2);
+end;
+
+procedure TTestMatrix3D.TestMultiplyFactor2;
+begin
+  M1:=TMatrix3D.Create( 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16);
+  M3:=M1*2;
+  AssertMatrix3D('Multiplied',[ 2, 4, 6, 8,
+                                10, 12, 14, 16,
+                                18, 20, 22, 24,
+                                26, 28, 30, 32],M3);
+end;
+
+procedure TTestMatrix3D.TestDiv;
+begin
+  M1:=TMatrix3D.Create( 2, 4, 6, 8,
+                        10, 12, 14, 16,
+                        18, 20, 22, 24,
+                        26, 28, 30, 32);
+  M3:=M1/2;
+  AssertMatrix3D('Divided',[ 1 ,  2,  3,  4,
+                                5 ,  6,  7,  8,
+                                9 , 10, 11, 12,
+                                13, 14, 15, 16 ],M3);
+
+end;
+
+procedure TTestMatrix3D.TestAdjoint;
+begin
+  // Test 74
+  M1:=TMatrix3D.CreateRotationZ(Pi/3);
+  M2:=M1.Adjoint;
+  AssertMatrix3D('Adjoint',[0.5000,-0.8660, 0.0000, 0.0000,
+                            0.8660, 0.5000, 0.0000, 0.0000,
+                            0.0000, 0.0000, 1.0000, 0.0000,
+                            0.0000, 0.0000, 0.0000, 1.0000],M2);
+end;
+
+procedure TTestMatrix3D.TestDeterminant;
+begin
+  // Test 75
+  M1:=TMatrix3D.CreateRotationZ(Pi/3);
+  AssertEquals('Determinant',1.0,M1.Determinant);
+  // Test 76
+  M1:=TMatrix3D.Create(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16);
+  AssertEquals('Determinant',0,M1.Determinant);
+  M1:=TMatrix3D.Create(1,0,0,0, 0,1,0,0, 0,0,1,0, 0,0,0,1);
+  AssertEquals('Determinant',1,M1.Determinant);
+end;
+
+procedure TTestMatrix3D.TestEyePosition;
+begin
+  // Test 77
+  M1:=TMatrix3D.CreateRotationZ(Pi/3);
+  AssertPoint3D('EyePosition 1',0,0,0,M1.EyePosition);
+  // Test 78
+  M1:=TMatrix3D.CreateRotationZ(Pi/3)*TMatrix3D.CreateTranslation(Point3D(1,0,0));
+  AssertPoint3D('EyePosition 2',-0.50, 0.8660, 0,M1.EyePosition);
+end;
+
+procedure TTestMatrix3D.TestInverse;
+begin
+  // Test 77
+  M1:=TMatrix3D.CreateRotationZ(Pi/3);
+  M2:=TMatrix3D.CreateRotationZ(-Pi/3);
+  AssertMatrix3D('Inverse 1',M2,M1.Inverse);
+  AssertMatrix3D('Inverse 2',TMatrix3D.Identity,M1*M1.Inverse);
+
+end;
+
+procedure TTestMatrix3D.TestScale;
+begin
+  M1:=TMatrix3D.Create( 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16);
+  M3:=M1.Scale(2);
+  AssertMatrix3D('Multiplied',[ 2, 4, 6, 8,
+                                10, 12, 14, 16,
+                                18, 20, 22, 24,
+                                26, 28, 30, 32],M3);
+
+end;
+
+procedure TTestMatrix3D.TestToMatrix;
+
+Var
+  M : TMatrix;
+
+begin
+  M1:=TMatrix3D.CreateRotationZ(Pi/3);
+  M:=M1.ToMatrix;
+  AssertMatrix('To Matrix',TMatrix.CreateRotation(Pi/3),M);
+end;
+
+procedure TTestMatrix3D.TestTranspose;
+begin
+  M1:=TMatrix3D.Create( 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16);
+  M2:=TMatrix3D.Create( [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16]);
+  M3:=M1.Transpose;
+  AssertMatrix3D('Transposed',M2,M3);
+end;
+
+initialization
+   RegisterTest(TTestMatrix3D);
+end.
+

packages/rtl-objpas/tests/utcpoint.pp

@@ -0,0 +1,359 @@
+{
+    This file is part of the Free Pascal run time library.
+    Copyright (c) 2023 by Michael Van Canneyt
+    member of the Free Pascal development team
+
+    See the file COPYING.FPC, included in this distribution,
+    for details about the copyright.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+ **********************************************************************}
+
+unit utcpoint;
+
+{$mode ObjFPC}{$H+}
+
+interface
+
+uses
+  Classes, SysUtils, fpcunit, testregistry, types, utmathvectorbase, system.math.vectors;
+
+type
+    { TTestPoint3D }
+
+  TTestPoint3D = Class(TCMathVectorsBase)
+  Private
+    FP : Array[1..3] of TPoint3D;
+    procedure ClearPoints;
+    function GetP(aIndex: Integer): TPoint3D;
+    procedure SetP(aIndex: Integer; aValue: TPoint3D);
+  Protected
+    procedure SetUp; override;
+    procedure TearDown; override;
+    Property P1 : TPoint3D Index 1 Read GetP Write SetP;
+    Property P2 : TPoint3D Index 2 Read GetP Write SetP;
+    Property P3 : TPoint3D Index 3 Read GetP Write SetP;
+  Published
+    procedure TestHookUp;
+    Procedure TestZero;
+    Procedure TestCreate;
+    Procedure TestCreate2DPointAndZ;
+    Procedure TestCreateVector3D;
+    Procedure TestAdd;
+    Procedure TestMultiplyFactor;
+    Procedure TestMultiply;
+    Procedure TestSubtract;
+    Procedure TestEqual;
+    Procedure TestDivide;
+    Procedure TestUnEqual;
+    Procedure TestUnaryMinus;
+    Procedure TestAngleCosine;
+    Procedure TestCrossProduct;
+    Procedure TestDistance;
+    Procedure TestDotProduct;
+    Procedure TestEqualsTo;
+    Procedure TestLength;
+    Procedure TestMidPoint;
+    Procedure TestNormalize;
+    Procedure TestOffsetPoint;
+    Procedure TestOffsetDelta;
+    Procedure TestReflect;
+    Procedure TestRotate;
+  end;
+
+
+implementation
+{ TTestPoint3D }
+
+procedure TTestPoint3D.ClearPoints;
+
+var
+  I : integer;
+
+begin
+  For I:=1 to 3 do
+    begin
+    FP[I].X:=0;
+    FP[I].Y:=0;
+    FP[I].Z:=0;
+    end;
+end;
+
+function TTestPoint3D.GetP(aIndex: Integer): TPoint3D;
+begin
+  Result.X:=FP[AIndex].X;
+  Result.Y:=FP[AIndex].Y;
+  Result.Z:=FP[AIndex].Z;
+end;
+
+procedure TTestPoint3D.SetP(aIndex: Integer; aValue: TPoint3D);
+begin
+  FP[AIndex].X:=aValue.X;
+  FP[AIndex].Y:=aValue.Y;
+  FP[AIndex].Z:=aValue.Z;
+end;
+
+procedure TTestPoint3D.SetUp;
+begin
+  inherited SetUp;
+  ClearPoints;
+end;
+
+procedure TTestPoint3D.TearDown;
+begin
+  inherited TearDown;
+end;
+
+procedure TTestPoint3D.TestHookUp;
+var
+  I : Integer;
+begin
+  For I:=1 to 3 do
+    begin
+    AssertEquals('Point3D '+intTostr(i)+'.X',0.0,FP[I].X);
+    AssertEquals('Point3D '+intTostr(i)+'.Y',0.0,FP[I].Y);
+    AssertEquals('Point3D '+intTostr(i)+'.Z',0.0,FP[I].Z);
+    end;
+end;
+
+procedure TTestPoint3D.TestZero;
+begin
+  P1:=TPoint3D.Zero;
+  AssertPoint3D('Zero',0,0,0,P1);
+end;
+
+procedure TTestPoint3D.TestCreate;
+begin
+  P1:=TPoint3D.Create(1,2,3);
+  AssertPoint3D('Create',1,2,3,P1);
+end;
+
+procedure TTestPoint3D.TestCreate2DPointAndZ;
+begin
+  P1:=TPoint3D.Create(PointF(1,2),3);
+  AssertPoint3D('Create',1,2,3,P1);
+end;
+
+procedure TTestPoint3D.TestCreateVector3D;
+begin
+  P1:=TPoint3D.Create([1,2,3,4]);
+  AssertPoint3D('Create',1,2,3,P1);
+end;
+
+procedure TTestPoint3D.TestAdd;
+begin
+  P1:=TPoint3D.Create(1,2,3);
+  P2:=TPoint3D.Create(5,4,3);
+  P3:=P2+P1;
+  AssertPoint3D('Add 1',6,6,6,P3);
+  P2:=TPoint3D.Create(7,8,9);
+  P3:=P2+P1;
+  AssertPoint3D('Add 2',8,10,12,P3);
+end;
+
+procedure TTestPoint3D.TestMultiplyFactor;
+begin
+  P1:=TPoint3D.Create(1,2,3);
+  P2:=P1*3;
+  AssertPoint3D('Factor 1',3,6,9,P2);
+  P2:=2*P1;
+  AssertPoint3D('Factor 2',2,4,6,P2);
+end;
+
+procedure TTestPoint3D.TestMultiply;
+begin
+  P1:=TPoint3D.Create(1,2,3);
+  P2:=TPoint3D.Create(5,4,3);
+  P3:=P2*P1;
+  AssertPoint3D('Create',5,8,9,P3);
+end;
+
+procedure TTestPoint3D.TestSubtract;
+begin
+  P1:=TPoint3D.Create(1,2,3);
+  P2:=TPoint3D.Create(5,4,3);
+  P3:=P2-P1;
+  AssertPoint3D('Subtract 1',4,2,0,P3);
+  P2:=TPoint3D.Create(7,8,9);
+  P3:=P2-P1;
+  AssertPoint3D('Subtract 2',6,6,6,P3);
+end;
+
+procedure TTestPoint3D.TestEqual;
+begin
+  P1:=TPoint3D.Create(1,2,3);
+  P2:=TPoint3D.Create(5,4,3);
+  AssertFalse('Equal NOk',P1=P2);
+  P2:=TPoint3D.Create(1,2,3);
+  AssertTrue('Equal OK',P1=P2);
+end;
+
+procedure TTestPoint3D.TestDivide;
+begin
+  P1:=TPoint3D.Create(3,6,9);
+  P2:=P1/3;
+  AssertPoint3D('Factor 1',1,2,3,P2);
+end;
+
+procedure TTestPoint3D.TestUnEqual;
+begin
+  P1:=TPoint3D.Create(1,2,3);
+  P2:=TPoint3D.Create(5,4,3);
+  AssertTrue('unEqual Ok',P1<>P2);
+  P2:=TPoint3D.Create(1,2,3);
+  AssertFalse('unEqual OK',P1<>P2);
+end;
+
+procedure TTestPoint3D.TestUnaryMinus;
+begin
+  P1:=TPoint3D.Create(1,2,3);
+  P2:=-P1;
+  AssertPoint3D('Unary minus',-1,-2,-3,P2);
+end;
+
+procedure TTestPoint3D.TestAngleCosine;
+begin
+  P1:=TPoint3D.Create(1,0,0);
+  P2:=TPoint3D.Create(0,1,0);
+  AssertEquals('Angle',0, P1.AngleCosine(P2));
+  P1:=TPoint3D.Create(1,0,0);
+  P2:=TPoint3D.Create(0,0,1);
+  AssertEquals('Angle',0, P1.AngleCosine(P2));
+  P1:=TPoint3D.Create(0,1,0);
+  P2:=TPoint3D.Create(0,0,1);
+  AssertEquals('Angle',0, P1.AngleCosine(P2));
+  P1:=TPoint3D.Create(0,1,0);
+  P2:=TPoint3D.Create(1,1,0);
+  AssertEquals('Angle',c45, P1.AngleCosine(P2));
+end;
+
+procedure TTestPoint3D.TestCrossProduct;
+begin
+  P1:=TPoint3D.Create(1,1,0);
+  P2:=TPoint3D.Create(2,2,0);
+  P3:=P2.CrossProduct(P1);
+  AssertPoint3D('T1',0,0,0,P3);
+  P1:=TPoint3D.Create(1,1,0);
+  P2:=TPoint3D.Create(2,2,0);
+  P3:=P2.CrossProduct(P1);
+  AssertPoint3D('T2',0,0,0,P3);
+end;
+
+procedure TTestPoint3D.TestDistance;
+begin
+  P1:=TPoint3D.Create(1,1,0);
+  P2:=TPoint3D.Create(1,2,0);
+  AssertEquals('Dist 1,',1,P1.Distance(P2));
+  P1:=TPoint3D.Create(0,1,1);
+  P2:=TPoint3D.Create(0,2,1);
+  AssertEquals('Dist 2,',1,P1.Distance(P2));
+  P1:=TPoint3D.Create(0,0,0);
+  P2:=TPoint3D.Create(1,1,0);
+  AssertEquals('Dist 3,',c45*2,P1.Distance(P2));
+end;
+
+procedure TTestPoint3D.TestDotProduct;
+begin
+  P1:=TPoint3D.Create(3,4,9);
+  P2:=TPoint3D.Create(3,4,9);
+  AssertEquals('Test 1',9+16+81,P1.DotProduct(P2));
+  P2:=TPoint3D.Create(1,1,0);
+  AssertEquals('Test 1',7,P1.DotProduct(P2));
+end;
+
+procedure TTestPoint3D.TestEqualsTo;
+begin
+  P1:=TPoint3D.Create(3,4,9);
+  P2:=TPoint3D.Create(3,4,9);
+  AssertTrue('Test 1',P1.EqualsTo(P2));
+  P2:=TPoint3D.Create(1,1,1);
+  AssertFalse('Test 2',P1.EqualsTo(P2));
+end;
+
+procedure TTestPoint3D.TestLength;
+begin
+  P1:=TPoint3D.Create(1,0,0);
+  AssertEquals('Dist 1,',1,P1.Length);
+  P1:=TPoint3D.Create(0,1,1);
+  AssertEquals('Dist 2,',c45*2,P1.Length);
+  P1:=TPoint3D.Create(1,1,1);
+  AssertEquals('Dist 3,',1.7321,P1.Length);
+end;
+
+procedure TTestPoint3D.TestMidPoint;
+begin
+  P1:=TPoint3D.Create(0,0,1);
+  P2:=TPoint3D.Create(0,2,1);
+  AssertPoint3D('Midpoint1',0,1,1,P1.MidPoint(P2));
+  P1:=TPoint3D.Create(1,0,0);
+  P2:=TPoint3D.Create(0,1,0);
+  AssertPoint3D('Midpoint2',0.5,0.5,0,P1.MidPoint(P2));
+end;
+
+procedure TTestPoint3D.TestNormalize;
+begin
+  P1:=TPoint3D.Create(1,1,1);
+  P2:=P1.Normalize;
+  AssertPoint3D('Normalize 1',0.5773,0.5773,0.5773,P2);
+  P1:=TPoint3D.Create(0,0,1);
+  P2:=P1.Normalize;
+  AssertPoint3D('Normalize 2',0,0,1,P2);
+  P1:=TPoint3D.Create(1,0,0);
+  P2:=P1.Normalize;
+  AssertPoint3D('Normalize 3',1,0,0,P2);
+end;
+
+procedure TTestPoint3D.TestOffsetPoint;
+
+Var
+  P : TPoint3D; // Cannot use property
+
+begin
+  P:=TPoint3D.Create(1,2,3);
+  P.Offset(Point3D(4,5,6));
+  AssertPoint3D('Off',5,7,9,P);
+end;
+
+procedure TTestPoint3D.TestOffsetDelta;
+
+Var
+  P : TPoint3D; // Cannot use property
+
+begin
+  P:=TPoint3D.Create(0,0,0);
+  P.Offset(1,2,3);
+  AssertPoint3D('Off',1,2,3,P);
+end;
+
+procedure TTestPoint3D.TestReflect;
+
+// P2 is normal of reflecting surface: X,Y plane.
+
+begin
+  P1:=TPoint3D.Create(0,0,1);
+  P2:=TPoint3D.Create(0,0,1);
+  P3:=P1.Reflect(P2);
+  AssertPoint3D('Reflect 1',0,0,-1,P3);
+  P1:=TPoint3D.Create(0,-1,-1);
+  P3:=P1.Reflect(P2);
+  AssertPoint3D('Reflect 2',0,-1,1,P3);
+end;
+
+procedure TTestPoint3D.TestRotate;
+begin
+  P1:=TPoint3D.Create(1,0,0);
+  P2:=TPoint3D.Create(0,0,1); // Z-axis
+  P3:=P1.Rotate(P2,Pi/2);
+  AssertPoint3D('Rotate 1',0,1,0,P3);
+  P2:=TPoint3D.Create(0,1,0); // Y-axis
+  P3:=P1.Rotate(P2,Pi/2);
+  AssertPoint3D('Rotate 2',0,0,-1,P3); // Z is pointing down!
+end;
+
+initialization
+  RegisterTests([TTestPoint3D]);
+end.
+

packages/rtl-objpas/tests/utcquaternion.pp

@@ -0,0 +1,177 @@
+{
+    This file is part of the Free Pascal run time library.
+    Copyright (c) 2023 by Michael Van Canneyt
+    member of the Free Pascal development team
+
+    See the file COPYING.FPC, included in this distribution,
+    for details about the copyright.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+ **********************************************************************}
+
+unit utcquaternion;
+
+{$mode ObjFPC}{$H+}
+
+interface
+
+uses
+  Classes, SysUtils, fpcunit, testregistry, utmathvectorbase, types, system.math.vectors;
+
+Type
+  { TTestQuaternion3D }
+
+  TTestQuaternion3D = class(TCMathVectorsBase)
+  Private
+    FQ : Array[1..3] of TQuaternion3D;
+    procedure Clear;
+    function GetQ(AIndex: Integer): TQuaternion3D;
+    procedure SetQ(AIndex: Integer; AValue: TQuaternion3D);
+  Protected
+    procedure SetUp; override;
+    procedure TearDown; override;
+    Property Q1 : TQuaternion3D Index 1 Read GetQ Write SetQ;
+    Property Q2 : TQuaternion3D Index 2 Read GetQ Write SetQ;
+    Property Q3 : TQuaternion3D Index 3 Read GetQ Write SetQ;
+  Published
+    procedure TestHookUp;
+    procedure TestCreateVectorAngle;
+    procedure TestCreateVectorVector;
+    procedure TestCreateYawPitchRoll;
+    procedure TestCreateMatrix;
+    procedure TestAssignToMatrix;
+    procedure TestMultiply;
+  end;
+
+implementation
+
+{ TTestQuaternion3D }
+
+procedure TTestQuaternion3D.Clear;
+begin
+
+end;
+
+function TTestQuaternion3D.GetQ(AIndex: Integer): TQuaternion3D;
+begin
+  Result:=FQ[aIndex];
+end;
+
+procedure TTestQuaternion3D.SetQ(AIndex: Integer; AValue: TQuaternion3D);
+begin
+  FQ[aIndex]:=aValue;
+end;
+
+procedure TTestQuaternion3D.SetUp;
+begin
+  inherited SetUp;
+  Clear;
+end;
+
+procedure TTestQuaternion3D.TearDown;
+begin
+  Clear;
+  inherited TearDown;
+end;
+
+procedure TTestQuaternion3D.TestHookUp;
+
+var
+  I,J : Integer;
+
+begin
+  For I:=1 to 3 do
+    begin
+    AssertEquals(Format('Q%d.Real',[I]),0.0,FQ[I].RealPart);
+    For J:=0 to 2 do
+      AssertEquals(Format('Q%d.ImagPart[%d]',[I,J]),0.0,FQ[I].ImagPart.V[J]);
+    end;
+end;
+
+Procedure TTestQuaternion3D.TestCreateVectorAngle;
+begin
+  // Test 202
+  Q1:=TQuaternion3D.Create(Point3D(1,0,0),Pi/2);
+  AssertQuaternion('Create 2', 0.7071, 0.7071,  0.0000, 0.0000, Q1);
+  // Test 203
+  Q1:=TQuaternion3D.Create(Point3D(0,1,0),Pi/2);
+  AssertQuaternion('Create 3', 0.7071, 0.0000, 0.7071, 0.0000, Q1);
+  // Test 204
+  Q1:=TQuaternion3D.Create(Point3D(0,0,1),Pi/2);
+  AssertQuaternion('Create 4', 0.7071, 0.0000, 0.0000, 0.7071, Q1);
+end;
+
+Procedure TTestQuaternion3D.TestCreateVectorVector;
+
+begin
+  // Test 201
+  Q1:=TQuaternion3D.Create(Point3D(1,0,0),Point3D(0,1,0));
+  AssertQuaternion('Create 1', 0.7071,  0.0000,  0.0000, 1.0000, Q1);
+end;
+
+Procedure TTestQuaternion3D.TestCreateYawPitchRoll;
+
+begin
+  // Test 205
+  Q1:=TQuaternion3D.Create(Pi/2,0,0);
+  AssertQuaternion('Create 1', 0.7071,  0.0000, 0.7071, 0.0000, Q1);
+  // Test 206
+  Q1:=TQuaternion3D.Create(0,Pi/2,0);
+  AssertQuaternion('Create 2', 0.7071,  0.7071, 0.0000, 0.0000, Q1);
+  // Test 207
+  Q1:=TQuaternion3D.Create(0,0,Pi/2);
+  AssertQuaternion('Create 3', 0.7071,  0.0000,  0.0000, 0.7071, Q1);
+end;
+
+procedure TTestQuaternion3D.TestCreateMatrix;
+
+var
+  M : TMatrix3D;
+
+begin
+  M:=TMatrix3D.CreateRotationX(pi/2);
+  Q1:=TQuaternion3D.Create(M);
+  AssertQuaternion('Create 1', 0.7071, 0.7071, 0.0000, 0.0000, Q1);
+  M:=TMatrix3D.CreateRotationY(pi/2);
+  Q1:=TQuaternion3D.Create(M);
+  AssertQuaternion('Create 2', 0.7071, 0.0000, 0.7071, 0.0000, Q1);
+  M:=TMatrix3D.CreateRotationZ(pi/2);
+  Q1:=TQuaternion3D.Create(M);
+  AssertQuaternion('Create 2', 0.7071, 0.0000, 0.0000, 0.7071, Q1);
+end;
+
+procedure TTestQuaternion3D.TestAssignToMatrix;
+
+var
+  M,M2 : TMatrix3D;
+
+begin
+  M:=TMatrix3D.CreateRotationX(pi/2);
+  Q1:=TQuaternion3D.Create(M);
+  M2:=Q1;
+  AssertMatrix3D('Assign',M,M2);
+end;
+
+procedure TTestQuaternion3D.TestMultiply;
+begin
+  // test 208
+  Q1:=TQuaternion3D.Create(0,0,Pi/2);
+  Q2:=TQuaternion3D.Create(0,Pi/2,0);
+  Q3:=Q1*Q2;
+  AssertQuaternion('Multiply 1', 0.5,0.5,0.5,0.5, Q3);
+  // test 209
+  Q1:=TQuaternion3D.Create(Point3D(1,0,0),Pi/2);
+  Q2:=TQuaternion3D.Create(Point3D(0,1,0),Pi/2);
+  Q3:=Q1*Q2;
+  AssertQuaternion('Multiply 2', 0.5,0.5,0.5,0.5, Q3);
+end;
+
+
+Initialization
+  RegisterTest(TTestQuaternion3D);
+
+end.
+

packages/rtl-objpas/tests/utcvector.pp

@@ -0,0 +1,286 @@
+{
+    This file is part of the Free Pascal run time library.
+    Copyright (c) 2023 by Michael Van Canneyt
+    member of the Free Pascal development team
+
+    See the file COPYING.FPC, included in this distribution,
+    for details about the copyright.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+ **********************************************************************}
+unit utcvector;
+
+{$mode ObjFPC}{$H+}
+
+interface
+
+uses
+  Classes, SysUtils, fpcunit, testregistry, types, utmathvectorbase, system.math.vectors;
+
+Type
+
+  { TTestVector }
+
+  TTestVector = class(TCMathVectorsBase)
+  Private
+    FV : Array[1..3] of TVector;
+    procedure ClearVectors;
+    function GetV(AIndex: Integer): TVector;
+    procedure SetV(AIndex: Integer; AValue: TVector);
+  Protected
+    procedure SetUp; override;
+    procedure TearDown; override;
+    Property V1 : TVector Index 1 Read GetV Write SetV;
+    Property V2 : TVector Index 2 Read GetV Write SetV;
+    Property V3 : TVector Index 3 Read GetV Write SetV;
+  Published
+    procedure TestHookUp;
+    Procedure TestZero;
+    Procedure TestCreate;
+    Procedure TestCreateW;
+    Procedure TestAssign;
+    Procedure TestAssignPointf;
+    Procedure TestAssignToPointf;
+    Procedure TestAdd;
+    Procedure TestMultiplyFactor;
+    Procedure TestDivide;
+    Procedure TestEqual;
+    Procedure TestNotEqual;
+    Procedure TestSubtract;
+    Procedure TestLength;
+    Procedure TestNormalize;
+    Procedure TestCrossProduct;
+    Procedure TestDotProduct;
+    Procedure TestToPointF;
+  end;
+
+implementation
+{ TTestVector }
+
+function TTestVector.GetV(AIndex: Integer): TVector;
+begin
+  Result:=FV[aIndex];
+end;
+
+procedure TTestVector.SetV(AIndex: Integer; AValue: TVector);
+begin
+  FV[aIndex]:=aValue;
+end;
+
+procedure TTestVector.ClearVectors;
+
+var
+  I : integer;
+
+begin
+  For I:=1 to 3 do
+    begin
+    FV[I].X:=0;
+    FV[I].Y:=0;
+    FV[I].W:=0;
+    end;
+end;
+
+
+
+procedure TTestVector.SetUp;
+begin
+  inherited SetUp;
+  ClearVectors;
+end;
+
+procedure TTestVector.TearDown;
+begin
+  inherited TearDown;
+end;
+
+procedure TTestVector.TestHookUp;
+var
+  I : Integer;
+begin
+  For I:=1 to 3 do
+    begin
+    AssertEquals('Vector '+intTostr(i)+'.X',0.0,FV[I].X);
+    AssertEquals('Vector '+intTostr(i)+'.Y',0.0,FV[I].Y);
+    AssertEquals('Vector '+intTostr(i)+'.W',0.0,FV[I].W);
+    end;
+end;
+
+procedure TTestVector.TestZero;
+begin
+  V1:=TVector.Zero;
+  AssertEquals('Vector.X',0.0,V1.X);
+  AssertEquals('Vector.Y',0.0,V1.Y);
+  AssertEquals('Vector.W',0.0,V1.W);
+end;
+
+procedure TTestVector.TestCreate;
+begin
+  V1:=TVector.Create(1,2);
+  AssertEquals('Vector.X',1.0,V1.X);
+  AssertEquals('Vector.Y',2.0,V1.Y);
+  AssertEquals('Vector.W',DefaultVectorWidth,V1.W);
+end;
+
+procedure TTestVector.TestCreateW;
+begin
+  V1:=TVector.Create(1,2,3);
+  AssertVector('Vector',1,2,3,V1);
+end;
+
+procedure TTestVector.TestAssign;
+begin
+  V2:=TVector.Create(1,2,3);
+  V1:=V2;
+  AssertVector('Assign',1,2,3,V1);
+end;
+
+procedure TTestVector.TestAssignPointf;
+
+var
+  P : TPointF;
+
+begin
+  P:=PointF(1,2);
+  V1:=P;
+  AssertVector('Vector',1,2,DefaultVectorWidth,V1);
+end;
+
+procedure TTestVector.TestAssignToPointf;
+Var
+  P : TPointF;
+
+begin
+  V1:=TVector.Create(1,2,3);
+  P:=V1;
+  AssertEquals('Assign 1',PointF(0.3333,0.6666),P);
+  V1:=TVector.Create(1,2,1);
+  P:=V1;
+  AssertEquals('Assign 2',PointF(1,2),P);
+  V1:=TVector.Create(1,2,0);
+  P:=V1;
+  AssertEquals('Assign 3',PointF(1,2),P);
+end;
+
+procedure TTestVector.TestAdd;
+
+begin
+  V1:=TVector.Create(1,2,3);
+  V2:=TVector.Create(6,5,4);
+  V3:=V1+V2;
+  AssertVector('Vector',7,7,7,V3);
+end;
+
+procedure TTestVector.TestMultiplyFactor;
+begin
+  V1:=TVector.Create(1,2,3);
+  V2:=V1*3;
+  AssertVector('Vector 1',3,6,9,V2);
+  V2:=3*V1;
+  AssertVector('Vector 2',3,6,9,V2);
+end;
+
+procedure TTestVector.TestDivide;
+begin
+  V1:=TVector.Create(1,2,3);
+  V2:=V1/3;
+  AssertVector('Vector 1',0.3333,0.6666,1,V2);
+end;
+
+procedure TTestVector.TestEqual;
+begin
+  V1:=TVector.Create(1,2,3);
+  V2:=TVector.Create(1,2,3);
+  AssertTrue('Equal 1',V1=V2);
+  V2:=TVector.Create(3,2,1);
+  AssertFalse('Equal 2',V1=V2);
+  V2:=TVector.Create(1+TEpsilon.Vector*0.99,2,3);
+  AssertTrue('Equal within precision',V1=V2);
+  V2:=TVector.Create(1+TEpsilon.Vector*1.1,2,3);
+  AssertFalse('Unequal outside precision',V1=V2);
+end;
+
+procedure TTestVector.TestNotEqual;
+begin
+  V1:=TVector.Create(1,2,3);
+  V2:=TVector.Create(1,2,3);
+  AssertFalse('Not Equal',V1<>V2);
+  V2:=TVector.Create(3,2,1);
+  AssertTrue('Equal',V1<>V2);
+  V2:=TVector.Create(1+TEpsilon.Vector*0.99,2,3);
+  AssertFalse('Equal within precision',V1<>V2);
+  V2:=TVector.Create(1+TEpsilon.Vector*1.1,2,3);
+  AssertTrue('Unequal outside precision',V1<>V2);
+end;
+
+procedure TTestVector.TestSubtract;
+begin
+  V1:=TVector.Create(1,2,3);
+  V2:=TVector.Create(6,5,4);
+  V3:=V2-V1;
+  AssertVector('Vector',5,3,1,V3);
+end;
+
+procedure TTestVector.TestLength;
+begin
+  V1:=TVector.Create(3,4,0);
+  AssertEquals('Length 1',5,V1.Length);
+  V1:=TVector.Create(3,4,1);
+  AssertEquals('Length 1',Sqrt(26),V1.Length,TEpsilon.Vector);
+end;
+
+procedure TTestVector.TestNormalize;
+begin
+  V1:=TVector.Create(3,4,0);
+  V2:=V1.Normalize;
+  AssertVector('No width',3/5,4/5,0,V2);
+  AssertEquals('Length 1',1,V2.Length,TEpsilon.Vector);
+  V1:=TVector.Create(3,4,1);
+  V2:=V1.Normalize;
+  AssertVector('No width',3/Sqrt(26),4/Sqrt(26),1/Sqrt(26),V2);
+  AssertEquals('Length 1',1,V2.Length,TEpsilon.Vector);
+end;
+
+procedure TTestVector.TestCrossProduct;
+begin
+  V1:=TVector.Create(1,1,0);
+  V2:=TVector.Create(2,2,0);
+  V3:=V2.CrossProduct(V1);
+  AssertVector('T1',0,0,0,V3);
+  V1:=TVector.Create(1,1,0);
+  V2:=TVector.Create(2,2,0);
+  V3:=V2.CrossProduct(V1);
+  AssertVector('T1',0,0,0,V3);
+end;
+
+procedure TTestVector.TestDotProduct;
+begin
+  V1:=TVector.Create(3,4,9);
+  V2:=TVector.Create(3,4,9);
+  AssertEquals('Test 1',9+16+81,V1.DotProduct(V2));
+  V2:=TVector.Create(1,1,0);
+  AssertEquals('Test 1',7,V1.DotProduct(V2));
+end;
+
+procedure TTestVector.TestToPointF;
+
+var
+  P : TPointF;
+
+begin
+  V1:=TVector.Create(1,2,3);
+  P:=V1.ToPointF;
+  AssertEquals('ToPointF 1',PointF(0.3333,0.6666),P);
+  V1:=TVector.Create(1,2,1);
+  P:=V1.ToPointF;
+  AssertEquals('ToPointF 2',PointF(1,2),P);
+end;
+
+
+initialization
+  RegisterTest(TTestVector);
+end.
+

packages/rtl-objpas/tests/utmathvectorbase.pas

@@ -0,0 +1,212 @@
+unit utmathvectorbase;
+
+{$mode objfpc}{$H+}
+
+interface
+
+uses
+  Classes, SysUtils, fpcunit, testutils, testregistry, Types, system.math.vectors;
+
+const
+  c45 = Sqrt(2)/2;  // cosine/sine 45°
+  c60 = Sqrt(3)/2;  // cosine 60°
+  c30 = 1/2;        // cosine 30°
+
+  s30 = C60; // sine 30°
+  s45 = c45; // sine 45°
+  s60 = c30; // sine 60°
+
+
+
+type
+  { TCMathVectorsBase }
+
+  TCMathVectorsBase = class(TTestCase)
+  protected
+    procedure SetUp; override;
+    procedure TearDown; override;
+    class Procedure AssertEquals(const Msg : String; aExpected,aActual : TVector); overload;
+    class Procedure AssertEquals(const Msg : String; aExpected,aActual : TPoint3D); overload;
+    class Procedure AssertEquals(const Msg : String; aExpected,aActual : TPointF); overload;
+    class Procedure AssertEquals(const Msg : String; aExpected,aActual : TVector3D); overload;
+    class Procedure AssertVector3D(const Msg : String; aExpectedX,aExpectedY,aExpectedZ,aExpectedW: Single; aActual : TVector3D);
+    class Procedure AssertVector(const Msg : String; aExpectedX,aExpectedY,aExpectedW: Single; aActual : TVector);
+    class Procedure AssertPoint3D(const Msg : String; aExpectedX,aExpectedY,aExpectedZ: Single; aActual : TPoint3D);
+    class procedure AssertMatrix(Const Msg : String; m11, m12, m13,  m21, m22, m23,  m31, m32, m33 : Single; aActual : TMatrix);
+    class procedure AssertMatrix(Const Msg : String; aExpected : Array of single; aActual : TMatrix);
+    class procedure AssertMatrix(Const Msg : String; aExpected, aActual : TMatrix);
+    class procedure AssertMatrix3D(Const Msg : String; m11, m12, m13, m14, m21, m22, m23, m24, m31, m32, m33, m34, m41, m42, m43, m44 : Single; aActual : TMatrix3D);
+    class procedure AssertMatrix3D(Const Msg : String; aExpected : Array of single; aActual : TMatrix3D);
+    class procedure AssertMatrix3D(Const Msg : String; aExpected, aActual : TMatrix3D);
+    class Procedure AssertPointF(const Msg : String; aExpectedX,aExpectedY: Single; aActual : TPointF);
+    class Procedure AssertQuaternion(const Msg : String; aExpectedReal,aExpectedImagX,aExpectedImagY,aExpectedImagZ: Single; aActual : TQuaternion3D);
+    class Procedure AssertQuaternion(const Msg : String; aExpectedReal: Single; aExpectedImag : TPoint3D; aActual : TQuaternion3D);
+  end;
+
+
+
+implementation
+
+procedure TCMathVectorsBase.SetUp;
+begin
+
+end;
+
+procedure TCMathVectorsBase.TearDown;
+begin
+
+end;
+
+class procedure TCMathVectorsBase.AssertEquals(const Msg: String; aExpected,
+  aActual: TVector);
+begin
+  AssertEquals(Msg+' X',aExpected.X,aActual.X,TEpsilon.Vector);
+  AssertEquals(Msg+' Y',aExpected.Y,aActual.Y,TEpsilon.Vector);
+  AssertEquals(Msg+' W',aExpected.W,aActual.W,TEpsilon.Vector);
+end;
+
+class procedure TCMathVectorsBase.AssertEquals(const Msg: String; aExpected,
+  aActual: TPoint3D);
+begin
+  AssertEquals(Msg+' X',aExpected.X,aActual.X,TEpsilon.Vector);
+  AssertEquals(Msg+' Y',aExpected.Y,aActual.Y,TEpsilon.Vector);
+  AssertEquals(Msg+' Z',aExpected.Z,aActual.Z,TEpsilon.Vector);
+end;
+
+class procedure TCMathVectorsBase.AssertEquals(const Msg: String; aExpected,
+  aActual: TPointF);
+begin
+  AssertEquals(Msg+' X',aExpected.X,aActual.X,TEpsilon.Vector);
+  AssertEquals(Msg+' Y',aExpected.Y,aActual.Y,TEpsilon.Vector);
+
+end;
+
+class procedure TCMathVectorsBase.AssertEquals(const Msg: String; aExpected,
+  aActual: TVector3D);
+begin
+  AssertEquals(Msg+' X',aExpected.X,aActual.X,TEpsilon.Vector);
+  AssertEquals(Msg+' Y',aExpected.Y,aActual.Y,TEpsilon.Vector);
+  AssertEquals(Msg+' Z',aExpected.Z,aActual.Z,TEpsilon.Vector);
+  AssertEquals(Msg+' W',aExpected.W,aActual.W,TEpsilon.Vector);
+end;
+
+class procedure TCMathVectorsBase.AssertVector3D(const Msg: String; aExpectedX,
+  aExpectedY, aExpectedZ, aExpectedW: Single; aActual: TVector3D);
+begin
+  AssertEquals(Msg+' X',aExpectedX,aActual.X,TEpsilon.Vector);
+  AssertEquals(Msg+' Y',aExpectedY,aActual.Y,TEpsilon.Vector);
+  AssertEquals(Msg+' Z',aExpectedZ,aActual.Z,TEpsilon.Vector);
+  AssertEquals(Msg+' W',aExpectedW,aActual.W,TEpsilon.Vector);
+end;
+
+class procedure TCMathVectorsBase.AssertVector(const Msg: String; aExpectedX,
+  aExpectedY, aExpectedW: Single; aActual: TVector);
+begin
+  AssertEquals(Msg+' X',aExpectedX,aActual.X,TEpsilon.Vector);
+  AssertEquals(Msg+' Y',aExpectedY,aActual.Y,TEpsilon.Vector);
+  AssertEquals(Msg+' W',aExpectedW,aActual.W,TEpsilon.Vector);
+end;
+
+class procedure TCMathVectorsBase.AssertPoint3D(const Msg: String; aExpectedX, aExpectedY, aExpectedZ: Single; aActual: TPoint3D);
+begin
+  AssertEquals(Msg+' X',aExpectedX,aActual.X,TEpsilon.Vector);
+  AssertEquals(Msg+' Y',aExpectedY,aActual.Y,TEpsilon.Vector);
+  AssertEquals(Msg+' Z',aExpectedZ,aActual.Z,TEpsilon.Vector);
+end;
+
+class procedure TCMathVectorsBase.AssertMatrix(const Msg: String; m11, m12, m13, m21, m22, m23, m31, m32, m33: Single;
+  aActual: TMatrix);
+begin
+  AssertEquals(Msg+' m11',m11,aActual.m11,TEpsilon.Vector);
+  AssertEquals(Msg+' m12',m12,aActual.m12,TEpsilon.Vector);
+  AssertEquals(Msg+' m13',m13,aActual.m13,TEpsilon.Vector);
+
+  AssertEquals(Msg+' m21',m21,aActual.m21,TEpsilon.Vector);
+  AssertEquals(Msg+' m22',m22,aActual.m22,TEpsilon.Vector);
+  AssertEquals(Msg+' m23',m23,aActual.m23,TEpsilon.Vector);
+
+  AssertEquals(Msg+' m31',m31,aActual.m31,TEpsilon.Vector);
+  AssertEquals(Msg+' m32',m32,aActual.m32,TEpsilon.Vector);
+  AssertEquals(Msg+' m33',m33,aActual.m33,TEpsilon.Vector);
+end;
+
+class procedure TCMathVectorsBase.AssertMatrix(const Msg: String; aExpected: array of single; aActual: TMatrix);
+begin
+  AssertEquals(Msg+' number of elements',9,Length(aExpected));
+  AssertMatrix(Msg,aExpected[0],aExpected[1],aExpected[2],
+                   aExpected[3],aExpected[4],aExpected[5],
+                   aExpected[6],aExpected[7],aExpected[8],aActual);
+end;
+
+class procedure TCMathVectorsBase.AssertMatrix(const Msg: String; aExpected, aActual: TMatrix);
+begin
+  With aExpected do
+    AssertMatrix(Msg,m11, m12, m13, m21, m22, m23, m31, m32, m33,aActual);
+end;
+
+class procedure TCMathVectorsBase.AssertMatrix3D(const Msg: String; m11, m12, m13, m14, m21, m22, m23, m24, m31, m32, m33,
+  m34, m41, m42, m43, m44 : Single; aActual: TMatrix3D);
+begin
+  AssertEquals(Msg+' m11',m11,aActual.m11,TEpsilon.Vector);
+  AssertEquals(Msg+' m12',m12,aActual.m12,TEpsilon.Vector);
+  AssertEquals(Msg+' m13',m13,aActual.m13,TEpsilon.Vector);
+  AssertEquals(Msg+' m14',m14,aActual.m14,TEpsilon.Vector);
+
+  AssertEquals(Msg+' m21',m21,aActual.m21,TEpsilon.Vector);
+  AssertEquals(Msg+' m22',m22,aActual.m22,TEpsilon.Vector);
+  AssertEquals(Msg+' m23',m23,aActual.m23,TEpsilon.Vector);
+  AssertEquals(Msg+' m24',m24,aActual.m24,TEpsilon.Vector);
+
+  AssertEquals(Msg+' m31',m31,aActual.m31,TEpsilon.Vector);
+  AssertEquals(Msg+' m32',m32,aActual.m32,TEpsilon.Vector);
+  AssertEquals(Msg+' m33',m33,aActual.m33,TEpsilon.Vector);
+  AssertEquals(Msg+' m34',m34,aActual.m34,TEpsilon.Vector);
+
+  AssertEquals(Msg+' m41',m41,aActual.m41,TEpsilon.Vector);
+  AssertEquals(Msg+' m42',m42,aActual.m42,TEpsilon.Vector);
+  AssertEquals(Msg+' m43',m43,aActual.m43,TEpsilon.Vector);
+  AssertEquals(Msg+' m44',m44,aActual.m44,TEpsilon.Vector);
+
+end;
+
+class procedure TCMathVectorsBase.AssertMatrix3D(const Msg: String; aExpected: array of single; aActual: TMatrix3D);
+begin
+  AssertEquals(Msg+' number of elements',16,Length(aExpected));
+  AssertMatrix3D(Msg,aExpected[0],aExpected[1],aExpected[2],aExpected[3],
+                   aExpected[4],aExpected[5],aExpected[6],aExpected[7],
+                   aExpected[8],aExpected[9], aExpected[10], aExpected[11],
+                   aExpected[12], aExpected[13], aExpected[14], aExpected[15], aActual);
+end;
+
+class procedure TCMathVectorsBase.AssertMatrix3D(const Msg: String; aExpected, aActual: TMatrix3D);
+begin
+  With aExpected do
+    AssertMatrix3D(Msg,m11, m12, m13, m14, m21, m22, m23, m24, m31, m32, m33, m34, m41, m42, m43, m44, aActual);
+end;
+
+class procedure TCMathVectorsBase.AssertPointF(const Msg: String; aExpectedX, aExpectedY: Single; aActual: TPointF);
+begin
+  AssertEquals(Msg+' X',aExpectedX,aActual.X,TEpsilon.Vector);
+  AssertEquals(Msg+' Y',aExpectedY,aActual.Y,TEpsilon.Vector);
+end;
+
+class procedure TCMathVectorsBase.AssertQuaternion(const Msg: String; aExpectedReal, aExpectedImagX, aExpectedImagY,
+  aExpectedImagZ: Single; aActual: TQuaternion3D);
+begin
+  AssertEquals(Msg+' RealPart',aExpectedReal,aActual.RealPart,TEpsilon.Vector);
+  AssertEquals(Msg+' ImagPart.X',aExpectedImagX,aActual.ImagPart.X,TEpsilon.Vector);
+  AssertEquals(Msg+' ImagPart.Y',aExpectedImagY,aActual.ImagPart.Y,TEpsilon.Vector);
+  AssertEquals(Msg+' ImagPart.Z',aExpectedImagZ,aActual.ImagPart.Z,TEpsilon.Vector);
+end;
+
+class procedure TCMathVectorsBase.AssertQuaternion(const Msg: String; aExpectedReal: Single; aExpectedImag: TPoint3D;
+  aActual: TQuaternion3D);
+begin
+  AssertEquals(Msg+' RealPart',aExpectedReal,aActual.RealPart,TEpsilon.Vector);
+  AssertEquals(Msg+' ImagPart.X',aExpectedImag.X,aActual.ImagPart.X,TEpsilon.Vector);
+  AssertEquals(Msg+' ImagPart.Y',aExpectedImag.Y,aActual.ImagPart.Y,TEpsilon.Vector);
+  AssertEquals(Msg+' ImagPart.Z',aExpectedImag.Z,aActual.ImagPart.Z,TEpsilon.Vector);
+end;
+
+end.
+