function matrix_ctor (M)

   M:identity ()
end

matrix = class (matrix_ctor)

function matrix:identity ()

   self.data = { { 1, 0, 0, 0, }, { 0, 1, 0, 0, }, { 0, 0, 1, 0, }, { 0, 0, 0, 1, }, }
end

function matrix:vector_product (V)

   local VNew = {}
   VNew[1] = self.data[1][1] * V[1] + self.data[2][1] * V[2] + self.data[3][1] * V[3] + self.data [4][1]
   VNew[2] = self.data[1][2] * V[1] + self.data[2][2] * V[2] + self.data[3][2] * V[3] + self.data [4][2]
   VNew[3] = self.data[1][3] * V[1] + self.data[2][3] * V[2] + self.data[3][3] * V[3] + self.data [4][3]
   return VNew
end

function matrix:vector_apply (V)

   local VNew = self:vector_product (V)
   V[1] = VNew[1]
   V[2] = VNew[2]
   V[3] = VNew[3]
end

function matrix:matrix_product (M)

   local MNew = matrix ()
   for X = 1, 4, 1 do
      for Y = 1, 4, 1 do
         MNew.data[X][Y] = self.data[1][Y] * M.data[X][1] +
                           self.data[2][Y] * M.data[X][2] +
                           self.data[3][Y] * M.data[X][3] +
                           self.data[4][Y] * M.data[X][4]
      end
   end
   return MNew
end

function matrix:matrix_multiply (M)

   local MNew = self:matrix_product (M)
   self.data = MNew.data
end

function matrix:init_scaling (X, Y, Z)

   if (not Y) then Y = X end
   if (not Z) then Z = X end
   self:identity ()
   self.data[1][1] = X
   self.data[2][2] = Y
   self.data[3][3] = Z
end

function matrix:scale (X, Y, Z)

   local MScale = matrix ()
   MScale:init_scaling (X, Y, Z)
   self:matrix_multiply (MScale)
end

function matrix:init_rotation (X, Y, Z, Angle)

   self:identity ()

   local L2 = X * X + Y * Y + Z * Z
   if (0 == Angle or 0 == L2) then return
   end

   local L = math.sqrt (L2)
   X = X / L
   Y = Y / L
   Z = Z / L

   local c = math.cos (Angle)
   local s = math.sin (Angle)
   local cc = 1 - c

   self.data[1][1] = (cc * X * X) + c;
   self.data[1][2] = (cc * X * Y) + (Z * s);
   self.data[1][3] = (cc * X * Z) - (Y * s);

   self.data[2][1] = (cc * X * Y) - (Z * s);
   self.data[2][2] = (cc * Y * Y) + c;
   self.data[2][3] = (cc * Z * Y) + (X * s);

   self.data[3][1] = (cc * X * Z) + (Y * s);
   self.data[3][2] = (cc * Y * Z) - (X * s);
   self.data[3][3] = (cc * Z * Z) + c;
end

function matrix:rotate (X, Y, Z, Angle)

   local MRotate = matrix ()
   MRotate:init_rotation (X, Y, Z, Angle)
   self:matrix_multiply (MRotate)
end

function matrix:init_translation (X, Y, Z)

   self:identity ()
   self.data[4][1] = X
   self.data[4][2] = Y
   self.data[4][3] = Z
end

function matrix:translate (X, Y, Z)

   local MTranslate = matrix ()
   MTranslate:init_translation (X, Y, Z)
   self:matrix_multiply (MTranslate)
end
