Source file : gl-math.adb
-------------------------------------------------------------------------
-- GL.Math
--
-- Copyright (c) Rod Kay 2016
-- AUSTRALIA
--
-- Permission granted to use this software, without any warranty,
-- for any purpose, provided this copyright note remains attached
-- and unmodified if sources are distributed further.
-------------------------------------------------------------------------
package body GL.Math is
use GL, REF;
-------------
-- Vectors --
-------------
function "*"(l: Double; v: Double_Vector_3D) return Double_Vector_3D is
begin
return (l*v(0),l*v(1),l*v(2));
end "*";
function "*"(v: Double_Vector_3D; l: Double) return Double_Vector_3D is
begin
return (l*v(0),l*v(1),l*v(2));
end "*";
function "+"(a,b: Double_Vector_3D) return Double_Vector_3D is
begin
return (a(0)+b(0),a(1)+b(1),a(2)+b(2));
end "+";
function "-"(a: Double_Vector_3D) return Double_Vector_3D is
begin
return (-a(0),-a(1),-a(2));
end "-";
function "-"(a,b: Double_Vector_3D) return Double_Vector_3D is
begin
return (a(0)-b(0),a(1)-b(1),a(2)-b(2));
end "-";
function "*"(a,b: Double_Vector_3D) return Double is -- Dot product.
begin
return a(0)*b(0)+a(1)*b(1)+a(2)*b(2);
end "*";
function "*"(a,b: Double_Vector_3D) return Double_Vector_3D is -- Cross product.
begin
return ( a(1)*b(2) - a(2)*b(1),
a(2)*b(0) - a(0)*b(2),
a(0)*b(1) - a(1)*b(0) );
end "*";
function Norm(a: Double_Vector_3D) return Double is
begin
return Sqrt(a(0)*a(0)+a(1)*a(1)+a(2)*a(2));
end Norm;
function Norm2(a: Double_Vector_3D) return Double is
begin
return a(0)*a(0)+a(1)*a(1)+a(2)*a(2);
end Norm2;
function Normalized(a: Double_Vector_3D) return Double_Vector_3D is
begin
return a * (1.0 / Norm(a));
end Normalized;
function Identical(a, b: Double_Vector_3D) return Boolean is
begin
return
Almost_zero(a(0)-b(0)) and then
Almost_zero(a(1)-b(1)) and then
Almost_zero(a(2)-b(2));
end Identical;
function Identical(a, b: RGB_Color) return Boolean is
begin
return
Almost_zero(a.red - b.red) and then
Almost_zero(a.green - b.green) and then
Almost_zero(a.blue - b.blue);
end Identical;
function Identical(a, b: RGBA_Color) return Boolean is
begin
return
Almost_zero(a.red - b.red) and then
Almost_zero(a.green - b.green) and then
Almost_zero(a.blue - b.blue) and then
Almost_zero(a.alpha - b.alpha);
end Identical;
-- Angles
--
function Angle (Point_1, Point_2, Point_3 : Double_Vector_3D) return Double
is
Vector_1 : constant Double_Vector_3D := Normalized (Point_1 - Point_2);
Vector_2 : constant Double_Vector_3D := Normalized (Point_3 - Point_2);
Cos_Theta : constant Double := Vector_1 * Vector_2;
begin
if Cos_Theta <= -1.0 then
return Ada.Numerics.Pi;
elsif Cos_Theta >= 1.0 then
return 0.0;
else
return Arccos (Cos_Theta);
end if;
end Angle;
function to_Degrees (Radians : Double) return Double
is
use Ada.Numerics;
begin
return Radians * 180.0 / Pi;
end to_Degrees;
function to_Radians (Degrees : Double) return Double
is
use Ada.Numerics;
begin
return Degrees * Pi / 180.0;
end to_Radians;
--------------
-- Matrices --
--------------
function "*"(A,B: Matrix_33) return Matrix_33 is
r: Double; AB: Matrix_33;
begin
for i in 1..3 loop
for j in 1..3 loop
r:= 0.0;
for k in 1..3 loop
r:= r + (A(i,k) * B(k,j));
end loop;
AB(i,j):= r;
end loop;
end loop;
return AB;
end "*";
function "*"(A: Matrix_33; x: Double_Vector_3D) return Double_Vector_3D is
r: Double;
Ax: Double_Vector_3D;
-- banana skin: Matrix has range 1..3, Vector 0..2 (GL)
begin
for i in 1..3 loop
r:= 0.0;
for j in 1..3 loop
r:= r + A(i,j) * x(j-1);
end loop;
Ax(i-1):= r;
end loop;
return Ax;
end "*";
function "*"(A: Matrix_44; x: Double_Vector_3D) return Double_Vector_3D is
r: Double;
type Vector_4D is array (0..3) of GL.Double;
x_4D : constant Vector_4D := (x(0), x(1), x(2) , 1.0);
Ax : Vector_4D;
-- banana skin: Matrix has range 1..3, Vector 0..2 (GL)
begin
for i in 1..4 loop
r:= 0.0;
for j in 1..4 loop
r:= r + A(i,j) * x_4D(j-1);
end loop;
Ax(i-1):= r;
end loop;
return (Ax (0), Ax (1), Ax(2));
end "*";
function "*"(A: Matrix_44; x: Double_Vector_3D) return Vector_4D is
r: Double;
x_4D : constant Vector_4D := (x(0), x(1), x(2) , 1.0);
Ax : Vector_4D;
-- banana skin: Matrix has range 1..3, Vector 0..2 (GL)
begin
for i in 1..4 loop
r:= 0.0;
for j in 1..4 loop
r:= r + A(i,j) * x_4D(j-1);
end loop;
Ax(i-1):= r;
end loop;
return Ax;
end "*";
function Transpose(A: Matrix_33) return Matrix_33 is
begin
return ( (A(1,1),A(2,1),A(3,1)),
(A(1,2),A(2,2),A(3,2)),
(A(1,3),A(2,3),A(3,3)));
end Transpose;
function Transpose(A: Matrix_44) return Matrix_44 is
begin
return ( (A(1,1),A(2,1),A(3,1),A(4,1)),
(A(1,2),A(2,2),A(3,2),A(4,2)),
(A(1,3),A(2,3),A(3,3),A(4,3)),
(A(1,4),A(2,4),A(3,4),A(4,4)));
end Transpose;
function Det(A: Matrix_33) return Double is
begin
return
A(1,1) * A(2,2) * A(3,3) +
A(2,1) * A(3,2) * A(1,3) +
A(3,1) * A(1,2) * A(2,3) -
A(3,1) * A(2,2) * A(1,3) -
A(2,1) * A(1,2) * A(3,3) -
A(1,1) * A(3,2) * A(2,3);
end Det;
function XYZ_rotation(ax,ay,az: Double) return Matrix_33 is
Mx, My, Mz: Matrix_33; c,s: Double;
begin
-- Around X.
c:= Cos( ax );
s:= Sin( ax );
Mx:= ( (1.0,0.0,0.0),
(0.0, c, -s),
(0.0, s, c) );
-- Around Y.
c:= Cos( ay );
s:= Sin( ay );
My:= ( ( c,0.0, -s),
(0.0,1.0,0.0),
( s,0.0, c) );
-- Around Z.
c:= Cos( az );
s:= Sin( az );
Mz:= ( ( c, -s,0.0),
( s, c,0.0),
(0.0,0.0,1.0) );
return Mz * My * Mx;
end XYZ_rotation;
function XYZ_rotation(v: Double_Vector_3D) return Matrix_33 is
begin
return XYZ_rotation(v(0),v(1),v(2));
end XYZ_rotation;
function Look_at(direction: Double_Vector_3D) return Matrix_33 is
v1, v2, v3: Double_Vector_3D;
begin
-- GL's look direction is the 3rd dimension (z).
v3:= Normalized(direction);
v2:= Normalized((v3(2),0.0,-v3(0)));
v1:= v2 * v3;
return
(((v1(0),v2(0),v3(0)),
(v1(1),v2(1),v3(1)),
(v1(2),v2(2),v3(2))
));
end Look_at;
function sub_Matrix (Self : in Matrix; start_Row, end_Row : in Positive;
start_Col, end_Col : in Positive) return Matrix
is
the_sub_Matrix : Matrix (1 .. end_Row - start_Row + 1,
1 .. end_Col - start_Col + 1);
begin
for Row in the_sub_Matrix'Range (1) loop
for Col in the_sub_Matrix'Range (2) loop
the_sub_Matrix (Row, Col) := Self (Row + start_Row - 1,
Col + start_Col - 1);
end loop;
end loop;
return the_sub_Matrix;
end sub_Matrix;
function Look_at (eye, center, up : Double_Vector_3D) return Matrix_33
is
forward : constant Double_Vector_3D := Normalized ((center (0) - eye (0), center (1) - eye (1), center (2) - eye (2)));
side : constant Double_Vector_3D := Normalized (forward * up);
new_up : constant Double_Vector_3D := side * forward;
begin
return (( side (0), side (1), side (2)),
( new_up (0), new_up (1), new_up (2)),
(-forward (0), -forward (1), -forward (2)));
end Look_at;
-- Following procedure is from Project Spandex, by Paul Nettle.
--
procedure Re_Orthonormalize(M: in out Matrix_33) is
dot1,dot2,vlen: Double;
begin
dot1:= M(1,1) * M(2,1) + M(1,2) * M(2,2) + M(1,3) * M(2,3);
dot2:= M(1,1) * M(3,1) + M(1,2) * M(3,2) + M(1,3) * M(3,3);
M(1,1) := M(1,1) - dot1 * M(2,1) - dot2 * M(3,1);
M(1,2) := M(1,2) - dot1 * M(2,2) - dot2 * M(3,2);
M(1,3) := M(1,3) - dot1 * M(2,3) - dot2 * M(3,3);
vlen:= 1.0 / Sqrt(M(1,1) * M(1,1) +
M(1,2) * M(1,2) +
M(1,3) * M(1,3));
M(1,1):= M(1,1) * vlen;
M(1,2):= M(1,2) * vlen;
M(1,3):= M(1,3) * vlen;
dot1:= M(2,1) * M(1,1) + M(2,2) * M(1,2) + M(2,3) * M(1,3);
dot2:= M(2,1) * M(3,1) + M(2,2) * M(3,2) + M(2,3) * M(3,3);
M(2,1) := M(2,1) - dot1 * M(1,1) - dot2 * M(3,1);
M(2,2) := M(2,2) - dot1 * M(1,2) - dot2 * M(3,2);
M(2,3) := M(2,3) - dot1 * M(1,3) - dot2 * M(3,3);
vlen:= 1.0 / Sqrt(M(2,1) * M(2,1) +
M(2,2) * M(2,2) +
M(2,3) * M(2,3));
M(2,1):= M(2,1) * vlen;
M(2,2):= M(2,2) * vlen;
M(2,3):= M(2,3) * vlen;
M(3,1):= M(1,2) * M(2,3) - M(1,3) * M(2,2);
M(3,2):= M(1,3) * M(2,1) - M(1,1) * M(2,3);
M(3,3):= M(1,1) * M(2,2) - M(1,2) * M(2,1);
end Re_Orthonormalize;
-- type Matrix_44 is array(0..3,0..3) of aliased Double; -- for GL.MultMatrix
-- pragma Convention(Fortran, Matrix_44); -- GL stores matrices columnwise
M: Matrix_44;
-- M is a global variable for a clean 'Access and for setting once 4th dim.
pM: constant GL.doublePtr:= M (1, 1)'Unchecked_Access;
procedure Multiply_GL_Matrix( A: Matrix_33 ) is
begin
for i in 1..3 loop
for j in 1..3 loop
M(i,j):= A(i,j);
-- Funny deformations...
-- if j=2 then
-- M(i-1,j-1):= 0.5 * A(i,j);
-- end if;
end loop;
end loop;
GL.MultMatrixd(pM);
end Multiply_GL_Matrix;
procedure Set_GL_Matrix( A: Matrix_33 ) is
begin
GL.LoadIdentity;
Multiply_GL_Matrix( A );
end Set_GL_Matrix;
-- Ada 95 Quality and Style Guide, 7.2.7:
-- Tests for
--
-- (1) absolute "equality" to 0 in storage,
-- (2) absolute "equality" to 0 in computation,
-- (3) relative "equality" to 0 in storage, and
-- (4) relative "equality" to 0 in computation:
--
-- abs X <= Float_Type'Model_Small -- (1)
-- abs X <= Float_Type'Base'Model_Small -- (2)
-- abs X <= abs X * Float_Type'Model_Epsilon -- (3)
-- abs X <= abs X * Float_Type'Base'Model_Epsilon -- (4)
function Almost_zero(x: Double) return Boolean is
begin
return abs x <= Double'Base'Model_Small;
end Almost_zero;
function Almost_zero(x: GL.Float) return Boolean is
begin
return abs x <= GL.Float'Base'Model_Small;
end Almost_zero;
begin
for i in 1..3 loop
M(i,4):= 0.0;
M(4,i):= 0.0;
end loop;
M(4,4):= 1.0;
end GL.Math;
GLOBE_3D: Ada library for real-time 3D rendering.
Ada programming.