nichkara/lace
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Add initial prototype.

Browse Source
This commit is contained in:
Rod Kay
2022-07-31 17:34:54 +10:00
commit 54a53b2ac0
1421 changed files with 358874 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

253
1-base/math/source/generic/pure/algebra/any_math-any_algebra-any_linear-any_d2.adb Normal file
View File

@@ -0,0 +1,253 @@
package body any_Math.any_Algebra.any_linear.any_d2
is
-----------
-- Vector_2
--
function Angle_between_pre_Norm (U, V : in Vector_2) return Radians
is
use Functions, Vectors;
Val : Real := U * V; -- Dot product.
begin
if val < -1.0 then val := -1.0; -- Clamp to avoid rounding errors. arcCos will
elsif val > 1.0 then val := 1.0; -- fail with values outside this range.
end if;
return arcCos (Val);
end Angle_between_pre_Norm;
function Midpoint (From, To : in Vector_2) return Vector_2
is
begin
return [(From (1) + To (1)) * 0.5,
(From (2) + To (2)) * 0.5];
end Midpoint;
function Distance (From, To : in Vector_2) return Real
is
begin
return abs (From - To);
end Distance;
function Interpolated (From, To : in Vector_2;
Percent : in unit_Percentage) return Vector_2
is
P : constant Real := to_Real (Percent);
S : constant Real := 1.0 - P;
begin
return [S * From (1) + P * To (1),
S * From (2) + P * To (2)];
end Interpolated;
-------------
-- Matrix_2x2
--
function to_Matrix (Row_1, Row_2 : in Vector_2) return Matrix_2x2
is
begin
return [[Row_1 (1), Row_1 (2)],
[Row_2 (1), Row_2 (2)]];
end to_Matrix;
function to_rotation_Matrix (Angle : in Radians ) return Matrix_2x2
is
use Functions;
begin
return [[ cos (Angle), sin (Angle)],
[-sin (Angle), cos (Angle)]];
end to_rotation_Matrix;
function up_Direction (Self : in Matrix_2x2) return Vector_2
is
begin
return Normalised (Row (Self, 2));
end up_Direction;
function right_Direction (Self : in Matrix_2x2) return Vector_2
is
begin
return Normalised (Row (Self, 1));
end right_Direction;
function to_rotation_Transform (Rotation : in Matrix_2x2) return Matrix_3x3
is
begin
return [[Rotation (1, 1), Rotation (1, 2), 0.0],
[Rotation (2, 1), Rotation (2, 2), 0.0],
[ 0.0, 0.0, 1.0]];
end to_rotation_Transform;
-------------
-- Transform
--
function to_Transform_2d (From : in Matrix_3x3) return Transform_2d
is
begin
return (Rotation => get_Rotation (From),
Translation => get_Translation (From));
end to_Transform_2d;
function to_Transform (From : in Transform_2d) return Matrix_3x3
is
begin
return to_rotation_Transform (From.Rotation) * to_translation_Transform (From.Translation);
end to_Transform;
function to_translation_Transform (Translation : Vector_2) return Matrix_3x3
is
begin
return [[ 1.0, 0.0, 0.0],
[ 0.0, 1.0, 0.0],
[Translation (1), Translation (2), 1.0]];
end to_translation_Transform;
function to_rotation_Transform (Angle : in Radians) return Matrix_3x3
is
use Functions;
begin
return [[ cos (Angle), sin (Angle), 0.0],
[-sin (Angle), cos (Angle), 0.0],
[ 0.0, 0.0, 1.0]];
end to_rotation_Transform;
function to_scale_Transform (Scale : in Vector_2) return Matrix_3x3
is
begin
return [[Scale (1), 0.0, 0.0],
[ 0.0, Scale (2), 0.0],
[ 0.0, 0.0, 1.0]];
end to_scale_Transform;
function to_Transform (Rotation : in Matrix_2x2;
Translation : in Vector_2) return Matrix_3x3
is
begin
return [[Rotation (1, 1), Rotation (1, 2), 0.0],
[Rotation (2, 1), Rotation (2, 2), 0.0],
[Translation (1), Translation (2), 1.0]];
end to_Transform;
function to_Transform_2d (Rotation : in Radians;
Translation : in Vector_2) return Transform_2d
is
begin
return (to_rotation_Matrix (Rotation),
Translation);
end to_Transform_2d;
function "*" (Left : in Vector_2; Right : in Transform_2d) return Vector_2
is
Pad : constant Vector_3 := [Left (1), Left (2), 1.0];
Result : constant Vector_3 := Pad * to_Transform (Right);
begin
return Vector_2 (Result (1 .. 2));
end "*";
function "*" (Left : in Vector_2; Right : in Matrix_3x3) return Vector_2
is
use Vectors;
Result : constant Vector := Vector (Left & 1.0) * Matrix (Right);
begin
return Vector_2 (Result (1 .. 2));
end "*";
function Invert (Transform : in Transform_2d) return Transform_2d
is
inverse_Rotation : constant Matrix_2x2 := Transpose (Transform.Rotation);
begin
return (Translation => inverse_Rotation * (-Transform.Translation),
Rotation => inverse_Rotation);
end Invert;
function inverse_Transform (Transform : in Transform_2d; Vector : in Vector_2) return Vector_2
is
V : constant Vector_2 := Vector - Transform.Translation;
begin
return Transpose (Transform.Rotation) * V;
end inverse_Transform;
function get_Rotation (Transform : in Matrix_3x3) return Matrix_2x2
is
begin
return [[Transform (1, 1), Transform (1, 2)],
[Transform (2, 1), Transform (2, 2)]];
end get_Rotation;
procedure set_Rotation (Transform : in out Matrix_3x3; To : in Matrix_2x2)
is
begin
Transform (1, 1) := To (1, 1);
Transform (1, 2) := To (1, 2);
Transform (1, 3) := 0.0;
Transform (2, 1) := To (2, 1);
Transform (2, 2) := To (2, 2);
Transform (2, 3) := 0.0;
end set_Rotation;
function get_Translation (Transform : in Matrix_3x3) return Vector_2
is
begin
return [Transform (3, 1),
Transform (3, 2)];
end get_Translation;
procedure set_Translation (Transform : in out Matrix_3x3; To : in Vector_2)
is
begin
Transform (3, 1) := To (1);
Transform (3, 2) := To (2);
Transform (3, 3) := 1.0;
end set_Translation;
end any_Math.any_Algebra.any_linear.any_d2;

55
1-base/math/source/generic/pure/algebra/any_math-any_algebra-any_linear-any_d2.ads Normal file
View File

@@ -0,0 +1,55 @@
generic
package any_Math.any_Algebra.any_linear.any_d2
is
pragma Pure;
-----------
-- Vector_2
--
function Interpolated (From, To : in Vector_2; Percent : in unit_Percentage) return Vector_2;
function Distance (From, To : in Vector_2) return Real;
function Midpoint (From, To : in Vector_2) return Vector_2;
function Angle_between_pre_Norm (U, V : in Vector_2) return Radians;
--
-- Given that the vectors 'U' and 'V' are already normalized, returns a positive angle between 0 and 180 degrees.
-------------
-- Matrix_2x2
--
function to_Matrix (Row_1,
Row_2 : in Vector_2) return Matrix_2x2;
function to_rotation_Matrix (Angle : in Radians) return Matrix_2x2;
function up_Direction (Self : in Matrix_2x2) return Vector_2;
function right_Direction (Self : in Matrix_2x2) return Vector_2;
------------
-- Transform
--
function to_Transform (Rotation : in Matrix_2x2;
Translation : in Vector_2) return Matrix_3x3;
function to_Transform (From : in Transform_2d) return Matrix_3x3;
function to_translation_Transform (Translation : in Vector_2) return Matrix_3x3;
function to_rotation_Transform (Rotation : in Matrix_2x2) return Matrix_3x3;
function to_rotation_Transform (Angle : in Radians ) return Matrix_3x3;
function to_scale_Transform (Scale : in Vector_2) return Matrix_3x3;
function to_Transform_2d (From : in Matrix_3x3) return Transform_2d;
function to_Transform_2d (Rotation : in Radians;
Translation : in Vector_2) return Transform_2d;
function "*" (Left : in Vector_2; Right : in Transform_2d) return Vector_2;
function "*" (Left : in Vector_2; Right : in Matrix_3x3) return Vector_2;
function Invert (Transform : in Transform_2d) return Transform_2d;
function inverse_Transform (Transform : in Transform_2d; Vector : in Vector_2) return Vector_2;
function get_Rotation (Transform : in Matrix_3x3) return Matrix_2x2;
procedure set_Rotation (Transform : in out Matrix_3x3; To : in Matrix_2x2);
function get_Translation (Transform : in Matrix_3x3) return Vector_2;
procedure set_Translation (Transform : in out Matrix_3x3; To : in Vector_2);
end any_Math.any_Algebra.any_linear.any_d2;

1027
1-base/math/source/generic/pure/algebra/any_math-any_algebra-any_linear-any_d3.adb Normal file
View File

File diff suppressed because it is too large Load Diff

226
1-base/math/source/generic/pure/algebra/any_math-any_algebra-any_linear-any_d3.ads Normal file
View File

@@ -0,0 +1,226 @@
generic
package any_math.any_Algebra.any_linear.any_d3
is
pragma Pure;
------------
--- Vector_3
--
function Distance (From, To : in Vector_3) return Real;
function Midpoint (From, To : in Vector_3) return Vector_3;
function Interpolated (From, To : in Vector_3;
Percent : in unit_Percentage) return Vector_3;
function Angle_between_pre_Norm (U, V : in Vector_3) return Radians;
--
-- Given that the vectors 'U' and 'V' are already normalized, returns a positive angle between 0 and 180 degrees.
function Angle (Point_1,
Point_2,
Point_3 : in Vector_3) return Radians;
--
-- Returns the angle between the vector Point_1 to Point_2 and the vector Point_3 to Point_2.
--------------
--- Matrix_3x3
--
z_Up_to_y_Up : constant Matrix_3x3; -- Provides a rotation which may be multiplied
y_Up_to_z_Up : constant Matrix_3x3; -- by a vector to change co-ordinate systems.
function to_Matrix (Row_1,
Row_2,
Row_3 : in Vector_3) return Matrix_3x3;
function forward_Direction (Matrix : in Matrix_3x3) return Vector_3;
function up_Direction (Matrix : in Matrix_3x3) return Vector_3;
function right_Direction (Matrix : in Matrix_3x3) return Vector_3;
procedure re_Orthonormalise (Matrix : in out Matrix_3x3);
-------------
--- Rotations
--
function x_Rotation_from (Angle : in Radians) return Matrix_3x3;
function y_Rotation_from (Angle : in Radians) return Matrix_3x3;
function z_Rotation_from (Angle : in Radians) return Matrix_3x3;
function xyz_Rotation (x_Angle,
y_Angle,
z_Angle : in Real) return Matrix_3x3;
function xyz_Rotation (Angles : in Vector_3) return Matrix_3x3;
function to_Rotation (Axis : in Vector_3;
Angle : in Real) return Matrix_3x3;
function to_Rotation (Axis_x,
Axis_y,
Axis_z : in Real;
Rotation : in Radians) return Matrix_3x3;
--
-- Returns a rotation matrix describing a given rotation about an axis.
-- (TODO: Make this obsolescent and use the vector Axis version instead.)
---------
--- Euler
--
type Euler is new Vector_3;
--
-- 1: Roll
-- 2: Pitch
-- 3: Yaw
function to_Rotation (Angles : in Euler) return Matrix_3x3;
--
-- The euler angles are used to produce a rotation matrix. The euler
-- angles are applied in ZYX order. That is, a vector is first rotated
-- about X, then Y and then Z.
-----------
--- General
--
function Look_at (Eye,
Center,
Up : in Vector_3) return Matrix_4x4;
function to_Viewport_Transform (Origin,
Extent : in Vector_2) return Matrix_4x4;
function to_Perspective (FoVy : in Degrees;
Aspect,
zNear,
zFar : in Real) return Matrix_4x4;
-------------
--- Transform
--
function to_Translation_Matrix (Translation : in Vector_3) return Matrix_4x4;
function to_Transform (Matrix : in Matrix_4x4) return Transform_3d;
function "*" (Left : in Transform_3d; Right : in Vector_3) return Vector_3;
function "*" (Left : in Vector_3; Right : in Transform_3d) return Vector_3;
function "*" (Left : in Transform_3d; Right : in Transform_3d) return Transform_3d;
function "*" (Left : in Vector_3; Right : in Matrix_4x4) return Vector_3;
function Invert (Transform : in Transform_3d) return Transform_3d;
function inverse_Transform (Transform : in Transform_3d; Vector : in Vector_3) return Vector_3;
--------------
--- Quaternion
--
procedure set_from_Matrix_3x3 (Quat : out Quaternion; Matrix : in Matrix_3x3);
function to_Quaternion (Matrix : in Matrix_3x3) return Quaternion;
function to_Matrix (Quat : in Quaternion) return Matrix_3x3;
function Norm (Quat : in Quaternion) return Real;
function Versor (Quat : in Quaternion) return Quaternion; -- Produces the unit quaternion of Quat.
function Farthest (Quat : in Quaternion; qd : in Quaternion) return Quaternion; -- TODO: Document this.
function Invert (Quat : in Quaternion) return Quaternion;
function Angle (Quat : in Quaternion) return Radians;
function Axis (Quat : in Quaternion) return Vector_3;
function "*" (Left, Right : in Quaternion) return Real; -- Dot product.
function "*" (Left, Right : in Quaternion) return Quaternion; -- Cross product.
function "+" (Left, Right : in Quaternion) return Quaternion;
function "-" (Left, Right : in Quaternion) return Quaternion;
function "-" (Quat : in Quaternion) return Quaternion;
function "*" (Left : in Quaternion; Right : in Vector_3) return Quaternion;
function "*" (Left : in Vector_3; Right : in Quaternion) return Quaternion;
function "*" (Left : in Quaternion; Right : in Real) return Quaternion;
function Interpolated (From,
To : in Quaternion;
Percent : in unit_Percentage) return Quaternion;
--
-- Return the quaternion which is the result of spherical linear interpolation (Slerp) between Initial and Final.
-- Percent is the ratio between 'From' and 'To' to interpolate.
-- If Percent = 0.0 the result is Initial.
-- If Percent = 100.0 the result is Final.
-- Interpolates assuming constant velocity.
------------
--- Vector_4
--
function "/" (Left, Right : in Vector_4) return Vector_4;
function max_Axis (Vector : in Vector_4) return Integer;
function closest_Axis (Vector : in Vector_4) return Integer;
function to_transform_Matrix (Transform : in Transform_3d) return Matrix_4x4;
function to_transform_Matrix (Rotation : in Matrix_3x3;
Translation : in Vector_3) return Matrix_4x4;
function to_rotate_Matrix (Rotation : in Matrix_3x3) return Matrix_4x4;
function to_translate_Matrix (Translation : in Vector_3) return Matrix_4x4;
function to_scale_Matrix (Scale : in Vector_3) return Matrix_4x4;
----------------------
--- Transform Matrices
--
function get_Rotation (Transform : in Matrix_4x4) return Matrix_3x3;
procedure set_Rotation (Transform : in out Matrix_4x4; To : in Matrix_3x3);
function get_Translation (Transform : in Matrix_4x4) return Vector_3;
procedure set_Translation (Transform : in out Matrix_4x4; To : in Vector_3);
function inverse_Rotation (Rotation : in Matrix_3x3) return Matrix_3x3;
function inverse_Transform (Transform : in Matrix_4x4) return Matrix_4x4;
--------------
--- un-Project
--
type Rectangle is
record
Min : Integers (1 .. 2); -- Bottom left corner.
Max : Integers (1 .. 2); -- Upper right corner.
end record;
function unProject (From : in Vector_3;
Model : in Matrix_4x4;
Projection : in Matrix_4x4;
Viewport : in Rectangle) return Vector_3;
--
-- Maps the 'From' window space coordinates into object space coordinates using Model, Projection and Viewport.
private
z_Up_to_y_Up : constant Matrix_3x3 := [[1.0, 0.0, 0.0],
[0.0, 0.0, 1.0],
[0.0, -1.0, 0.0]];
y_Up_to_z_Up : constant Matrix_3x3 := [[1.0, 0.0, 0.0],
[0.0, 0.0, -1.0],
[0.0, 1.0, 0.0]];
pragma Inline ("+");
pragma Inline ("-");
pragma Inline ("*");
end any_math.any_Algebra.any_linear.any_d3;

589
1-base/math/source/generic/pure/algebra/any_math-any_algebra-any_linear.adb Normal file
View File

@@ -0,0 +1,589 @@
with
ada.Characters.latin_1;
package body any_Math.any_Algebra.any_linear
is
-----------
--- Vectors
--
function Norm_squared (Self : in Vector) return Real
is
Norm_2 : Real := 0.0;
begin
for Each in Self'Range
loop
Norm_2 := Norm_2 + Self (Each) * Self (Each);
end loop;
return Norm_2;
end Norm_squared;
procedure normalise (Self : in out Vector)
is
use Vectors;
inverse_Norm : constant Real := 1.0 / abs Self;
begin
for Each in Self'Range
loop
Self (Each) := Self (Each) * inverse_Norm;
end loop;
end normalise;
function Normalised (Self : in Vector) return Vector
is
Result : Vector := Self;
begin
normalise (Result);
return Result;
end Normalised;
procedure normalise (Self : in out Vector_2)
is
inverse_Norm : constant Real := 1.0 / abs Self;
begin
Self := Self * inverse_Norm;
end normalise;
function Normalised (Self : in Vector_2) return Vector_2
is
inverse_Norm : constant Real := 1.0 / abs Self;
begin
return Self * inverse_Norm;
end Normalised;
procedure normalise (Self : in out Vector_3)
is
inverse_Norm : constant Real := 1.0 / abs Self;
begin
Self := Self * inverse_Norm;
end normalise;
function Normalised (Self : in Vector_3) return Vector_3
is
inverse_Norm : constant Real := 1.0 / abs Self;
begin
return Self * inverse_Norm;
end Normalised;
function Min (Left, Right : in Vector) return Vector
is
Min : Vector (Left'Range);
begin
pragma Assert (Left'Length = Right'Length);
for Each in Min'Range
loop
Min (Each) := Real'Min (Left (Each),
Right (Each));
end loop;
return Min;
end Min;
function Max (Left, Right : in Vector) return Vector
is
Max : Vector (Left'Range);
begin
pragma Assert (Left'Length = Right'Length);
for Each in Max'Range
loop
Max (Each) := Real'Max (Left (Each),
Right (Each));
end loop;
return Max;
end Max;
function scaled (Self : in Vector; By : in Vector) return Vector
is
Result : Vector (Self'Range);
begin
for Each in Result'Range
loop
Result (Each) := Self (Each) * By (Each);
end loop;
return Result;
end scaled;
------------
--- Matrices
--
function to_Matrix (Row_1, Row_2, Row_3 : in Vector_3) return Matrix_3x3
is
begin
return [[Row_1 (1), Row_1 (2), Row_1 (3)],
[Row_2 (1), Row_2 (2), Row_2 (3)],
[Row_3 (1), Row_3 (2), Row_3 (3)]];
end to_Matrix;
function Min (Self : in Matrix) return Real
is
Min : Real := Real'Last;
begin
for each_Row in Self'Range (1)
loop
for each_Col in Self'Range (2)
loop
Min := Real'Min (Min,
Self (each_Row, each_Col));
end loop;
end loop;
return Min;
end Min;
function Max (Self : in Matrix) return Real
is
Max : Real := Real'First;
begin
for each_Row in Self'Range (1)
loop
for each_Col in Self'Range (2)
loop
Max := Real'Max (Max,
Self (each_Row, each_Col));
end loop;
end loop;
return Max;
end Max;
function Image (Self : in Matrix) return String
is
Image : String (1 .. 1024 * 1024); -- Handles one megabyte image, excess is truncated.
Count : Standard.Natural := 0;
procedure add (Text : in String)
is
begin
Image (Count + 1 .. Count + text'Length) := Text;
Count := Count + text'Length;
end add;
begin
add ("(");
for Row in self'Range (1)
loop
add ([1 => ada.Characters.latin_1.LF]);
if Row /= self'First (1)
then
add (", ");
end if;
for Col in self'Range (2)
loop
if Col /= self'First (2)
then
add (", ");
end if;
add (Real'Image (Self (Row, Col)));
end loop;
end loop;
add (")");
return Image (1 .. Count);
exception
when others =>
return Image (1 .. Count);
end Image;
function is_Square (Self : in Matrix) return Boolean
is
begin
return Self'Length (1) = Self'Length (2);
end is_Square;
function sub_Matrix (Self : in Matrix; start_Row, end_Row : in Index;
start_Col, end_Col : in Index) return Matrix
is
sub_Matrix : Matrix (1 .. end_Row - start_Row + 1,
1 .. end_Col - start_Col + 1);
begin
for each_Row in sub_Matrix'Range (1)
loop
for each_Col in sub_Matrix'Range (2)
loop
sub_Matrix (each_Row, each_Col) := Self (each_Row + start_Row - 1,
each_Col + start_Col - 1);
end loop;
end loop;
return sub_Matrix;
end sub_Matrix;
function Identity (Size : in Index := 3) return Matrix
is
Result : Matrix (1 .. Size, 1 .. Size);
begin
for Row in 1 .. Size
loop
for Col in 1 .. Size
loop
if Row = Col
then Result (Row, Col) := 1.0;
else Result (Row, Col) := 0.0;
end if;
end loop;
end loop;
return Result;
end Identity;
procedure invert (Self : in out Matrix)
is
use Vectors;
begin
Self := Inverse (Self);
end invert;
---------------
--- Quaternions
--
function to_Quaternion (axis_X,
axis_Y,
axis_Z : in Real;
Angle : in Real) return Quaternion
is
Result : Quaternion;
L : Real := axis_X * axis_X + axis_Y * axis_Y + axis_Z * axis_Z;
begin
if L > 0.0
then
declare
use Functions;
half_Angle : constant Real := Angle * 0.5;
begin
Result.R := Cos (half_Angle);
L := Sin (half_Angle) * (1.0 / SqRt (L));
Result.V (1) := axis_X * L;
Result.V (2) := axis_Y * L;
Result.V (3) := axis_Z * L;
end;
else
Result.R := L;
Result.V (1) := 0.0;
Result.V (2) := 0.0;
Result.V (3) := 0.0;
end if;
return Result;
end to_Quaternion;
function to_Quaternion (Axis : in Vector_3;
Angle : in Real) return Quaternion
is
Result : Quaternion;
L : Real := Axis * Axis;
begin
if L > 0.0
then
declare
use Functions;
half_Angle : constant Real := Angle * 0.5;
begin
Result.R := Cos (half_Angle);
L := Sin (half_Angle) * (1.0 / SqRt (L));
Result.V := Axis * L;
end;
else
Result.R := L;
Result.V := [0.0, 0.0, 0.0];
end if;
return Result;
end to_Quaternion;
function "*" (Self : in Quaternion;
By : in Quaternion) return Quaternion
is
x : constant := 1;
y : constant := 2;
z : constant := 3;
A : Quaternion renames Self;
B : Quaternion renames By;
AtBt : constant Real := A.R * B.R;
AxBx : constant Real := A.V (x) * B.V (x);
AyBy : constant Real := A.V (y) * B.V (y);
AzBz : constant Real := A.V (z) * B.V (z);
AtBx : constant Real := A.R * B.V (x);
AxBt : constant Real := A.V (x) * B.R;
AyBz : constant Real := A.V (y) * B.V (z);
AzBy : constant Real := A.V (z) * B.V (y);
AtBy : constant Real := A.R * B.V (y);
AxBz : constant Real := A.V (x) * B.V (z);
AyBt : constant Real := A.V (y) * B.R;
AzBx : constant Real := A.V (z) * B.V (x);
AtBz : constant Real := A.R * B.V (z);
AxBy : constant Real := A.V (x) * B.V (y);
AyBx : constant Real := A.V (y) * B.V (x);
AzBt : constant Real := A.V (z) * B.R;
begin
return (R => AtBt - AxBx - AyBy - AzBz,
V => [AtBx + AxBt + AyBz - AzBy,
AtBy - AxBz + AyBt + AzBx,
AtBz + AxBy - AyBx + AzBt]);
end "*";
function Unit (Self : in Quaternion) return Quaternion
is
begin
return to_Quaternion ( to_Vector (Self)
/ abs to_Vector (Self));
end Unit;
function infinitesimal_Rotation_from (Self : in Quaternion;
angular_Velocity : in Vector_3) return Quaternion
is
i_Rotation : Quaternion;
begin
i_Rotation.R := 0.5 * (- angular_Velocity (1) * Self.V (1)
- angular_Velocity (2) * Self.V (2)
- angular_Velocity (3) * Self.V (3));
i_Rotation.V (1) := 0.5 * ( angular_Velocity (1) * Self.R
+ angular_Velocity (2) * Self.V (3)
- angular_Velocity (3) * Self.V (2));
i_Rotation.V (2) := 0.5 * (- angular_Velocity (1) * Self.V (3)
+ angular_Velocity (2) * Self.R
+ angular_Velocity (3) * Self.V (1));
i_Rotation.V (3) := 0.5 * ( angular_Velocity (1) * Self.V (2)
- angular_Velocity (2) * Self.V (1)
+ angular_Velocity (3) * Self.R);
return i_Rotation;
end infinitesimal_Rotation_from;
function euler_Angles (Self : in Quaternion) return Vector_3 -- 'Self' can be a non-normalised quaternion.
is
use Functions;
w : Real renames Self.R;
x : Real renames Self.V (1);
y : Real renames Self.V (2);
z : Real renames Self.V (3);
the_Angles : Vector_3;
Bank : Real renames the_Angles (1);
Heading : Real renames the_Angles (2);
Attitude : Real renames the_Angles (3);
sqw : constant Real := w * w;
sqx : constant Real := x * x;
sqy : constant Real := y * y;
sqz : constant Real := z * z;
unit : constant Real := sqx + sqy + sqz + sqw; -- If normalised then is 1.0 else is a correction factor.
test : constant Real := x * y + z * w;
begin
if test > 0.499 * unit
then -- Singularity at north pole.
Heading := 2.0 * arcTan (x, w);
Attitude := Pi / 2.0;
Bank := 0.0;
return the_Angles;
end if;
if test < -0.499 * unit
then -- Singularity at south pole.
Heading := -2.0 * arcTan (x, w);
Attitude := -Pi / 2.0;
Bank := 0.0;
return the_Angles;
end if;
Heading := arcTan (2.0 * y * w - 2.0 * x * z, sqx - sqy - sqz + sqw);
Bank := arcTan (2.0 * x * w - 2.0 * y * z, -sqx + sqy - sqz + sqw);
Attitude := arcSin (2.0 * test / unit);
return the_Angles;
end euler_Angles;
function to_Quaternion (Self : in Matrix_3x3) return Quaternion
is
use Functions;
TR : Real;
S : Real;
Result : Quaternion;
begin
TR := Self (1, 1) + Self (2, 2) + Self (3, 3);
if TR >= 0.0
then
S := SqRt (TR + 1.0);
Result.R := 0.5 * S;
S := 0.5 * (1.0 / S);
Result.V (1) := (Self (3, 2) - Self (2, 3)) * S;
Result.V (2) := (Self (1, 3) - Self (3, 1)) * S;
Result.V (3) := (Self (2, 1) - Self (1, 2)) * S;
return Result;
end if;
-- Otherwise, find the largest diagonal element and apply the appropriate case.
--
declare
function case_1_Result return Quaternion
is
begin
S := SqRt (Self (1, 1) - (Self (2, 2) + Self (3, 3)) + 1.0);
Result.V (1) := 0.5 * S;
S := 0.5 * (1.0 / S);
Result.V (2) := (Self (1, 2) + Self (2, 1)) * S;
Result.V (3) := (Self (3, 1) + Self (1, 3)) * S;
Result.R := (Self (3, 2) - Self (2, 3)) * S;
return Result;
end case_1_Result;
function case_2_Result return Quaternion
is
begin
S := SqRt (Self (2, 2) - (Self (3, 3) + Self (1, 1)) + 1.0);
Result.V (2) := 0.5 * S;
S := 0.5 * (1.0 / S);
Result.V (3) := (Self (2, 3) + Self (3, 2)) * S;
Result.V (1) := (Self (1, 2) + Self (2, 1)) * S;
Result.R := (Self (1, 3) - Self (3, 1)) * S;
return Result;
end case_2_Result;
function case_3_Result return Quaternion
is
begin
S := SqRt (Self (3, 3) - (Self (1, 1) + Self (2, 2)) + 1.0);
Result.V (3) := 0.5 * S;
S := 0.5 * (1.0 / S);
Result.V (1) := (Self (3, 1) + Self (1, 3)) * S;
Result.V (2) := (Self (2, 3) + Self (3, 2)) * S;
Result.R := (Self (2, 1) - Self (1, 2)) * S;
return Result;
end case_3_Result;
pragma Inline (case_1_Result);
pragma Inline (case_2_Result);
pragma Inline (case_3_Result);
begin
if Self (2, 2) > Self (1, 1)
then
if Self (3, 3) > Self (2, 2)
then
return case_3_Result;
end if;
return case_2_Result;
end if;
if Self (3, 3) > Self (1, 1)
then
return case_3_Result;
end if;
return case_1_Result;
end;
end to_Quaternion;
function Conjugate (Self : in Quaternion) return Quaternion
is
begin
return (Self.R, -Self.V);
end conjugate;
procedure normalise (Self : in out Quaternion)
is
begin
Self := Normalised (Self);
end normalise;
function Normalised (Self : in Quaternion) return Quaternion
is
begin
return to_Quaternion (Vector_4 (Normalised (Vector (to_Vector (Self)))));
end Normalised;
end any_Math.any_Algebra.any_linear;

103
1-base/math/source/generic/pure/algebra/any_math-any_algebra-any_linear.ads Normal file
View File

@@ -0,0 +1,103 @@
generic
package any_Math.any_Algebra.any_linear
is
pragma Pure;
----------
-- Vector
--
function Norm_squared (Self : in Vector) return Real; -- Length squared.
function Normalised (Self : in Vector) return Vector;
procedure Normalise (Self : in out Vector);
function Normalised (Self : in Vector_2) return Vector_2;
procedure Normalise (Self : in out Vector_2);
function Normalised (Self : in Vector_3) return Vector_3;
procedure Normalise (Self : in out Vector_3);
function Min (Left, Right : in Vector) return Vector;
function Max (Left, Right : in Vector) return Vector;
function Scaled (Self : in Vector; By : in Vector) return Vector;
----------
-- Matrix
--
function to_Matrix (Row_1,
Row_2,
Row_3 : in Vector_3) return Matrix_3x3;
function Identity (Size : in Index := 3) return Matrix;
function Min (Self : in Matrix) return Real;
function Max (Self : in Matrix) return Real;
function Image (Self : in Matrix) return String;
procedure invert (Self : in out Matrix);
function is_Square (Self : in Matrix) return Boolean;
function sub_Matrix (Self : in Matrix;
start_Row, end_Row : in Index;
start_Col, end_Col : in Index) return Matrix;
---------------
-- Quaternion
--
function to_Quaternion (axis_X,
axis_Y,
axis_Z : in Real;
Angle : in Real) return Quaternion;
--
-- Returns a quaternion defined by a rotation about an axis.
-- (TODO: rid this and use Vector_3 version instead.)
function to_Quaternion (Axis : in Vector_3;
Angle : in Real) return Quaternion;
--
-- Returns a quaternion defined by a rotation about an axis.
function to_Quaternion (Self : in Matrix_3x3) return Quaternion;
function "*" (Self : in Quaternion; By : in Quaternion) return Quaternion;
--
-- Grassmann product.
function Unit (Self : in Quaternion) return Quaternion;
function Conjugate (Self : in Quaternion) return Quaternion;
--
-- (TODO: only for unit quaternions.)
function euler_Angles (Self : in Quaternion) return Vector_3;
function infinitesimal_Rotation_from
(Self : in Quaternion; angular_Velocity : in Vector_3) return Quaternion;
--
-- An infinitesimal rotation may be multiplied by a duration and then added to the original attitude
-- to produce the attitude at the given time.
function Normalised (Self : in Quaternion) return Quaternion;
procedure normalise (Self : in out Quaternion);
private
pragma Inline ("*");
pragma Inline_Always (Norm_squared);
pragma Inline_Always (Normalise);
end any_Math.any_Algebra.any_linear;

7
1-base/math/source/generic/pure/algebra/any_math-any_algebra.ads Normal file
View File

@@ -0,0 +1,7 @@
generic
package any_Math.any_Algebra
is
pragma Pure;
end any_Math.any_Algebra;

8
1-base/math/source/generic/pure/analysis/any_math-any_analysis.ads Normal file
View File

@@ -0,0 +1,8 @@
generic
package any_Math.any_Analysis
is
pragma Pure;
pragma Optimize (Time);
end any_Math.any_Analysis;

8
1-base/math/source/generic/pure/arithmetic/any_math-any_arithmetic.ads Normal file
View File

@@ -0,0 +1,8 @@
generic
package any_Math.any_Arithmetic
is
pragma Pure;
pragma Optimize (Time);
end any_Math.any_Arithmetic;

192
1-base/math/source/generic/pure/geometry/any_math-any_geometry-any_d2-any_hexagon.adb Normal file
View File

@@ -0,0 +1,192 @@
package body any_Math.any_Geometry.any_d2.any_Hexagon
is
-------------
--- vertex_Id
--
function prior_Vertex (to_Vertex : in vertex_Id) return vertex_Id
is
begin
if To_Vertex = 1
then return 6;
else return to_Vertex - 1;
end if;
end prior_Vertex;
function next_Vertex (to_Vertex : in vertex_Id) return vertex_Id
is
begin
if to_Vertex = 6
then return 1;
else return to_Vertex + 1;
end if;
end next_Vertex;
-----------
--- Hexagon
--
function to_Hexagon (circumRadius : in Real) return Item
is
begin
return (circumRadius => circumRadius);
end to_Hexagon;
function circumRadius (Self : in Item) return Real
is
begin
return Self.circumRadius;
end circumRadius;
function Area (Self : in Item) return Real
is
begin
return 3.0 * R (Self) * inRadius (Self);
end Area;
function Perimeter (Self : in Item) return Real
is
begin
return 6.0 * side_Length (Self);
end Perimeter;
function Angle (Self : in Item with Unreferenced;
at_Vertex : in vertex_Id with Unreferenced) return Radians
is
begin
return to_Radians (120.0);
end Angle;
function minimal_Diameter (Self : in Item) return Real
is
begin
return 2.0 * inRadius (Self);
end minimal_Diameter;
function maximal_Diameter (Self : in Item) return Real
is
begin
return 2.0 * Self.circumRadius;
end maximal_Diameter;
function inRadius (Self : in Item) return Real
is
use Functions;
begin
return cos (to_Radians (30.0)) * R (Self);
end inRadius;
function side_Length (Self : in Item) return Real
is
begin
return Self.circumRadius;
end side_Length;
function Site (Self : in Item; of_Vertex : in vertex_Id) return any_d2.Site
is
use Functions;
Angle : constant Radians := to_Radians (60.0 * Degrees (of_Vertex - 1));
begin
return any_d2.Site' (1 => Self.circumRadius * cos (Angle),
2 => Self.circumRadius * sin (Angle));
end Site;
function horizontal_Distance (Self : in Item) return Real
is
begin
return Width (Self) * 3.0 / 4.0;
end horizontal_Distance;
function vertical_Distance (Self : in Item) return Real
is
begin
return Height (Self);
end vertical_Distance;
--------
--- Grid
--
function to_Grid (Rows, Cols : in Positive;
circumRadius : in Real) return Grid
is
Hex : constant Item := (circumRadius => circumRadius);
inRadius : constant Real := any_Hexagon.inRadius (Hex);
maximal_Diameter : constant Real := any_Hexagon.maximal_Diameter (Hex);
minimal_Diameter : constant Real := any_Hexagon.minimal_Diameter (Hex);
Result : Grid (Rows, Cols);
begin
Result.circumRadius := circumRadius;
for Row in 1 .. Rows
loop
for Col in 1 .. Cols
loop
Result.Centers (Row, Col) := [circumRadius + Real (Col - 1) * (maximal_Diameter - 0.5 * circumRadius),
inRadius + Real (Row - 1) * minimal_Diameter];
if Col mod 2 = 0 -- Even column.
then
Result.Centers (Row, Col) (2) := @ + inRadius;
end if;
end loop;
end loop;
return Result;
end to_Grid;
function hex_Center (Grid : in any_Hexagon.Grid; Coords : in Coordinates) return any_d2.Site
is
begin
return Grid.Centers (Coords.Row, Coords.Col);
end hex_Center;
function vertex_Site (Self : in Grid; hex_Id : in any_Hexagon.Coordinates;
Which : in any_Hexagon.vertex_Id) return any_d2.Site
is
Hex : constant Item := (circumRadius => Self.circumRadius);
begin
return hex_Center (Self, hex_Id)
+ Site (Hex, of_Vertex => Which);
end vertex_Site;
end any_Math.any_Geometry.any_d2.any_Hexagon;

109
1-base/math/source/generic/pure/geometry/any_math-any_geometry-any_d2-any_hexagon.ads Normal file
View File

@@ -0,0 +1,109 @@
generic
package any_Math.any_Geometry.any_d2.any_Hexagon with Pure
--
-- Models a regular, flat-topped hexagon.
--
-- https://en.wikipedia.org/wiki/Hexagon
--
--
-- 5 6
-- ---
-- 4/ \1
-- \ /
-- ---
-- 3 2
--
is
-------------
--- vertex_Id
--
subtype vertex_Id is any_Geometry.vertex_Id range 1 .. 6;
function prior_Vertex (to_Vertex : in vertex_Id) return vertex_Id;
function next_Vertex (to_Vertex : in vertex_Id) return vertex_Id;
--------
--- Item
--
type Item is private;
function to_Hexagon (circumRadius : in Real) return Item;
function maximal_Diameter (Self : in Item) return Real;
function minimal_Diameter (Self : in Item) return Real; -- 'd'
function circumRadius (Self : in Item) return Real; -- 'r'
function inRadius (Self : in Item) return Real; -- 'r'
function Area (Self : in Item) return Real;
function Perimeter (Self : in Item) return Real;
function Width (Self : in Item) return Real renames maximal_Diameter;
function Height (Self : in Item) return Real renames minimal_Diameter;
function side_Length (Self : in Item) return Real;
function Site (Self : in Item; of_Vertex : in vertex_Id) return any_d2.Site;
function Angle (Self : in Item; at_Vertex : in vertex_Id) return Radians;
function R (Self : in Item) return Real renames circumRadius;
function D (Self : in Item) return Real renames maximal_Diameter;
function t (Self : in Item) return Real renames side_Length;
function horizontal_Distance (Self : in Item) return Real; -- The distance between adjacent
function vertical_Distance (Self : in Item) return Real; -- hexagon centers.
--------
--- Grid
--
-- Origin is at the top left corner.
-- X increases to the right.
-- Y increases downwards.
--
type Grid (Rows : Positive;
Cols : Positive) is private;
type Coordinates is
record
Row, Col: Positive;
end record;
function to_Grid (Rows, Cols : in Positive;
circumRadius : in Real) return Grid;
function hex_Center (Grid : in any_Hexagon.Grid; Coords : in Coordinates) return any_d2.Site;
--
-- Returns the centre of the hexagon at the given co-ordinates.
function vertex_Site (Self : in Grid; hex_Id : in any_Hexagon.Coordinates;
Which : in any_Hexagon.vertex_Id) return any_d2.Site;
private
type Item is
record
circumRadius : Real;
end record;
type Grid (Rows : Positive;
Cols : Positive) is
record
circumRadius : Real;
Centers : any_d2.Grid (1 .. Rows,
1 .. Cols);
end record;
end any_Math.any_Geometry.any_d2.any_Hexagon;

584
1-base/math/source/generic/pure/geometry/any_math-any_geometry-any_d2.adb Normal file
View File

@@ -0,0 +1,584 @@
package body any_Math.any_Geometry.any_d2
is
---------
-- Sites
--
function Distance (From, To : Site) return Real
is
use Functions;
begin
return SqRt ( (To (1) - From (1)) ** 2
+ (To (2) - From (2)) ** 2);
end Distance;
function to_Polar (Self : in Site) return polar_Site
is
use any_Math.complex_Reals;
the_Complex : constant Complex := compose_from_Cartesian (Self (1),
Self (2));
begin
return (Angle => Argument (the_Complex),
Extent => Modulus (the_Complex));
end to_Polar;
function to_Site (Self : in polar_Site) return Site
is
use any_Math.complex_Reals;
the_Complex : constant Complex := compose_from_Polar (Modulus => Self.Extent,
Argument => Self.Angle);
begin
return [the_Complex.Re,
the_Complex.Im];
end to_Site;
function Angle (Self : in Site) return Radians
is
use any_Math.complex_Reals;
the_Complex : constant Complex := compose_from_Cartesian (Self (1),
Self (2));
begin
return Argument (the_Complex);
end Angle;
function Extent (Self : in Site) return Real
is
use any_Math.complex_Reals;
the_Complex : constant Complex := compose_from_Cartesian (Self (1),
Self (2));
begin
return Modulus (the_Complex);
end Extent;
---------
-- Lines
--
function to_Line (Anchor : in Site;
Angle : in Radians) return Line
is
use Functions;
begin
return (Kind => anchored_Gradient,
Anchor => Anchor,
Gradient => Tan (Angle)); -- TODO: What about infinite gradient ? ie 90 and 270 degrees ?
end to_Line;
function to_Line (Site_1,
Site_2 : in Site) return Line
is
begin
return (Kind => two_Points,
Sites => [Site_1,
Site_2]);
end to_Line;
function X_of (Self : in Line; Y : in Real) return Real
is
begin
return
(Y - Self.Anchor (2)) / Self.Gradient
+ Self.Anchor (1);
end X_of;
function Y_of (Self : in Line; X : in Real) return Real
is
begin
return
Self.Gradient * (X - Self.Anchor (1))
+ Self.Anchor (2);
end Y_of;
function Gradient (Self : in Line) return Real
is
Run : constant Real := Self.Sites (2)(1) - Self.Sites (1)(1);
begin
if Run = 0.0
then
return Real'Last;
else
return (Self.Sites (2) (2) - Self.Sites (1) (2))
/ Run;
end if;
end Gradient;
----------
-- Bounds
--
function to_bounding_Box (Self : Sites) return bounding_Box
is
Result : bounding_Box := null_Bounds;
begin
for Each in Self'Range
loop
Result.Lower (1) := Real'Min (Result.Lower (1), Self (Each)(1));
Result.Lower (2) := Real'Min (Result.Lower (2), Self (Each)(2));
Result.Upper (1) := Real'Max (Result.Upper (1), Self (Each)(1));
Result.Upper (2) := Real'Max (Result.Upper (2), Self (Each)(2));
end loop;
return Result;
end to_bounding_Box;
function Extent (Self : in bounding_Box; Dimension : in Index) return Real
is
begin
return Self.Upper (Dimension) - Self.Lower (Dimension);
end Extent;
function "or" (Left : in bounding_Box; Right : in Site) return bounding_Box
is
Result : bounding_Box;
begin
for i in Right'Range
loop
if Right (i) < Left.Lower (i)
then Result.Lower (i) := Right (i);
else Result.Lower (i) := Left.Lower (i);
end if;
if Right (i) > Left.Upper (i)
then Result.Upper (i) := Right (i);
else Result.Upper (i) := Left.Upper (i);
end if;
end loop;
return Result;
end "or";
function "or" (Left : in bounding_Box; Right : in bounding_Box) return bounding_Box
is
Result : bounding_Box := Left or Right.Lower;
begin
Result := Result or Right.Upper;
return Result;
end "or";
function "+" (Left : in bounding_Box; Right : in Vector_2) return bounding_Box
is
begin
return (Left.Lower + Right,
Left.Upper + Right);
end "+";
function Image (Self : bounding_Box) return String
is
begin
return "(Lower => " & Image (Self.Lower)
& ", Upper => " & Image (Self.Upper) & ")";
end Image;
------------
-- Triangles
--
procedure check (Self : in Triangle)
is
begin
if Self.Vertices (1) = Self.Vertices (2)
or Self.Vertices (1) = Self.Vertices (3)
or Self.Vertices (2) = Self.Vertices (3)
then
raise Degenerate;
end if;
declare
L1 : constant Line := to_Line (Self.Vertices (1), Self.Vertices (2));
L2 : constant Line := to_Line (Self.Vertices (2), Self.Vertices (3));
L3 : constant Line := to_Line (Self.Vertices (3), Self.Vertices (1));
M1 : constant Real := Gradient (L1);
M2 : constant Real := Gradient (L2);
M3 : constant Real := Gradient (L3);
begin
if M1 = M2
or M1 = M3
or M2 = M3
then
raise Colinear with
" G1: " & Image (M1)
& " G2: " & Image (M2)
& " G3: " & Image (M3);
end if;
end;
end check;
pragma Unreferenced (check);
-- function Area (Self : in Triangle) return Real
-- --
-- -- This is an implementation of Heron's formula.
-- -- It is numerically unstable with very small angles.
-- --
-- is
-- use Functions;
--
-- A : constant Real := Distance (Self.Vertices (1), Self.Vertices (2));
-- B : constant Real := Distance (Self.Vertices (2), Self.Vertices (3));
-- C : constant Real := Distance (Self.Vertices (3), Self.Vertices (1));
--
-- S : constant Real := (A + B + C) / 2.0; -- Semi-perimeter.
--
-- begin
-- return Real (SqRt (S * (S - A) * (S - B) * (S - C))); -- Herons formula.
-- end Area;
function Area (Self : in Triangle) return Real
--
-- This is a numerically stable implementation of Heron's formula.
-- See: https://en.wikipedia.org/wiki/Heron%27s_formula#Numerical_stability.
--
is
use Functions;
a : Real := Distance (Self.Vertices (1), Self.Vertices (2));
b : Real := Distance (Self.Vertices (2), Self.Vertices (3));
c : Real := Distance (Self.Vertices (3), Self.Vertices (1));
D : Real;
begin
-- Sort the lengths such that a >= b >= c.
--
if c > b then swap (b, c); end if;
if a < b then swap (a, b); end if;
if b < c then swap (b, c); end if;
D := (a + (b + c))
* (c - (a - b))
* (c + (a - b))
* (a + (b - c));
if D <= 0.0
then
return 0.0;
end if;
return 0.25 * SqRt (D);
end Area;
function Perimeter (Self : Triangle) return Real
is
begin
return
Distance (Self.Vertices (1), Self.Vertices (2))
+ Distance (Self.Vertices (2), Self.Vertices (3))
+ Distance (Self.Vertices (3), Self.Vertices (1));
end Perimeter;
function prior_Vertex (Self : in Triangle; to_Vertex : in Positive) return Site
is
begin
if to_Vertex = 1
then return Self.Vertices (3);
else return Self.Vertices (to_Vertex - 1);
end if;
end prior_Vertex;
function next_Vertex (Self : in Triangle; to_Vertex : in Positive) return Site
is
begin
if to_Vertex = 3
then return Self.Vertices (1);
else return Self.Vertices (to_Vertex + 1);
end if;
end next_Vertex;
function Angle (Self : in Triangle; at_Vertex : in Positive) return Radians
is
use Functions;
a : constant Real := Distance (next_Vertex (Self, to_vertex => at_Vertex),
prior_Vertex (Self, to_vertex => at_Vertex));
b : constant Real := Distance (Self.Vertices (at_Vertex), next_Vertex (Self, to_vertex => at_Vertex));
c : constant Real := Distance (Self.Vertices (at_Vertex), prior_Vertex (Self, to_vertex => at_Vertex));
cos_A : constant Real := (b**2 + c**2 - a**2) / (2.0 * b * c);
begin
if cos_A < -1.0 then return to_Radians (180.0);
elsif cos_A > 1.0 then return 0.0;
else return arcCos (cos_A);
end if;
end Angle;
----------
-- Circles
--
function Area (Self : Circle) return Real
is
begin
return Pi * Self.Radius**2;
end Area;
function Perimeter (Self : Circle) return Real
is
begin
return 2.0 * Pi * Self.Radius;
end Perimeter;
-----------
-- Polygons
--
function Centroid (Self : in Polygon) return Site
is
Result : Site := Origin_2d;
begin
for i in 1 .. Self.Vertex_Count
loop
Result := Result + Self.Vertices (i);
end loop;
Result := Result / Real (Self.Vertex_Count);
return Result;
end Centroid;
procedure center (Self : in out Polygon)
is
Center : constant Site := Centroid (Self);
begin
for i in 1 .. Self.Vertex_Count
loop
Self.Vertices (i) := Self.Vertices (i) - Center;
end loop;
end center;
function prior_Vertex (Self : in Polygon; to_Vertex : in Positive) return Site
is
begin
if To_Vertex = 1
then return Self.Vertices (Self.Vertex_Count);
else return Self.Vertices (to_Vertex - 1);
end if;
end prior_Vertex;
function next_Vertex (Self : in Polygon; to_Vertex : in Positive) return Site
is
begin
if to_Vertex = Self.Vertex_Count
then return Self.Vertices (1);
else return Self.Vertices (to_Vertex + 1);
end if;
end next_Vertex;
function is_Triangle (Self : in Polygon) return Boolean
is
begin
return Self.Vertex_Count = 3;
end is_Triangle;
function is_Clockwise (Self : in Polygon) return Boolean
is
i : constant Site := Self.Vertices (1);
j : constant Site := Self.Vertices (1);
k : constant Site := Self.Vertices (1);
z : Real := (j (1) - i (1))
* (k (2) - j (2));
begin
z := z - (j (2) - i (2))
* (k (1) - j (1));
return z < 0.0;
end is_Clockwise;
function is_Convex (Self : in Polygon) return Boolean
is
negative_Found,
positive_Found : Boolean := False;
begin
if is_Triangle (Self)
then
return True; -- All triangles are convex.
end if;
for i in 1 .. Self.Vertex_Count
loop
declare
k0 : constant Site := Self.Vertices (i);
function get_k1 return Site
is
begin
if i = Self.Vertex_Count
then return Self.Vertices (1);
else return Self.Vertices (i + 1);
end if;
end get_k1;
k1 : constant Site := get_k1;
function get_k2 return Site
is
begin
if i = Self.Vertex_Count - 1 then return Self.Vertices (1);
elsif i = Self.Vertex_Count then return Self.Vertices (2);
else return Self.Vertices (i + 2);
end if;
end get_k2;
k2 : constant Site := get_k2;
function get_Crossproduct return Real
is
dx1 : constant Real := k1 (1) - k0 (1);
dy1 : constant Real := k1 (2) - k0 (2);
dx2 : constant Real := k2 (1) - k1 (1);
dy2 : constant Real := k2 (2) - k1 (2);
begin
return dx1 * dy2 - dy1 * dx2;
end get_Crossproduct;
Crossproduct : constant Real := get_Crossproduct;
begin
if Crossproduct > 0.0
then
if negative_Found
then
return False;
end if;
positive_Found := True;
elsif Crossproduct < 0.0
then
if positive_Found
then
return False;
end if;
negative_Found := True;
end if;
end;
end loop;
return True;
end is_Convex;
function Area (Self : Polygon) return Real
is
Result : Real := 0.0;
begin
for i in 2 .. Self.Vertex_Count - 1
loop
Result := Result + Area (Triangle' (Vertices => [Self.Vertices (1),
Self.Vertices (i),
Self.Vertices (i + 1)]));
end loop;
return Result;
end Area;
function Perimeter (Self : Polygon) return Real
is
Result : Real := Distance (Self.Vertices (1),
Self.Vertices (Self.Vertex_Count));
begin
for i in 1 .. Self.Vertex_Count - 1
loop
Result := Result + Distance (Self.Vertices (i),
Self.Vertices (i + 1));
end loop;
return Result;
end Perimeter;
function Angle (Self : in Polygon; at_Vertex : in Positive) return Radians
is
Tri : constant Triangle := (vertices => [Self.Vertices (at_Vertex),
next_Vertex (Self, at_Vertex),
prior_Vertex (Self, at_Vertex)]);
begin
return Angle (Tri, 1);
end Angle;
function Image (Self : in Polygon) return String
is
pragma Unreferenced (Self);
begin
return "Polygon image (TODO)";
end Image;
end any_Math.any_Geometry.any_d2;

176
1-base/math/source/generic/pure/geometry/any_math-any_geometry-any_d2.ads Normal file
View File

@@ -0,0 +1,176 @@
generic
package any_Math.any_Geometry.any_d2
--
-- Provides a namespace and core types for 2D geometry.
--
is
pragma Pure;
---------
-- Sites
--
-- Cartesian
--
subtype Site is Vector_2; -- 2D cartesian coordinates.
type Sites is array (Positive range <>) of Site;
type Grid is array (Positive range <>,
Positive range <>) of Site;
function Distance (From, To : Site) return Real;
-- Polar
--
type polar_Site is -- 2D polar coordinates.
record
Angle : Radians;
Extent : Real;
end record;
function to_Polar (Self : in Site) return polar_Site;
function to_Site (Self : in polar_Site) return Site;
function Angle (Self : in Site) return Radians;
function Extent (Self : in Site) return Real;
---------
-- Lines
--
type Line is private;
function to_Line (Anchor : in Site;
Angle : in Radians) return Line;
function to_Line (Site_1,
Site_2 : in Site) return Line;
function X_of (Self : in Line; Y : in Real) return Real;
function Y_of (Self : in Line; X : in Real) return Real;
function Gradient (Self : in Line) return Real;
----------
-- Bounds
--
type bounding_Box is
record
Lower,
Upper : Site;
end record;
null_Bounds : constant bounding_Box;
function to_bounding_Box (Self : Sites) return bounding_Box;
function "or" (Left : in bounding_Box; Right : in Site) return bounding_Box;
--
-- Returns the bounds expanded to include the vector.
function "or" (Left : in bounding_Box; Right : in bounding_Box) return bounding_Box;
--
-- Returns the bounds expanded to include both Left and Right.
function "+" (Left : in bounding_Box; Right : in Vector_2) return bounding_Box;
--
-- Returns the bounds translated by the vector.
function Extent (Self : in bounding_Box; Dimension : in Index) return Real;
function Image (Self : in bounding_Box) return String;
----------
-- Circles
--
type Circle is
record
Radius : Real;
end record;
function Area (Self : Circle) return Real;
function Perimeter (Self : Circle) return Real;
-----------
-- Polygons
--
type Polygon (Vertex_Count : Positive) is
record
Vertices : Sites (1 .. Vertex_Count);
end record;
function Area (Self : in Polygon) return Real; -- Polygon must be convex.
function Perimeter (Self : in Polygon) return Real;
function Angle (Self : in Polygon; at_Vertex : in Positive) return Radians;
function is_Triangle (Self : in Polygon) return Boolean;
function is_Convex (Self : in Polygon) return Boolean;
function is_Clockwise (Self : in Polygon) return Boolean;
function Centroid (Self : in Polygon) return Site;
procedure center (Self : in out Polygon);
function prior_Vertex (Self : in Polygon; to_Vertex : in Positive) return Site;
function next_Vertex (Self : in Polygon; to_Vertex : in Positive) return Site;
function Image (Self : in Polygon) return String;
------------
-- Triangles
--
type Triangle is
record
Vertices : Sites (1 .. 3);
end record;
function Area (Self : in Triangle) return Real;
function Perimeter (Self : in Triangle) return Real;
function Angle (Self : in Triangle; at_Vertex : in Positive) return Radians;
function prior_Vertex (Self : in Triangle; to_Vertex : in Positive) return Site;
function next_Vertex (Self : in Triangle; to_Vertex : in Positive) return Site;
Degenerate,
Colinear : exception;
private
type Line_Format is (anchored_Gradient, two_Points);
type Line (Kind : Line_Format := Line_Format'First) is
record
case Kind is
when anchored_Gradient =>
Anchor : Site;
Gradient : Real;
when two_Points =>
Sites : any_d2.Sites (1 .. 2);
end case;
end record;
null_Bounds : constant bounding_Box := (lower => [Real'Last, Real'Last],
upper => [Real'First, Real'First]);
end any_Math.any_Geometry.any_d2;

546
1-base/math/source/generic/pure/geometry/any_math-any_geometry-any_d3-any_modeller-any_forge.adb Normal file
View File

@@ -0,0 +1,546 @@
with
ada.Text_IO,
ada.Strings.unbounded,
ada.Strings.Maps;
package body any_Math.any_Geometry.any_d3.any_Modeller.any_Forge
is
function to_Box_Model (half_Extents : in Vector_3 := [0.5, 0.5, 0.5]) return a_Model
is
pragma Unreferenced (half_Extents);
Modeller : any_Modeller.item;
begin
Modeller.add_Triangle ([0.0, 0.0, 0.0],
[1.0, 0.0, 0.0],
[1.0, 1.0, 0.0]);
Modeller.add_Triangle ([1.0, 1.0, 0.0],
[0.0, 1.0, 0.0],
[0.0, 0.0, 0.0]);
-- TODO: Add the rest.
return Modeller.Model;
end to_Box_Model;
function to_Capsule_Model (Length : in Real := 1.0;
Radius : in Real := 0.5) return a_Model
is
use Functions;
quality_Level : constant Positive := 4;
sides_Count : constant Positive := Positive (quality_Level * 4); -- Number of sides to the cylinder (divisible by 4).
type Edge is -- 'Barrel' edge.
record
Fore : Site;
Aft : Site;
end record;
type Edges is array (Positive range 1 .. sides_Count) of Edge;
type arch_Edges is array (Positive range 1 .. quality_Level) of Sites (1 .. sides_Count);
tmp,
ny, nz,
start_nx,
start_ny : Real;
a : constant Real := Pi * 2.0 / Real (sides_Count);
ca : constant Real := Cos (a);
sa : constant Real := Sin (a);
L : constant Real := Length * 0.5;
the_Edges : Edges;
Modeller : any_Modeller.item;
begin
-- Define cylinder body.
--
ny := 1.0;
nz := 0.0; -- Normal vector = (0, ny, nz)
for Each in Edges'Range
loop
the_Edges (Each).Fore (1) := ny * Radius;
the_Edges (Each).Fore (2) := nz * Radius;
the_Edges (Each).Fore (3) := L;
the_Edges (Each).Aft (1) := ny * Radius;
the_Edges (Each).Aft (2) := nz * Radius;
the_Edges (Each).Aft (3) := -L;
-- Rotate ny, nz.
--
tmp := ca * ny - sa * nz;
nz := sa * ny + ca * nz;
ny := tmp;
end loop;
for Each in Edges'Range
loop
if Each /= Edges'Last
then
Modeller.add_Triangle (the_Edges (Each) .Fore,
the_Edges (Each) .Aft,
the_Edges (Each + 1).Aft);
Modeller.add_Triangle (the_Edges (Each + 1).Aft,
the_Edges (Each + 1).Fore,
the_Edges (Each) .Fore);
else
Modeller.add_Triangle (the_Edges (Each) .Fore,
the_Edges (Each) .Aft,
the_Edges (edges'First).Aft);
Modeller.add_Triangle (the_Edges (edges'First).Aft,
the_Edges (edges'First).Fore,
the_Edges (Each) .Fore);
end if;
end loop;
-- Define fore cylinder cap.
--
declare
the_arch_Edges : arch_Edges;
begin
start_nx := 0.0;
start_ny := 1.0;
for each_Hoop in 1 .. quality_Level
loop
-- Get start_n2 = rotated start_n.
--
declare
start_nx2 : constant Real := ca * start_nx + sa * start_ny;
start_ny2 : constant Real := -sa * start_nx + ca * start_ny;
begin
-- Get n = start_n and n2 = start_n2.
--
ny := start_ny;
nz := 0.0;
declare
nx2 : constant Real := start_nx2;
ny2 : Real := start_ny2;
nz2 : Real := 0.0;
begin
for Each in 1 .. sides_Count
loop
the_arch_Edges (each_Hoop)(Each) (1) := ny2 * Radius;
the_arch_Edges (each_Hoop)(Each) (2) := nz2 * Radius;
the_arch_Edges (each_Hoop)(Each) (3) := nx2 * Radius + L;
-- Rotate n, n2.
--
tmp := ca * ny - sa * nz;
nz := sa * ny + ca * nz;
ny := tmp;
tmp := ca * ny2 - sa * nz2;
nz2 := sa * ny2 + ca * nz2;
ny2 := tmp;
end loop;
end;
start_nx := start_nx2;
start_ny := start_ny2;
end;
end loop;
for Each in 1 .. sides_Count
loop
if Each /= sides_Count
then
Modeller.add_Triangle (the_Edges (Each) .Fore,
the_Edges (Each + 1).Fore,
the_arch_Edges (1) (Each));
else
Modeller.add_Triangle (the_Edges (Each).Fore,
the_Edges (1) .Fore,
the_arch_Edges (1) (Each));
end if;
if Each /= sides_Count
then
Modeller.add_Triangle (the_Edges (Each + 1).Fore,
the_arch_Edges (1) (Each + 1),
the_arch_Edges (1) (Each));
else
Modeller.add_Triangle (the_Edges (1).Fore,
the_arch_Edges (1) (1),
the_arch_Edges (1) (Each));
end if;
end loop;
for each_Hoop in 1 .. quality_Level - 1
loop
for Each in 1 .. sides_Count
loop
declare
function next_Hoop_Vertex return Positive
is
begin
if Each = sides_Count then return 1;
else return Each + 1;
end if;
end next_Hoop_Vertex;
begin
Modeller.add_Triangle (the_arch_Edges (each_Hoop) (Each),
the_arch_Edges (each_Hoop) (next_Hoop_Vertex),
the_arch_Edges (each_Hoop + 1) (Each));
if each_Hoop /= quality_Level - 1
then
Modeller.add_Triangle (the_arch_Edges (each_Hoop) (next_Hoop_Vertex),
the_arch_Edges (each_Hoop + 1) (next_Hoop_Vertex),
the_arch_Edges (each_Hoop + 1) (Each));
end if;
end;
end loop;
end loop;
end;
-- Define aft cylinder cap.
--
declare
the_arch_Edges : arch_Edges;
begin
start_nx := 0.0;
start_ny := 1.0;
for each_Hoop in 1 .. quality_Level
loop
declare
-- Get start_n2 = rotated start_n.
--
start_nx2 : constant Real := ca * start_nx - sa * start_ny;
start_ny2 : constant Real := sa * start_nx + ca * start_ny;
begin
-- Get n = start_n and n2 = start_n2.
--
ny := start_ny;
nz := 0.0;
declare
nx2 : constant Real := start_nx2;
ny2 : Real := start_ny2;
nz2 : Real := 0.0;
begin
for Each in 1 .. sides_Count
loop
the_arch_Edges (each_Hoop) (Each) (1) := ny2 * Radius;
the_arch_Edges (each_Hoop) (Each) (2) := nz2 * Radius;
the_arch_Edges (each_Hoop) (Each) (3) := nx2 * Radius - L;
-- Rotate n, n2
--
tmp := ca * ny - sa * nz;
nz := sa * ny + ca * nz;
ny := tmp;
tmp := ca * ny2 - sa * nz2;
nz2 := sa * ny2 + ca * nz2;
ny2 := tmp;
end loop;
end;
start_nx := start_nx2;
start_ny := start_ny2;
end;
end loop;
for Each in 1 .. sides_Count
loop
if Each /= sides_Count
then
Modeller.add_Triangle (the_Edges (Each).Aft,
the_arch_Edges (1) (Each),
the_Edges (Each + 1).Aft);
else
Modeller.add_Triangle (the_Edges (Each).Aft,
the_arch_Edges (1) (Each),
the_Edges (1).Aft);
end if;
if Each /= sides_Count
then
Modeller.add_Triangle (The_Edges (Each + 1).Aft,
the_arch_Edges (1) (Each),
the_arch_Edges (1) (Each + 1));
else
Modeller.add_Triangle (the_Edges (1).Aft,
the_arch_Edges (1) (Each),
the_arch_Edges (1) (1));
end if;
end loop;
for each_Hoop in 1 .. quality_Level - 1
loop
for Each in 1 .. sides_Count
loop
declare
function next_Hoop_Vertex return Positive
is
begin
if Each = sides_Count then return 1;
else return Each + 1;
end if;
end next_hoop_Vertex;
begin
Modeller.add_Triangle (the_arch_Edges (each_Hoop) (Each),
the_arch_Edges (each_Hoop + 1) (Each),
the_arch_Edges (each_Hoop) (next_Hoop_Vertex));
if each_Hoop /= quality_Level - 1
then
Modeller.add_Triangle (the_arch_Edges (each_Hoop) (next_hoop_Vertex),
the_arch_Edges (each_Hoop + 1) (Each),
the_arch_Edges (each_Hoop + 1) (next_Hoop_Vertex));
end if;
end;
end loop;
end loop;
end;
return Modeller.Model;
end to_capsule_Model;
-- Polar to euclidian shape models.
--
function to_Radians (From : in Latitude) return Radians
is
begin
return Radians (From) * Pi / 180.0;
end to_Radians;
function to_Radians (From : in Longitude) return Radians
is
begin
return Radians (From) * Pi / 180.0;
end to_Radians;
function polar_Model_from (model_Filename : in String) return polar_Model -- TODO: Handle different file formats.
is
use Functions,
ada.Text_IO,
ada.Strings.unbounded;
the_File : File_type;
the_Text : unbounded_String;
begin
open (the_File, in_File, model_Filename);
while not end_of_File (the_File)
loop
append (the_Text, get_Line (the_File) & " ");
end loop;
declare
text_Length : constant Natural := Length (the_Text);
First : Positive := 1;
function get_Real return Real
is
use ada.Strings,
ada.Strings.Maps;
real_Set : constant Character_Set := to_Set (Span => (Low => '0',
High => '9'))
or to_Set ('-' & '.');
Last : Positive;
Result : Real;
begin
find_Token (the_Text, Set => real_Set,
From => First,
Test => Inside,
First => First,
Last => Last);
Result := Real'Value (Slice (the_Text,
Low => First,
High => Last));
First := Last + 1;
return Result;
end get_Real;
Lat : Latitude;
Long : Longitude;
Value : Integer;
Distance : Real;
Scale : constant Real := 10.0; -- TODO: Add a 'Scale' parameter.
the_Model : polar_Model;
begin
while First < text_Length
loop
Value := Integer (get_Real);
exit when Value = 360;
Long := Longitude (Value);
Lat := Latitude (get_Real);
Distance := get_Real;
the_Model (Long) (Lat).Site (1) := Scale * Distance * Cos (to_Radians (Lat)) * Sin (to_Radians (Long));
the_Model (Long) (Lat).Site (2) := Scale * Distance * Sin (to_Radians (Lat));
the_Model (Long) (Lat).Site (3) := Scale * Distance * Cos (to_Radians (Lat)) * Cos (to_Radians (Long));
end loop;
return the_Model;
end;
end polar_Model_from;
function mesh_Model_from (Model : in polar_Model) return a_Model
is
the_raw_Model : polar_Model := Model;
the_mesh_Model : a_Model (site_Count => 2522,
tri_Count => 73 * (16 * 4 + 6));
the_longitude : Longitude := 0;
the_latitude : Latitude ;
the_Vertex : Positive := 1;
the_Triangle : Positive := 1;
the_North_Pole : Positive;
the_South_Pole : Positive;
function Sum (the_Longitude : in Longitude; Increment : in Integer) return Longitude
is
Result : Integer := Integer (the_Longitude) + Increment;
begin
if Result >= 360
then
Result := Result - 360;
end if;
return longitude (Result);
end Sum;
begin
the_mesh_Model.Sites (the_Vertex) := (the_raw_model (0) (-90).Site);
the_North_Pole := the_Vertex;
the_raw_Model (0) (-90).Id := the_Vertex;
the_Vertex := the_Vertex + 1;
the_mesh_Model.Sites (the_Vertex) := (the_raw_model (0) (90).Site);
the_south_Pole := the_Vertex;
the_raw_Model (0) (90).Id := the_Vertex;
the_Vertex := the_Vertex + 1;
loop
the_latitude := -90;
loop
if the_Latitude = -90
then
the_raw_Model (the_Longitude) (the_Latitude).Id := the_North_Pole;
elsif the_Latitude = 90
then
the_raw_Model (the_Longitude) (the_Latitude).Id := the_South_Pole;
else
the_mesh_Model.Sites (the_Vertex) := the_raw_model (the_Longitude) (the_Latitude).Site;
the_raw_Model (the_Longitude) (the_Latitude).Id := the_Vertex;
the_Vertex := the_Vertex + 1;
end if;
exit when the_Latitude = 90;
the_Latitude := the_Latitude + 5;
end loop;
exit when the_Longitude = 355;
the_Longitude := the_Longitude + 5;
end loop;
the_Longitude := 0;
loop
the_mesh_Model.Triangles (the_Triangle) := [1 => the_North_Pole,
2 => the_raw_Model (Sum (the_Longitude, 5)) (-85).Id,
3 => the_raw_Model ( the_Longitude ) (-85).Id];
the_Triangle := the_Triangle + 1;
the_mesh_Model.Triangles (the_Triangle) := [1 => the_South_Pole,
2 => the_raw_Model (the_Longitude) (85).Id,
3 => the_raw_Model (Sum (the_Longitude, 5)) (85).Id];
the_Triangle := the_Triangle + 1;
the_Latitude := -85;
loop
the_mesh_Model.Triangles (the_Triangle) := [1 => the_raw_Model ( the_Longitude) (the_Latitude ).Id,
2 => the_raw_Model (Sum (the_Longitude, 5)) (the_Latitude ).Id,
3 => the_raw_Model ( the_Longitude) (the_Latitude + 5).Id];
the_Triangle := the_Triangle + 1;
the_mesh_Model.Triangles (the_Triangle) := [1 => the_raw_Model (the_Longitude) (the_Latitude + 5).Id,
2 => the_raw_Model (Sum (the_Longitude, 5)) (the_Latitude ).Id,
3 => the_raw_Model (Sum (the_Longitude, 5)) (the_Latitude + 5).Id];
the_Triangle := the_Triangle + 1;
the_Latitude := the_Latitude + 5;
exit when the_Latitude = 85;
end loop;
exit when the_Longitude = 355;
the_Longitude := the_Longitude + 5;
end loop;
the_mesh_Model.Triangles (the_Triangle) := [1 => the_North_Pole,
2 => the_raw_Model (5) (-85).Id,
3 => the_raw_Model (0) (-85).Id];
the_Triangle := the_Triangle + 1;
the_mesh_Model.Triangles (the_Triangle) := [1 => the_South_Pole,
2 => the_raw_Model (0) (85).Id,
3 => the_raw_Model (5) (85).Id];
the_Triangle := the_Triangle + 1;
the_latitude := -85;
loop
the_mesh_Model.Triangles (the_Triangle) := [1 => the_raw_Model (0) (the_Latitude ).Id,
2 => the_raw_Model (5) (the_Latitude ).Id,
3 => the_raw_Model (0) (the_Latitude + 5).Id];
the_Triangle := the_Triangle + 1;
the_mesh_Model.Triangles (the_Triangle) := [1 => the_raw_Model (0) (the_Latitude + 5).Id,
2 => the_raw_Model (5) (the_Latitude ).Id,
3 => the_raw_Model (5) (the_Latitude + 5).Id];
the_Triangle := the_Triangle + 1;
the_Latitude := the_Latitude + 5;
exit when the_Latitude = 85;
end loop;
return the_mesh_Model;
end mesh_Model_from;
end any_Math.any_Geometry.any_d3.any_Modeller.any_Forge;

42
1-base/math/source/generic/pure/geometry/any_math-any_geometry-any_d3-any_modeller-any_forge.ads Normal file
View File

@@ -0,0 +1,42 @@
generic
package any_Math.any_Geometry.any_d3.any_Modeller.any_Forge
--
-- Provides constructors for several geometry primitives.
--
is
function to_Box_Model (half_Extents : in Vector_3 := [0.5, 0.5, 0.5]) return a_Model;
function to_Capsule_Model (Length : in Real := 1.0;
Radius : in Real := 0.5) return a_Model;
type Latitude is range -90 .. 90;
type Longitude is range 0 .. 359;
no_Id : constant := Positive'Last;
type Vertex is
record
Id : Positive := no_Id;
Site : any_Geometry.any_d3.Site;
end record;
type longitude_Line is array (Latitude) of Vertex;
type polar_Model is array (Longitude) of longitude_Line;
type Vertices is array (Positive range <>) of Vertex;
type Triangle is array (Positive range 1 .. 3) of Positive;
type Triangles is array (Positive range <>) of Triangle;
-- type mesh_Model (num_Vertices : Positive;
-- num_Triangles : Positive) is
-- record
-- Vertices : mesh.Vertices (1 .. num_Vertices);
-- Triangles : mesh.Triangles (1 .. num_Triangles);
-- end record;
function polar_Model_from (Model_Filename : in String) return polar_model;
function mesh_Model_from (Model : in polar_Model) return a_Model; -- mesh_Model;
end any_Math.any_Geometry.any_d3.any_Modeller.any_Forge;

158
1-base/math/source/generic/pure/geometry/any_math-any_geometry-any_d3-any_modeller.adb Normal file
View File

@@ -0,0 +1,158 @@
with
ada.Strings.Hash;
package body any_Math.any_Geometry.any_d3.any_Modeller
is
use ada.Containers;
function Hash (Site : in my_Vertex) return ada.Containers.Hash_type
is
use ada.Strings;
begin
return Hash ( Site (1)'Image
& Site (2)'Image
& Site (3)'Image);
end Hash;
function demand_Index (Self : in out Item;
for_Vertex : in my_Vertex) return Natural
--
-- If the vertex exists in the map, return the associated index.
-- Otherwise add the new vertex and return it's index.
--
is
use Vertex_Maps_of_Index;
Cursor : constant Vertex_Maps_of_Index.Cursor := Self.Index_Map.find (for_Vertex);
begin
if has_Element (Cursor)
then
return Element (Cursor);
end if;
Self.Vertices.append (Vertex (for_Vertex));
declare
new_Index : constant Natural := Natural (Self.Vertices.Length);
begin
Self.Index_Map.insert (for_Vertex, new_Index);
return new_Index;
end;
end demand_Index;
function "<" (Left, Right : in Index_Triangle) return Boolean
is
begin
if Left (1) < Right (1) then return True; end if;
if Left (1) > Right (1) then return False; end if;
if Left (2) < Right (2) then return True; end if;
if Left (2) > Right (2) then return False; end if;
if Left (3) < Right (3) then return True; end if;
return False;
end "<";
procedure add_Triangle (Self : in out Item; Vertex_1, Vertex_2, Vertex_3 : in Site)
is
vertex_1_Index : constant Natural := demand_Index (Self, my_Vertex (Vertex_1));
vertex_2_Index : constant Natural := demand_Index (Self, my_Vertex (Vertex_2));
vertex_3_Index : constant Natural := demand_Index (Self, my_Vertex (Vertex_3));
new_Triangle : constant index_Triangle := [vertex_1_Index, vertex_2_Index, vertex_3_Index];
new_Triangle_rotated_1 : constant index_Triangle := [vertex_3_Index, vertex_1_Index, vertex_2_Index];
new_Triangle_rotated_2 : constant index_Triangle := [vertex_2_Index, vertex_3_Index, vertex_1_Index];
begin
if new_Triangle (1) = new_Triangle (2)
or else new_Triangle (1) = new_Triangle (3)
or else new_Triangle (2) = new_Triangle (3)
then
null; -- Discard collapsed triangle.
else
if Self.Triangles.contains (new_triangle)
or else Self.Triangles.contains (new_triangle_rotated_1)
or else Self.Triangles.contains (new_triangle_rotated_2)
then
null; -- Triangle is already present.
else
Self.Triangles.include (new_Triangle);
end if;
end if;
end add_Triangle;
procedure clear (Self : in out Item)
is
begin
Self.Triangles.clear;
Self.Vertices .clear;
Self.Index_Map.clear;
end clear;
function Triangle_Count (Self : in Item) return Natural
is
begin
return Natural (Self.Triangles.Length);
end triangle_Count;
function Model (Self : in Item) return a_Model
is
Result : a_Model := (Site_Count => Integer (Self.Vertices.Length),
Tri_Count => Integer (Self.Triangles.Length),
Sites => <>,
Triangles => <>);
begin
for i in 1 .. Index (Result.site_Count)
loop
Result.Sites (i) := Self.Vertices.Element (i);
end loop;
declare
use Index_Triangle_Sets;
Cursor : Index_Triangle_Sets.Cursor := Self.Triangles.First;
begin
for i in 1 .. Result.Tri_Count
loop
Result.Triangles (i) := Element (Cursor);
next (Cursor);
end loop;
end;
return Result;
end Model;
function bounding_Sphere_Radius (Self : in out Item) return Real
is
use Functions;
begin
if Self.bounding_Sphere_Radius = Real'First
then
for Each of Self.Vertices
loop
Self.bounding_sphere_Radius := Real'Max (Self.bounding_sphere_Radius,
SqRt ( Each (1) * Each (1)
+ Each (2) * Each (2)
+ Each (3) * Each (3)));
end loop;
end if;
return Self.bounding_sphere_Radius;
end bounding_sphere_Radius;
end any_Math.any_Geometry.any_d3.any_Modeller;

81
1-base/math/source/generic/pure/geometry/any_math-any_geometry-any_d3-any_modeller.ads Normal file
View File

@@ -0,0 +1,81 @@
private
with
ada.Containers.Vectors,
ada.Containers.hashed_Maps,
ada.Containers.ordered_Sets;
generic
package any_Math.any_Geometry.any_d3.any_Modeller
is
type Item is tagged private;
type View is access all Item;
--------------
-- Attributes
--
procedure add_Triangle (Self : in out Item; Vertex_1,
Vertex_2,
Vertex_3 : in Site);
function Triangle_Count (Self : in Item) return Natural;
function Model (Self : in Item) return a_Model;
function bounding_Sphere_Radius (Self : in out Item) return Real;
--
-- Caches the radius on 1st call.
--------------
-- Operations
--
procedure clear (Self : in out Item);
private
subtype Vertex is Site;
type my_Vertex is new Vertex;
--------------
-- Containers
--
function Hash (Site : in my_Vertex) return ada.Containers.Hash_type;
package Vertex_Maps_of_Index is new ada.Containers.hashed_Maps (my_Vertex,
Natural,
Hash,
"=");
subtype Vertex_Map_of_Index is Vertex_Maps_of_Index.Map;
package Vertex_Vectors is new Ada.Containers.Vectors (Positive, Vertex);
subtype Vertex_Vector is Vertex_Vectors.Vector;
subtype Index_Triangle is any_Geometry.Triangle;
function "<" (Left, Right : in Index_Triangle) return Boolean;
package Index_Triangle_Sets is new ada.Containers.ordered_Sets (Element_Type => Index_Triangle,
"<" => "<",
"=" => "=");
subtype Index_Triangle_Set is Index_Triangle_Sets.Set;
------------
-- Modeller
--
type Item is tagged
record
Triangles : Index_Triangle_Set;
Vertices : Vertex_Vector;
Index_Map : Vertex_Map_of_Index;
bounding_Sphere_Radius : Real := Real'First;
end record;
end any_Math.any_Geometry.any_d3.any_Modeller;

119
1-base/math/source/generic/pure/geometry/any_math-any_geometry-any_d3.adb Normal file
View File

@@ -0,0 +1,119 @@
package body any_Math.any_Geometry.any_d3
is
--------
-- Plane
--
procedure normalise (the_Plane : in out Plane)
is
use Functions;
inverse_Magnitude : constant Real := 1.0 / SqRt ( the_Plane (1) * the_Plane (1)
+ the_Plane (2) * the_Plane (2)
+ the_Plane (3) * the_Plane (3));
begin
the_Plane (1) := the_Plane (1) * inverse_Magnitude;
the_Plane (2) := the_Plane (2) * inverse_Magnitude;
the_Plane (3) := the_Plane (3) * inverse_Magnitude;
the_Plane (4) := the_Plane (4) * inverse_Magnitude;
end normalise;
function Image (the_Model : in a_Model) return String
is
begin
return
"(Site_Count =>" & Integer'Image (the_Model.Site_Count) & ","
& " Tri_Count =>" & Integer'Image (the_Model. Tri_Count) & ")";
exception
when others =>
return "<TODO>";
end Image;
----------
-- Bounds
--
function to_bounding_Box (Self : Sites) return bounding_Box
is
Bounds : bounding_Box := null_Bounds;
begin
for Each in Self'Range
loop
Bounds.Lower (1) := Real'Min (Bounds.Lower (1), Self (Each)(1));
Bounds.Lower (2) := Real'Min (Bounds.Lower (2), Self (Each)(2));
Bounds.Lower (3) := Real'Min (Bounds.Lower (3), Self (Each)(3));
Bounds.Upper (1) := Real'Max (Bounds.Upper (1), Self (Each)(1));
Bounds.Upper (2) := Real'Max (Bounds.Upper (2), Self (Each)(2));
Bounds.Upper (3) := Real'Max (Bounds.Upper (3), Self (Each)(3));
end loop;
return Bounds;
end to_bounding_Box;
function Extent (Self : in bounding_Box; Dimension : in Index) return Real
is
begin
return Self.Upper (Dimension) - Self.Lower (Dimension);
end Extent;
function "or" (Left : in bounding_Box; Right : in Site) return bounding_Box
is
Result : bounding_Box;
begin
for i in Right'Range
loop
if Right (i) < Left.Lower (i)
then Result.Lower (i) := Right (i);
else Result.Lower (i) := Left.Lower (i);
end if;
if Right (i) > Left.Upper (i)
then Result.Upper (i) := Right (i);
else Result.Upper (i) := Left.Upper (i);
end if;
end loop;
return Result;
end "or";
function "or" (Left : in bounding_Box;
Right : in bounding_Box) return bounding_Box
is
Result : bounding_Box := Left or Right.Lower;
begin
Result := Result or Right.Upper;
return Result;
end "or";
function "+" (Left : in bounding_Box; Right : in Vector_3) return bounding_Box
is
begin
return (Left.Lower + Right,
Left.Upper + Right);
end "+";
function Image (Self : bounding_Box) return String
is
begin
return "(lower => " & Image (Self.Lower)
& ", upper => " & Image (Self.Upper) & ")";
end Image;
end any_Math.any_Geometry.any_d3;

78
1-base/math/source/generic/pure/geometry/any_math-any_geometry-any_d3.ads Normal file
View File

@@ -0,0 +1,78 @@
generic
package any_Math.any_Geometry.any_d3
--
-- Provides a namespace and core types for 3D geometry.
--
is
pragma Pure;
--------------
-- Core Types
--
subtype Site is Vector_3;
type Sites is array (Positive range <>) of Site;
type a_Model (Site_Count : Positive;
Tri_Count : Positive) is
record
Sites : any_d3 .Sites (1 .. Site_Count);
Triangles : any_Geometry.Triangles (1 .. Tri_Count);
end record;
function Image (the_Model : in a_Model) return String;
---------
-- Planes
--
type Plane is new Vector_4; -- A general plane equation.
procedure normalise (the_Plane : in out Plane);
----------
-- Bounds
--
type bounding_Box is
record
Lower,
Upper : Site;
end record;
null_Bounds : constant bounding_Box;
function to_bounding_Box (Self : Sites) return bounding_Box;
function "or" (Left : in bounding_Box; Right : in Site) return bounding_Box;
--
-- Returns the bounds expanded to include the vector.
function "or" (Left : in bounding_Box; Right : in bounding_Box) return bounding_Box;
--
-- Returns the bounds expanded to include both Left and Right.
function "+" (Left : in bounding_Box; Right : in Vector_3) return bounding_Box;
--
-- Returns the bounds translated by the vector.
function Extent (Self : in bounding_Box; Dimension : in Index) return Real;
function Image (Self : in bounding_Box) return String;
private
null_Bounds : constant bounding_Box := (Lower => [Real'Last, Real'Last, Real'Last],
Upper => [Real'First, Real'First, Real'First]);
end any_Math.any_Geometry.any_d3;

59
1-base/math/source/generic/pure/geometry/any_math-any_geometry.adb Normal file
View File

@@ -0,0 +1,59 @@
package body any_Math.any_Geometry
is
function Image (Self : in Triangle) return String
is
begin
return "(" & Vertex_Id'Image (Self (1)) & ","
& Vertex_Id'Image (Self (2)) & ","
& Vertex_Id'Image (Self (3)) & ")";
end Image;
function Image (Self : in Triangles) return String
is
Result : String (1 .. 1024);
Last : Standard.Natural := 0;
begin
for Each in Self'Range
loop
declare
Id_Image : constant String := Image (Self (Each));
begin
Result (Last + 1 .. Last + Id_Image'Length) := Id_Image;
Last := Last + Id_Image'Length;
end;
end loop;
return Result (1 .. Last);
exception
when Constraint_Error =>
declare
Ellipsis : constant String := " ...";
begin
Result (Result'Last - ellipsis'Length + 1 .. Result'Last) := ellipsis;
return Result (1 .. Last);
end;
end Image;
function Image (Self : in Model) return String
is
begin
return Self.Triangles.Image;
end Image;
function Image (Self : in Model_Triangles) return String
is
begin
return "Triangle_Count =>" & standard.Positive'Image (Self.Triangle_Count)
& Image (Self.Triangles);
end Image;
end any_Math.any_Geometry;

66
1-base/math/source/generic/pure/geometry/any_math-any_geometry.ads Normal file
View File

@@ -0,0 +1,66 @@
generic
package any_Math.any_Geometry
--
-- Provides a namespace and core types for geometry.
--
is
pragma Pure;
subtype Vertex_Id is Index;
type Vertex_Ids is array (Index range <>) of Vertex_Id;
subtype Triangle is Vertex_Ids (1 .. 3);
type Triangles is array (Index range <>) of Triangle;
function Image (Self : in Triangle) return String;
function Image (Self : in Triangles) return String;
--------
-- Model
--
type Model_Options is tagged null record;
default_Model_Options : constant Model_Options;
type Model_Triangles (Triangle_Count : Index) is tagged
record
Triangles : any_Geometry.Triangles (1 .. Triangle_Count);
end record;
function Image (Self : in Model_Triangles) return String;
type Model is abstract tagged
record
Triangles : access Model_Triangles'Class;
end record;
function Image (Self : in Model) return String;
----------------
-- Geometry Item
--
type Item is abstract tagged private;
procedure destroy (Self : in out Item) is abstract;
procedure expand (Self : access Item; By : in Real) is abstract;
private
type Item is abstract tagged
record
null;
end record;
default_Model_Options : constant Model_Options := (others => <>);
end any_Math.any_Geometry;

30
1-base/math/source/generic/pure/geometry/trigonometry/any_math-any_fast_rotation.adb Normal file
View File

@@ -0,0 +1,30 @@
with
cached_Rotation;
package body any_Math.any_fast_Rotation
is
function to_Matrix_2x2 (m11, m12,
m21, m22 : Real) return Matrix_2x2
is
begin
return (1 => (m11, m12),
2 => (m21, m22));
end to_Matrix_2x2;
package the_Cache is new cached_Rotation (Float_type => any_Math.Real,
Matrix_2x2_type => any_Math.Matrix_2x2,
float_elementary_Functions => any_math.Functions,
to_Matrix_2x2 => to_Matrix_2x2,
slot_Count => 10_000);
function to_Rotation (Angle : in Real) return access constant Matrix_2x2
is
begin
return the_Cache.to_Rotation (Angle);
end to_Rotation;
end any_Math.any_fast_Rotation;

11
1-base/math/source/generic/pure/geometry/trigonometry/any_math-any_fast_rotation.ads Normal file
View File

@@ -0,0 +1,11 @@
generic
package any_Math.any_fast_Rotation
is
function to_Rotation (Angle : in Real) return access constant Matrix_2x2;
private
pragma Inline_Always (to_Rotation);
end any_Math.any_fast_Rotation;

10
1-base/math/source/generic/pure/geometry/trigonometry/any_math-any_fast_trigonometry.ads Normal file
View File

@@ -0,0 +1,10 @@
with
cached_Trigonometry;
generic
package any_math.any_fast_Trigonometry
is
package Default is new cached_Trigonometry (Float_type => any_Math.Real,
slot_Count => 10_000);
end any_math.any_fast_Trigonometry;

43
1-base/math/source/generic/pure/geometry/trigonometry/cached_rotation.adb Normal file
View File

@@ -0,0 +1,43 @@
package body cached_Rotation
is
use ada.Numerics,
float_elementary_Functions;
the_Cache : array (0 .. slot_Count - 1) of aliased Matrix_2x2_type;
Pi_x_2 : constant := Pi * 2.0;
last_slot_Index : constant Float_type := Float_type (slot_Count - 1);
index_Factor : constant Float_type := last_slot_Index / Pi_x_2;
function to_Rotation (Angle : in Float_type) return access constant Matrix_2x2_type
is
the_Index : standard.Integer := standard.Integer (Angle * index_Factor) mod slot_Count;
begin
if the_Index < 0
then
the_index := the_Index + slot_Count;
end if;
return the_Cache (the_Index)'Access;
end to_Rotation;
begin
for Each in the_Cache'Range
loop
declare
Angle : constant Float_type := ( Float_type (Each) / Float_type (slot_Count - 1)
* Pi_x_2);
C : constant Float_type := Cos (Angle);
S : constant Float_type := Sin (Angle);
begin
the_Cache (Each) := to_Matrix_2x2 (C, -S,
S, C);
end;
end loop;
end cached_Rotation;

30
1-base/math/source/generic/pure/geometry/trigonometry/cached_rotation.ads Normal file
View File

@@ -0,0 +1,30 @@
with
ada.Numerics.generic_elementary_Functions;
generic
type Float_type is digits <>;
type Matrix_2x2_type is private;
with package float_elementary_Functions is new ada.Numerics.generic_elementary_Functions (Float_type);
with function to_Matrix_2x2 (m11, m12,
m21, m22 : Float_type) return Matrix_2x2_type;
slot_Count : Standard.Positive;
package cached_Rotation
--
-- Caches 2x2 rotation matrices of angles for speed at the cost of precision.
--
is
pragma Optimize (Time);
function to_Rotation (Angle : in Float_type) return access constant Matrix_2x2_type;
private
pragma Inline_Always (to_Rotation);
end cached_Rotation;

75
1-base/math/source/generic/pure/geometry/trigonometry/cached_trigonometry.adb Normal file
View File

@@ -0,0 +1,75 @@
with
ada.Numerics.generic_elementary_Functions;
package body cached_Trigonometry
is
Sin_Cache : array (0 .. slot_Count - 1) of Float_Type;
Cos_Cache : array (0 .. slot_Count - 1) of Float_Type;
Pi_x_2 : constant := ada.Numerics.Pi * 2.0;
last_slot_Index : constant Float_Type := Float_Type (slot_Count - 1);
index_Factor : constant Float_Type := last_slot_Index / Pi_x_2;
function Cos (Angle : in Float_Type) return Float_Type
is
Index : standard.Integer := standard.Integer (Angle * index_Factor) mod slot_Count;
begin
if Index < 0 then
Index := Index + slot_Count;
end if;
return Cos_Cache (Index);
end Cos;
function Sin (Angle : in Float_Type) return Float_Type
is
Index : standard.Integer := standard.Integer (Angle * index_Factor) mod slot_Count;
begin
if Index < 0 then
Index := Index + slot_Count;
end if;
return Sin_Cache (Index);
end Sin;
procedure get (Angle : in Float_Type; the_Cos : out Float_Type;
the_Sin : out Float_Type)
is
Index : standard.Integer := standard.Integer (Angle * index_Factor) mod slot_Count;
begin
if Index < 0 then
Index := Index + slot_Count;
end if;
the_Sin := Sin_Cache (Index);
the_Cos := Cos_Cache (Index);
end get;
-- TODO: Tan, arcCos, etc
package Functions is new Ada.Numerics.generic_elementary_Functions (Float_Type);
begin
for Each in cos_Cache'Range
loop
cos_Cache (Each) := Functions.cos ( Float_Type (Each) / Float_Type (slot_Count - 1)
* Pi_x_2);
end loop;
for Each in sin_Cache'Range
loop
sin_Cache (Each) := Functions.sin ( Float_Type (Each) / Float_Type (slot_Count - 1)
* Pi_x_2);
end loop;
end cached_Trigonometry;

29
1-base/math/source/generic/pure/geometry/trigonometry/cached_trigonometry.ads Normal file
View File

@@ -0,0 +1,29 @@
generic
type Float_type is digits <>;
slot_Count : standard.Positive;
package cached_Trigonometry
--
-- Caches trig functions for speed at the cost of precision.
--
is
pragma Optimize (Time);
function Cos (Angle : in Float_type) return Float_type;
function Sin (Angle : in Float_type) return Float_type;
procedure get (Angle : in Float_type; the_Cos : out Float_type;
the_Sin : out Float_type);
-- TODO: tan, arccos, etc
private
pragma Inline_Always (Cos);
pragma Inline_Always (Sin);
pragma Inline_Always (Get);
end cached_Trigonometry;