1. package body GL.Math is 
  2.  
  3.    use Ada.Numerics; 
  4.    use GL, REF; 
  5.  
  6.    ------------- 
  7.    -- Vectors -- 
  8.    ------------- 
  9.  
  10.    function "*" (l : Double; v : Double_Vector_3D) return Double_Vector_3D is (l * v (0), l * v (1), l * v (2)); 
  11.  
  12.    function "*" (v : Double_Vector_3D; l : Double) return Double_Vector_3D is (l * v (0), l * v (1), l * v (2)); 
  13.  
  14.    function "+" (a, b : Double_Vector_3D) return Double_Vector_3D is (a (0) + b (0), a (1) + b (1), a (2) + b (2)); 
  15.  
  16.    function "-" (a : Double_Vector_3D) return Double_Vector_3D is (-a (0), -a (1), -a (2)); 
  17.  
  18.    function "-" (a, b : Double_Vector_3D) return Double_Vector_3D is (a (0) - b (0), a (1) - b (1), a (2) - b (2)); 
  19.  
  20.    function "*" (a, b : Double_Vector_3D) return Double is (a (0) * b (0) + a (1) * b (1) + a (2) * b (2)); 
  21.  
  22.    function "*" (a, b : Double_Vector_3D) return Double_Vector_3D is (a (1) * b (2) - a (2) * b (1), 
  23.                                                                       a (2) * b (0) - a (0) * b (2), 
  24.                                                                       a (0) * b (1) - a (1) * b (0)); 
  25.  
  26.    function Norm (a : Double_Vector_3D) return Double is (Sqrt (a (0) * a (0) + a (1) * a (1) + a (2) * a (2))); 
  27.  
  28.    function Norm2 (a : Double_Vector_3D) return Double is (a (0) * a (0) + a (1) * a (1) + a (2) * a (2)); 
  29.  
  30.    function Normalized (a : Double_Vector_3D) return Double_Vector_3D is (a * (1.0 / Norm (a))); 
  31.  
  32.    -- Angles 
  33.    -- 
  34.  
  35.    function Angle (Point_1, Point_2, Point_3  : Double_Vector_3D) return Double is 
  36.  
  37.       Vector_1   : constant Double_Vector_3D := Normalized (Point_1 - Point_2); 
  38.       Vector_2   : constant Double_Vector_3D := Normalized (Point_3 - Point_2); 
  39.       Cos_Theta  : constant Double      := Vector_1 * Vector_2; 
  40.  
  41.    begin 
  42.       case Cos_Theta >= 1.0 is 
  43.          when True  => return Pi; 
  44.          when False => return Arccos (Cos_Theta); 
  45.       end case; 
  46.    end Angle; 
  47.  
  48.    function to_Degrees (Radians  : Double) return Double is (Radians * 180.0 / Pi); 
  49.  
  50.    function to_Radians (Degrees  : Double) return Double is (Degrees * Pi / 180.0); 
  51.  
  52.    -------------- 
  53.    -- Matrices -- 
  54.    -------------- 
  55.  
  56.    function "*" (A, B : Matrix_33) return Matrix_33 is 
  57.  
  58.       AB : Matrix_33; 
  59.  
  60.    begin 
  61.       for i in Matrix_33'Range (1) loop 
  62.          for j in Matrix_33'Range (2) loop 
  63.             declare 
  64.                r : Double := 0.0; 
  65.             begin 
  66.                for k in Matrix_33'Range (1) loop -- same as Matrix_33'Range (2) 
  67.                   r := r + (A (i, k) * B (k, j)); 
  68.                end loop; 
  69.                AB (i, j) := r; 
  70.             end; 
  71.          end loop; 
  72.       end loop; 
  73.       return AB; 
  74.    end "*"; 
  75.  
  76.    function "*" (A : Matrix_33; x : Double_Vector_3D) return Double_Vector_3D is 
  77.  
  78.       Ax : Double_Vector_3D; 
  79.  
  80.    begin 
  81.       for i in Matrix_33'Range (1) loop 
  82.          declare 
  83.             r : Double := 0.0; 
  84.          begin 
  85.             for j in Matrix_33'Range (2) loop 
  86.                r := r + A (i, j) * x (j - Matrix_33'First (2) + Double_Vector_3D'First); 
  87.             end loop; 
  88.             Ax (i - Matrix_33'First (1) + Double_Vector_3D'First) := r; 
  89.          end; 
  90.       end loop; 
  91.       return Ax; 
  92.    end "*"; 
  93.  
  94.    function "*" (A : Matrix_44; x : Double_Vector_3D) return Double_Vector_3D is 
  95.  
  96.       type Vector_4D_GL_Double is array (0 .. 3) of GL.Double; 
  97.       x_4D : constant Vector_4D_GL_Double := (x (0), x (1), x (2), 1.0); 
  98.       Ax   :          Vector_4D_GL_Double; 
  99.  
  100.       -- banana skin : Matrix has range 1 .. 3, Vector 0 .. 2 (GL) 
  101.    begin 
  102.       for i in Matrix_44'Range (1) loop 
  103.          declare 
  104.             r : Double := 0.0; 
  105.          begin 
  106.             for j in Matrix_44'Range (2) loop 
  107.                r := r + A (i, j) * x_4D (j - Matrix_44'First (2) + Vector_4D_GL_Double'First); 
  108.             end loop; 
  109.             Ax (i - Matrix_44'First (1) + Vector_4D_GL_Double'First) := r; 
  110.          end; 
  111.       end loop; 
  112.       return (Ax (0), Ax (1), Ax (2)); 
  113.    end "*"; 
  114.  
  115.    function "*" (A : Matrix_44; x : Double_Vector_3D) return Vector_4D is 
  116.  
  117.       x_4D  : constant Vector_4D := (x (0), x (1), x (2), 1.0); 
  118.       Ax    : Vector_4D; 
  119.  
  120.    begin 
  121.       for i in Matrix_44'Range (1) loop 
  122.          declare 
  123.             r : Double := 0.0; 
  124.          begin 
  125.             for j in Matrix_44'Range (2) loop 
  126.                r := r + A (i, j) * x_4D (j - Matrix_44'First (2) + Vector_4D'First); 
  127.             end loop; 
  128.             Ax (i - Matrix_44'First (2) + Vector_4D'First) := r; 
  129.          end; 
  130.       end loop; 
  131.       return Ax; 
  132.    end "*"; 
  133.  
  134.    function Transpose (A : Matrix_33) return Matrix_33 is ((A (1, 1), A (2, 1), A (3, 1)), 
  135.                                                            (A (1, 2), A (2, 2), A (3, 2)), 
  136.                                                            (A (1, 3), A (2, 3), A (3, 3))); 
  137.  
  138.    function Transpose (A : Matrix_44) return Matrix_44 is ((A (1, 1), A (2, 1), A (3, 1), A (4, 1)), 
  139.                                                            (A (1, 2), A (2, 2), A (3, 2), A (4, 2)), 
  140.                                                            (A (1, 3), A (2, 3), A (3, 3), A (4, 3)), 
  141.                                                            (A (1, 4), A (2, 4), A (3, 4), A (4, 4))); 
  142.  
  143.    function Det (A : Matrix_33) return Double is   (A (1, 1) * A (2, 2) * A (3, 3) 
  144.                                                     + A (2, 1) * A (3, 2) * A (1, 3) 
  145.                                                     + A (3, 1) * A (1, 2) * A (2, 3) 
  146.                                                     - A (3, 1) * A (2, 2) * A (1, 3) 
  147.                                                     - A (2, 1) * A (1, 2) * A (3, 3) 
  148.                                                     - A (1, 1) * A (3, 2) * A (2, 3)); 
  149.  
  150.    function XYZ_rotation (ax, ay, az : Double) return Matrix_33 is 
  151.  
  152.       Mx, My, Mz : Matrix_33; c, s : Double; 
  153.  
  154.    begin 
  155.       -- Around X 
  156.       c := Cos (ax); 
  157.       s := Sin (ax); 
  158.       Mx := ((1.0, 0.0, 0.0), 
  159.              (0.0,   c,  -s), 
  160.              (0.0,   s,   c)); 
  161.       -- Around Y 
  162.       c := Cos (ay); 
  163.       s := Sin (ay); 
  164.       My := ((c,   0.0,  -s), 
  165.              (0.0, 1.0, 0.0), 
  166.              (s,   0.0,   c)); 
  167.       -- Around Z 
  168.       c := Cos (az); 
  169.       s := Sin (az); 
  170.       Mz := ((c,    -s, 0.0), 
  171.              (s,     c, 0.0), 
  172.              (0.0, 0.0, 1.0)); 
  173.  
  174.       return Mz * My * Mx; 
  175.  
  176.    end XYZ_rotation; 
  177.  
  178.    function XYZ_rotation (v : Double_Vector_3D) return Matrix_33 is (XYZ_rotation (v (0), v (1), v (2))); 
  179.  
  180.    function Look_at (direction : Double_Vector_3D) return Matrix_33 is 
  181.  
  182.       v1, v2, v3 : Double_Vector_3D; 
  183.  
  184.    begin 
  185.       -- GL's look direction is the 3rd dimension (z) 
  186.       v3 := Normalized (direction); 
  187.       v2 := Normalized ((v3 (2), 0.0, -v3 (0))); 
  188.       v1 := v2 * v3; 
  189.       return ((v1 (0), v2 (0), v3 (0)), 
  190.               (v1 (1), v2 (1), v3 (1)), 
  191.               (v1 (2), v2 (2), v3 (2))); 
  192.    end Look_at; 
  193.  
  194.    function Sub_Matrix (Self                : Matrix; 
  195.                         start_Row, end_Row, 
  196.                         start_Col, end_Col  : Positive) return Matrix is 
  197.  
  198.       the_sub_Matrix  : Matrix (1 .. end_Row - start_Row + 1, 
  199.                                 1 .. end_Col - start_Col + 1); 
  200.  
  201.    begin 
  202.       for Row in the_sub_Matrix'Range (1) loop 
  203.          for Col in the_sub_Matrix'Range (2) loop 
  204.             the_sub_Matrix (Row, Col) := Self (Row + start_Row - 1, 
  205.                                                Col + start_Col - 1); 
  206.          end loop; 
  207.       end loop; 
  208.  
  209.       return the_sub_Matrix; 
  210.    end Sub_Matrix; 
  211.  
  212.    function Look_at (eye, center, up  : Double_Vector_3D) return Matrix_33 is 
  213.  
  214.       forward  : constant Double_Vector_3D := Normalized ((center (0) - eye (0),  center (1) - eye (1),  center (2) - eye (2))); 
  215.       side     : constant Double_Vector_3D := Normalized (forward * up); 
  216.       new_up   : constant Double_Vector_3D := side * forward; 
  217.  
  218.    begin 
  219.       return ((side     (0),  side    (1),   side    (2)), 
  220.               (new_up   (0),  new_up  (1),   new_up  (2)), 
  221.               (-forward (0), -forward (1), -forward (2))); 
  222.    end Look_at; 
  223.  
  224.    -- Following procedure is from Project Spandex, by Paul Nettle 
  225.    procedure Re_Orthonormalize (M : in out Matrix_33) is 
  226.  
  227.       dot1, dot2, vlen : Double; 
  228.  
  229.    begin 
  230.       dot1 := M (1, 1) * M (2, 1) + M (1, 2) * M (2, 2) + M (1, 3) * M (2, 3); 
  231.       dot2 := M (1, 1) * M (3, 1) + M (1, 2) * M (3, 2) + M (1, 3) * M (3, 3); 
  232.  
  233.       M (1, 1) := M (1, 1) - dot1 * M (2, 1) - dot2 * M (3, 1); 
  234.       M (1, 2) := M (1, 2) - dot1 * M (2, 2) - dot2 * M (3, 2); 
  235.       M (1, 3) := M (1, 3) - dot1 * M (2, 3) - dot2 * M (3, 3); 
  236.  
  237.       vlen := 1.0 / Sqrt (M (1, 1) * M (1, 1) + 
  238.                             M (1, 2) * M (1, 2) + 
  239.                             M (1, 3) * M (1, 3)); 
  240.  
  241.       M (1, 1) := M (1, 1) * vlen; 
  242.       M (1, 2) := M (1, 2) * vlen; 
  243.       M (1, 3) := M (1, 3) * vlen; 
  244.  
  245.       dot1 := M (2, 1) * M (1, 1) + M (2, 2) * M (1, 2) + M (2, 3) * M (1, 3); 
  246.       dot2 := M (2, 1) * M (3, 1) + M (2, 2) * M (3, 2) + M (2, 3) * M (3, 3); 
  247.  
  248.       M (2, 1) := M (2, 1) - dot1 * M (1, 1) - dot2 * M (3, 1); 
  249.       M (2, 2) := M (2, 2) - dot1 * M (1, 2) - dot2 * M (3, 2); 
  250.       M (2, 3) := M (2, 3) - dot1 * M (1, 3) - dot2 * M (3, 3); 
  251.  
  252.       vlen := 1.0 / Sqrt (M (2, 1) * M (2, 1) + 
  253.                             M (2, 2) * M (2, 2) + 
  254.                             M (2, 3) * M (2, 3)); 
  255.  
  256.       M (2, 1) := M (2, 1) * vlen; 
  257.       M (2, 2) := M (2, 2) * vlen; 
  258.       M (2, 3) := M (2, 3) * vlen; 
  259.  
  260.       M (3, 1) := M (1, 2) * M (2, 3) - M (1, 3) * M (2, 2); 
  261.       M (3, 2) := M (1, 3) * M (2, 1) - M (1, 1) * M (2, 3); 
  262.       M (3, 3) := M (1, 1) * M (2, 2) - M (1, 2) * M (2, 1); 
  263.    end Re_Orthonormalize; 
  264.  
  265.    --    type Matrix_44 is array (0 .. 3, 0 .. 3) of aliased Double; -- for GL.MultMatrix 
  266.    --    pragma Convention (Fortran, Matrix_44);  -- GL stores matrices columnwise 
  267.  
  268.    M : Matrix_44; 
  269.    -- M is a global variable for a clean 'Access and for setting once 4th dim 
  270.    pM : constant GL.doublePtr := M (1, 1)'Unchecked_Access; 
  271.  
  272.    procedure Multiply_GL_Matrix (A : Matrix_33) is 
  273.  
  274.    begin 
  275.       for i in Matrix_33'Range (1) loop 
  276.          for j in Matrix_33'Range (2) loop 
  277.             M (i, j) := A (i, j); 
  278.             --  Funny deformations .. . 
  279.             --  if j=2 then 
  280.             --    M (i - 1, j - 1) := 0.5 * A (i, j); 
  281.             --  end if; 
  282.          end loop; 
  283.       end loop; 
  284.       GL.MultMatrixd (pM); 
  285.    end Multiply_GL_Matrix; 
  286.  
  287.    procedure Set_GL_Matrix (A : Matrix_33) is 
  288.  
  289.    begin 
  290.       GL.LoadIdentity; 
  291.       Multiply_GL_Matrix (A); 
  292.    end Set_GL_Matrix; 
  293.  
  294.    -- Ada 95 Quality and Style Guide, 7.2.7: 
  295.    -- Tests for 
  296.    -- 
  297.    -- (1) absolute "equality" to 0 in storage, 
  298.    -- (2) absolute "equality" to 0 in computation, 
  299.    -- (3) relative "equality" to 0 in storage, and 
  300.    -- (4) relative "equality" to 0 in computation: 
  301.    -- 
  302.    --  abs X <= Float_Type'Model_Small                      -- (1) 
  303.    --  abs X <= Float_Type'Base'Model_Small                 -- (2) 
  304.    --  abs X <= abs X * Float_Type'Model_Epsilon            -- (3) 
  305.    --  abs X <= abs X * Float_Type'Base'Model_Epsilon       -- (4) 
  306.  
  307.    function Almost_zero (x : Double) return Boolean is (abs x <= Double'Base'Model_Small); 
  308.  
  309.    function Almost_zero (x : GL.C_Float) return Boolean is (abs x <= GL.C_Float'Base'Model_Small); 
  310.  
  311. begin 
  312.    for i in 1 .. 3 loop 
  313.       M (i, 4) := 0.0; 
  314.       M (4, i) := 0.0; 
  315.    end loop; 
  316.    M (4, 4) := 1.0; 
  317. end GL.Math;