\section{RSAGL.Matrix}

\begin{code}

module RSAGL.Math.Matrix
    (Matrix,
     matrix,
     transformHomogenous,
     rowMajorForm,
     colMajorForm,
     rowAt,
     matrixAt,
     identity_matrix,
     translationMatrix,
     rotationMatrix,
     scaleMatrix,
     xyzMatrix,
     matrixAdd,
     matrixMultiply,
     matrixTranspose,
     matrixInverse,
     determinant,
     matrixInversePrim,
     matrixTransposePrim,
     matrixInverseTransposePrim,
     determinantPrim)
    where

import Data.List as List
import RSAGL.Math.Angle
import RSAGL.Math.Vector
import Data.Vec as Vec
import RSAGL.Types
\end{code}

A 4-by-4 matrix with cached inverse, transpose, inverse transpose, and determinant.

\begin{code}
data Matrix = Matrix { matrix_data :: !(Mat44 RSdouble),
                       matrix_inverse :: Matrix,
                       matrix_transpose :: Matrix,
                       matrix_determinant :: RSdouble }

instance Eq Matrix where
    x == y = matrix_data x == matrix_data y

instance Show Matrix where
    show m = show $ rowMajorForm m

rowMajorForm :: Matrix -> [[RSdouble]]
rowMajorForm = matToLists . matrix_data

colMajorForm :: Matrix -> [[RSdouble]]
colMajorForm = List.transpose . rowMajorForm
\end{code}

rowAt answers the nth row of a matrix.

\begin{code}
rowAt :: Matrix -> Int -> [RSdouble]
rowAt m n = (rowMajorForm m) !! n
\end{code}

matrixAt answers the (i'th,j'th) element of a matrix.

\begin{code}
matrixAt :: Matrix -> (Int,Int) -> RSdouble
matrixAt m (i,j) = rowAt m i !! j
\end{code}

\subsection{Constructing matrices}

matrix constructs a matrix from row major list form.  (Such a list form can be formatted correctly in 
monospaced font and haskell syntax, so that it looks like a matrix as it would be normally written.)

\begin{code}
matrix :: [[RSdouble]] -> Matrix
matrix = uncheckedMatrix . matFromLists

-- | Generate a column matrix of length 4, perform an affine transformation on it, and produce the resulting value.
{-# INLINE transformHomogenous #-}
transformHomogenous :: RSdouble -> RSdouble -> RSdouble -> RSdouble -> (RSdouble -> RSdouble -> RSdouble -> a) -> Matrix -> a
transformHomogenous x y z w f m = f x' y' z'
    where (x':.y':.z':._) = multmv (matrix_data m) (x:.y:.z:.w:.())

uncheckedMatrix :: Mat44 RSdouble -> Matrix
uncheckedMatrix dats = m
    where m_inverse = (matrixInversePrim m) { matrix_inverse = m, matrix_transpose = m_inverse_transpose, matrix_determinant = recip m_det }
          m_transpose = (matrixTransposePrim m) { matrix_inverse = m_inverse_transpose, matrix_transpose = m, matrix_determinant = m_det }
          m_inverse_transpose = (matrixTransposePrim m_inverse) { matrix_inverse = m_transpose, matrix_transpose = m_inverse, matrix_determinant = recip m_det }
          m_det = determinantPrim m
          m = Matrix { matrix_data=dats,
                       matrix_inverse = m_inverse,
                       matrix_transpose = m_transpose,
                       matrix_determinant = m_det }
\end{code}

identityMatrix constructs the n by n identity matrix

\begin{code}
identity_matrix :: Matrix
identity_matrix = Matrix {
    matrix_data = identity,
    matrix_inverse = identity_matrix,
    matrix_transpose = identity_matrix,
    matrix_determinant = 1 }
\end{code}

\begin{code}
translationMatrix :: Vector3D -> Matrix
translationMatrix (Vector3D x y z) = uncheckedMatrix $ (1:.0:.0:.x:.()) :.
                                                       (0:.1:.0:.y:.()) :.
                                                       (0:.0:.1:.z:.()) :.
                                                       (0:.0:.0:.1:.()) :. ()
\end{code}

\begin{code}
rotationMatrix :: Vector3D -> Angle -> Matrix
rotationMatrix vector angle = let s = sine angle
				  c = cosine angle
				  c' = 1 - c
				  (Vector3D x y z) = vectorNormalize vector
				  in uncheckedMatrix $ ((c+c'*x*x)   :. (c'*y*x - s*z) :. (c'*z*x + s*y) :. 0 :. ()) :.
					               ((c'*x*y+s*z) :. (c+c'*y*y)     :. (c'*z*y - s*x) :. 0 :. ()) :.
					               ((c'*x*z-s*y) :. (c'*y*z+s*x)   :. (c+c'*z*z)     :. 0 :. ()) :.
					               (0            :. 0              :. 0              :. 1 :. ()) :. ()
\end{code}

\begin{code}
scaleMatrix :: Vector3D -> Matrix
scaleMatrix (Vector3D x y z) = uncheckedMatrix $ (x :. 0 :. 0 :. 0) :.
				                 (0 :. y :. 0 :. 0) :.
				                 (0 :. 0 :. z :. 0) :.
				                 (0 :. 0 :. 0 :. 1) :. ()
\end{code}

\texttt{xyzMatrix} constructs the matrix in which the x y and z axis are transformed to point in the direction of the specified
vectors.

\begin{code}
xyzMatrix :: Vector3D -> Vector3D -> Vector3D -> Matrix
xyzMatrix (Vector3D x1 y1 z1) (Vector3D x2 y2 z2) (Vector3D x3 y3 z3) =
    uncheckedMatrix $ (x1 :. x2 :. x3 :. 0) :.
                      (y1 :. y2 :. y3 :. 0) :.
                      (z1 :. z2 :. z3 :. 0) :.
                      (0  :. 0  :. 0  :. 1.0) :. ()
\end{code}

\subsection{Matrix arithmetic}

\begin{code}
matrixAdd :: Matrix -> Matrix -> Matrix
matrixAdd m n = uncheckedMatrix $ Vec.zipWith (+) (matrix_data m) (matrix_data n)
\end{code}

\label{howtoUseMatrixMultiplyImpl}

\begin{code}
matrixMultiply :: Matrix -> Matrix -> Matrix
matrixMultiply m n = uncheckedMatrix $ multmm (matrix_data m) (matrix_data n)

matrixTranspose :: Matrix -> Matrix
matrixTranspose = matrix_transpose
\end{code}

\subsection{Inverse and determinant of a matrix}

\begin{code}
matrixInverse :: Matrix -> Matrix
matrixInverse = matrix_inverse
\end{code}

\begin{code}
determinant :: Matrix -> RSdouble
determinant = matrix_determinant
\end{code}

\begin{code}
matrixInverseTransposePrim :: Matrix -> Matrix
matrixInverseTransposePrim = uncheckedMatrix . Vec.transpose . fst . invertAndDet . matrix_data

matrixInversePrim :: Matrix -> Matrix
matrixInversePrim = uncheckedMatrix . fst . invertAndDet . matrix_data

determinantPrim :: Matrix -> RSdouble
determinantPrim = det . matrix_data

matrixTransposePrim :: Matrix -> Matrix
matrixTransposePrim = uncheckedMatrix . Vec.transpose . matrix_data
\end{code}