indexing
	description: " A 4x4 Matrix with basic operations"
	date: "$Date$"
	revision: "$Revision$"

class
	EM_MATRIX44_REF

inherit
	DOUBLE_MATH
		export
			{NONE} all
		redefine
			out,
			default_create,
			copy,
			is_equal
		end
	ANY
		redefine
			out,
			default_create,
			copy,
			is_equal
		end

create
	default_create,
	make_empty,
	make_unit,
	--make_from_tuple,
	make,
	make_from_quaternion,
	make_from_scalar,
	make_from_translation,
	make_from_rotation,
	make_from_reference

feature -- Access
	set_element( value: like element ; c: INTEGER ; r: INTEGER) is
			-- Sets the element in the c'th column and the r'th row
		require
			1<=c and c<=4
			1<=r and r<=4
		do
			area.put( value, (r-1)*4 + (c-1) )
		end

	set_diagonal (v: EM_VECTOR4D) is
			-- Set diagonal to `v'
		do
			set_element (v.x, 1, 1)
			set_element (v.y, 2, 2)
			set_element (v.z, 3, 3)
			set_element (v.w, 4, 4)
		ensure
			set: diagonal = v
		end


	set_row (v: EM_VECTOR4D; i: INTEGER) is
			-- Set `i'th row to `v'
		require
			valid_index: 1 <= i and i <= 4
		do
			set_element (v.x, 1, i)
			set_element (v.y, 2, i)
			set_element (v.z, 3, i)
			set_element (v.w, 4, i)
		ensure
			set: row (i) = v
		end

	set_column (v: EM_VECTOR4D; i: INTEGER) is
			-- Set `i'th column to `v'
		require
			valid_index: 1 <= i and i <= 4
		do
			set_element (v.x, i, 1)
			set_element (v.y, i, 2)
			set_element (v.z, i, 3)
			set_element (v.w, i, 4)
		ensure
			set: column (i) = v
		end

	element alias "[]" ( c: INTEGER; r: INTEGER): DOUBLE assign set_element is
			-- returns the element in the c'th column and the r'th row
		require
			1<=c and c<=4
			1<=r and r<=4
		do
			result := area.item( 4*(r-1) + (c-1) )
		end

	diagonal: EM_VECTOR4D is
			-- Diagonal of current matrix as vector
		do
			Result.set (
				element (1, 1),
				element (2, 2),
				element (3, 3),
				element (4, 4)
			)
		end


	row (i: INTEGER): EM_VECTOR4D is
			-- `i'th row of current matrix as vector
		require
			valid_index: 1 <= i and i <= 4
		do
			Result.set (
				element (1, i),
				element (2, i),
				element (3, i),
				element (4, i)
			)
		end

	column (i: INTEGER): EM_VECTOR4D is
			-- `i'th column of current matrix as vector
		require
			valid_index: 1 <= i and i <= 4
		do
			Result.set (
				element (i, 1),
				element (i, 2),
				element (i, 3),
				element (i, 4)
			)
		end


feature {EM_MATRIX44_REF} -- Access
	area: SPECIAL[DOUBLE]
feature {NONE} -- Initialization
	make_from_reference (v: EM_MATRIX44_REF ) is
			-- Initialize `Current' with `v.area'.
		require
			v_not_void: v /= Void
		do
			make_empty
			area.copy_data( v.area, 0, 0, 16 )
		ensure
			item_set: area.is_equal ( v.area )
		end

feature -- Conversion

	to_reference: EM_MATRIX44_REF is
			-- Associated reference of Current
		do
			create Result
			Result.area.copy_data ( area, 0, 0, 16 )
		ensure
			to_reference_not_void: Result /= Void
		end


feature -- Initialize
	default_create is
			-- Make an ampty Matrix
		do
			make_unit
		end

	make_empty is
			-- Makes an empty Matrix
		do
			create area.make(16)
		ensure
			area.count = 16
		end

	make_from_vector_tuple( a: TUPLE[EM_VECTOR4D, EM_VECTOR4D, EM_VECTOR4D, EM_VECTOR4D] ) is
			-- Creates the matrix from an array (4x4 would be optimal) of column vectors
		local
			vector_ref: EM_VECTOR4D_REF
		do
			make_empty
			vector_ref ?= a.reference_item (1)
			check
				vector_ref/=void
			end

			set_element( vector_ref.x, 1,1)
			set_element( vector_ref.y, 2,1)
			set_element( vector_ref.z, 3,1)
			set_element( vector_ref.w, 4,1)

			vector_ref ?= a.reference_item (2)
			check
				vector_ref/=void
			end

			set_element( vector_ref.x, 1,2)
			set_element( vector_ref.y, 2,2)
			set_element( vector_ref.z, 3,2)
			set_element( vector_ref.w, 4,2)

			vector_ref ?= a.reference_item (3)
			check
				vector_ref/=void
			end

			set_element( vector_ref.x, 1,3)
			set_element( vector_ref.y, 2,3)
			set_element( vector_ref.z, 3,3)
			set_element( vector_ref.w, 4,3)

			vector_ref ?= a.reference_item (4)
			check
				vector_ref/=void
			end

			set_element( vector_ref.x, 1,4)
			set_element( vector_ref.y, 2,4)
			set_element( vector_ref.z, 3,4)
			set_element( vector_ref.w, 4,4)
		end

	make( a_11: DOUBLE; a_12: DOUBLE; a_13: DOUBLE; a_14: DOUBLE;
	      a_21: DOUBLE; a_22: DOUBLE; a_23: DOUBLE; a_24: DOUBLE;
	      a_31: DOUBLE; a_32: DOUBLE; a_33: DOUBLE; a_34: DOUBLE;
	      a_41: DOUBLE; a_42: DOUBLE; a_43: DOUBLE; a_44: DOUBLE ) is
			-- Creates a Matrix from the given elements
		do
			make_empty
			set_element( a_11, 1, 1 )
			set_element( a_21, 2, 1 )
			set_element( a_31, 3, 1 )
			set_element( a_41, 4, 1 )
			set_element( a_12, 1, 2 )
			set_element( a_22, 2, 2 )
			set_element( a_32, 3, 2 )
			set_element( a_42, 4, 2 )
			set_element( a_13, 1, 3 )
			set_element( a_23, 2, 3 )
			set_element( a_33, 3, 3 )
			set_element( a_43, 4, 3 )
			set_element( a_14, 1, 4 )
			set_element( a_24, 2, 4 )
			set_element( a_34, 3, 4 )
			set_element( a_44, 4, 4 )
		end

	make_unit is
			-- Make a unit matrix
		do
			make_from_scalar( 1 )
		end

	make_from_quaternion( q: EM_QUATERNION ) is
			-- Make a matrix from a Quaternion
		local
			x: DOUBLE
			y: DOUBLE
			z: DOUBLE
			w: DOUBLE
		do
			q.normalize
			x := q.v.x
			y := q.v.y
			z := q.v.z
			w := q.w

			make_empty
			make ( 1 - 2 * (y*y + z*z), 2 * (x*y - w*z)    , 2 * (x*z + w*y)     , 0,
			       2 * (x*y + w*z)    , 1 - 2 * (x*x + z*z), 2 * (y*z - w*x)     , 0,
			       2 * (x*z - y*w)    , 2 * (y*z + w*x)    , 1 - 2 * (x*x + y*y) , 0,
			       0                  , 0                  , 0                   , 1)
		end

	make_from_scalar( s: DOUBLE ) is
			-- Make a scalar matrix
		do
			make_empty
			set_element( s, 1, 1 )
			set_element( s, 2, 2 )
			set_element( s, 3, 3 )
			set_element( s, 4, 4 )
		end

	make_from_translation( v: EM_VECTOR3D ) is
			-- Make a scalar matrix
		do
			make_unit
			set_element( v.x, 1, 4 )
			set_element( v.y, 2, 4 )
			set_element( v.z, 3, 4 )
			set_element( 1, 4, 4 )
		end

	make_from_rotation( v: EM_VECTOR3D; angle: DOUBLE ) is
			-- Make a scalar matrix
		local
			s: DOUBLE
			c: DOUBLE
			vv: EM_VECTOR3D
			x: DOUBLE
			y: DOUBLE
			z: DOUBLE
		do
			s := sine ( angle )
			c := cosine ( angle )
			vv := v.normalized
			x := vv.x
			y := vv.y
			z := vv.z

			make_empty

			set_element( (x*x*(1-c)) + c, 1,1)
			set_element( (y*x*(1-c)) + z*s, 1,2)
			set_element( (z*x*(1-c)) - y*s, 1,3)
			set_element( 0, 1,4)

			set_element( (x*y*(1-c)) - z*s, 2,1)
			set_element( (y*y*(1-c)) + c, 2,2)
			set_element( (z*y*(1-c)) + x*s, 2,3)
			set_element( 0, 2,4)

			set_element( (x*z*(1-c)) + y*s, 3,1)
			set_element( (y*z*(1-c)) - x*s, 3,2)
			set_element( (z*z*(1-c)) + c, 3,3)
			set_element( 0, 3,4)

			set_element( 0, 4,1)
			set_element( 0, 4,2)
			set_element( 0, 4,3)
			set_element( 1, 4,4)
		end

feature -- Conversion
	to_c: ANY is
			-- Address of actual sequence of values,
			-- for passing to external (non-Eiffel) routines.
		do
			result := area
		end

	to_pointer: POINTER is
			-- Address of actual sequence of values,
			-- for passing to external (non-Eiffel) routines.
		do
			result := $area
		end

	to_c_real: ANY is
			-- Address of actual sequence of values,
			-- for passing to external (non-Eiffel) routines.
		obsolete
			"Do not use this function its dangerous, use a real matrix instead"
		local
			tmp: ARRAY[ REAL ]
			i: INTEGER
		do
			create tmp.make( 0, 15 )
			from i:=0 until i>15 loop
				tmp[i] := area[i]
			end
			result := tmp.to_c
		end



feature -- Commands
	multiply (other: like current) is
			-- Multiply two matrices
		do
			copy( current * other )
		end

	infix "*" (other: like current): like Current is
			-- Multiply two matrices
		local
			i: INTEGER
			j: INTEGER
			k: INTEGER
			t: DOUBLE
		do
		-- ISE BUG: I need to create this expanded class or I'll get a "call on void target"
		-- Its quite strange... :(
			create result
			from i := 1 until i > 4 loop -- i'th row
				from j := 1 until j > 4 loop -- j'th column
					t := 0
					from k := 1 until k > 4 loop
						t := t + element(i, k)*other.element (k, j)
						k := k + 1
					end
					result.set_element( t, i, j )
					j := j + 1
				end
				i := i + 1
			end
		end

	infix "-" (other: like current): like Current is
			-- Multiply two matrices
		local
			i: INTEGER
			j: INTEGER
		do
			create result
			from i := 1 until i > 4 loop -- i'th row
				from j := 1 until j > 4 loop -- j'th column
					result[i,j] := element(i,j)-other[i,j]
					j := j + 1
				end
				i := i + 1
			end
		end

	infix "+" (other: like current): like Current is
			-- Multiply two matrices
		local
			i: INTEGER
			j: INTEGER
		do
			create result
			from i := 1 until i > 4 loop -- i'th row
				from j := 1 until j > 4 loop -- j'th column
					result[i,j] := element(i,j) + other[i,j]
					j := j + 1
				end
				i := i + 1
			end
		end

	mult (other: EM_VECTOR4D): EM_VECTOR4D is
			-- Multiply a matrix and a vector
		do
			result.set( element(1,1)*other.x + element(1,2)*other.y + element(1,3)*other.z + element(1,4)*other.w,
						element(2,1)*other.x + element(2,2)*other.y + element(2,3)*other.z + element(2,4)*other.w,
						element(3,1)*other.x + element(3,2)*other.y + element(3,3)*other.z + element(3,4)*other.w,
						element(4,1)*other.x + element(4,2)*other.y + element(4,3)*other.z + element(4,4)*other.w)
		end

	scale (factor: DOUBLE) is
			-- Scales the 'Current' matrix by a 'factor'.
		local
			i: INTEGER
		do
			from i := area.lower until i > area.upper loop -- i'th element
				area[i] := area[i] * factor
				i := i + 1
			end
		end


	scaled (factor: DOUBLE): like Current is
			-- returns a the a scaled copy of the current matrix
		do
			result := twin
			result.scale( factor )
		end


	swap_rows( a, b: INTEGER) is
			-- Swap two rows
		local
			tmp: DOUBLE
		do
			tmp := element(a,1); set_element( element(b,1), a,1); set_element( tmp, a,1)
			tmp := element(a,2); set_element( element(b,2), a,2); set_element( tmp, a,2)
			tmp := element(a,3); set_element( element(b,3), a,3); set_element( tmp, a,3)
			tmp := element(a,4); set_element( element(b,4), a,4); set_element( tmp, a,4)
		end

	inverse is
			-- return the inverse, using gauss
			-- needs det!=0
		require
			det.abs > 0.00000001 -- tolerance
		local
			tmp: DOUBLE
			a,b: EM_MATRIX44
			i,j,jj: INTEGER
			jmax: INTEGER
		do
			a := current
			b.set_unit
			from i:=1 until i>4 loop
				jmax := i
				from j := i+1 until j>4 loop
					if (a[i,jmax]).abs<(a[i,j]).abs then
						jmax := i
					end
					j := j + 1
				end
				check
					(a[i,jmax]).abs>0.00000001
				end
				-- swap the rows
				a.swap_rows(i, jmax)
				b.swap_rows(i, jmax)

				tmp := 1/a[i,i]
				a[1,i] := a[1,i]*tmp
				a[2,i] := a[2,i]*tmp
				a[3,i] := a[3,i]*tmp
				a[4,i] := a[4,i]*tmp

				b[1,i] := b[1,i]*tmp
				b[2,i] := b[2,i]*tmp
				b[3,i] := b[3,i]*tmp
				b[4,i] := b[4,i]*tmp

				from j:=1 until j>4 loop
					if j/=i then
						tmp := a[i,j]
						from jj:=1 until jj>4 loop
							a[jj,j] := a[jj,j] - a[jj,i]*tmp
							b[jj,j] := b[jj,j] - b[jj,i]*tmp
							jj := jj + 1
						end
					end
					j := j + 1
				end

				i:=i+1
			end

			area := b.area
		ensure
			-- This is too much for ES
--			matrix_inversed: (((old current.deep_twin) * current) - (create {EM_MATRIX44_REF}.make_unit)).column_norm <= 0.001
		end


--	column_norm: DOUBLE is //TODO: brauchts die?
--			-- Calculate the 1-norm of the matrix
--		local
--			d:DOUBLE
--		do
--			result := 0
--			d := element(1,1).abs+element(2,1).abs+element(3,1).abs+element(4,1).abs
--			if d>result then
--				result := d
--			end
--			d := element(1,2).abs+element(2,2).abs+element(3,2).abs+element(4,2).abs
--			if d>result then
--				result := d
--			end
--			d := element(1,3).abs+element(2,3).abs+element(3,3).abs+element(4,3).abs
--			if d>result then
--				result := d
--			end
--			d := element(1,4).abs+element(2,4).abs+element(3,4).abs+element(4,4).abs
--			if d>result then
--				result := d
--			end
--		end


	inversed: like current is
			-- return the inverse
			-- needs det!=0
		require
			det.abs > 0.0000000001 -- tolerance
		do
			result := twin
			result.inverse
		end

	det: DOUBLE is
			-- Calculate the determinant
		do
			result := element(1,1) * det3x3( element(2,2), element(2,3), element(2,4), element(3,2), element(3,3), element(3,4), element(4,2), element(4,3), element(4,4))
					- element(2,1) * det3x3( element(1,2), element(1,3), element(1,4), element(3,2), element(3,3), element(3,4), element(4,2), element(4,3), element(4,4))
					+ element(3,1) * det3x3( element(1,2), element(1,3), element(1,4), element(2,2), element(2,3), element(2,4), element(4,2), element(4,3), element(4,4))
					- element(4,1) * det3x3( element(1,2), element(1,3), element(1,4), element(2,2), element(2,3), element(2,4), element(3,2), element(3,3), element(3,4));
		end

	srt_inversed: like current is
			-- assuming this is a homogenous SRT transformation matrix, it returns the inversed transformation matrix
		do
			result := twin
			result.srt_inverse
		end

	srt_inverse is
			-- inverses the matrix assuming it is a homogenous SRT transformation matrix.
		local
			a: EM_MATRIX44
			position:EM_VECTOR4D
			x_axis, y_axis, z_axis: EM_VECTOR3D
			x_scale, y_scale, z_scale: DOUBLE
		do
			a:=current

			x_axis.set(a.element(1,1), a.element(1,2), a.element(1,3))
			y_axis.set(a.element(2,1), a.element(2,2), a.element(2,3))
			z_axis.set(a.element(3,1), a.element(3,2), a.element(3,1))
			x_scale := x_axis.length
			y_scale := y_axis.length
			z_scale := z_axis.length

			position := -a.column(4)
			position.set_element (1,4)
			a.set_element(0.0,4,1)
			a.set_element(0.0,4,2)
			a.set_element(0.0,4,3)

			a.set_row(a.row(1)*(1/x_scale),1)
			a.set_row(a.row(2)*(1/y_scale),2)
			a.set_row(a.row(3)*(1/z_scale),3)

			a.transpose

			position := a.mult(position)

			a.set_row(a.row(1)*(1/x_scale),1)
			a.set_row(a.row(2)*(1/y_scale),2)
			a.set_row(a.row(3)*(1/z_scale),3)

			a.set_column(position,4)
			a.set_element(1,4,4)

			area := a.area
		end

	transpose is
			-- transposes the matrix
		local
		row1,row2,row3,row4 : EM_VECTOR4D
		do
			row1 := row(1)
			row2 := row(2)
			row3 := row(3)
			row4 := row(4)
			set_column(row1,1)
			set_column(row2,2)
			set_column(row3,3)
			set_column(row4,4)
		end

feature {NONE} --Commands
	det2x2( a: DOUBLE; b: DOUBLE; c:DOUBLE; d: DOUBLE ):DOUBLE is
			-- Calculate the determinant of a 2x2 matrix
		do
			result := a*d - b*c
		end

	det3x3( a1: DOUBLE; a2: DOUBLE; a3: DOUBLE; b1: DOUBLE; b2: DOUBLE; b3: DOUBLE; c1:DOUBLE; c2:DOUBLE; c3:DOUBLE ):DOUBLE is
			-- Calculate the determinant of a 2x2 matrix
		do
			result := a1 * det2x2( b2, b3, c2, c3 )
        			- b1 * det2x2( a2, a3, c2, c3 )
        			+ c1 * det2x2( a2, a3, b2, b3 )
		end

feature -- Output

	out: STRING is
			-- Output the matrix to a string
		do
			Result := "! " + element(1,1).out + " 	 " + element(1,2).out + " 	 " + element(1,3).out + " 	 " + element(1,4).out + " !%N" +
					  "! " + element(2,1).out + " 	 " + element(2,2).out + " 	 " + element(2,3).out + " 	 " + element(2,4).out + " !%N" +
					  "! " + element(3,1).out + " 	 " + element(3,2).out + " 	 " + element(3,3).out + " 	 " + element(3,4).out + " !%N" +
					  "! " + element(4,1).out + " 	 " + element(4,2).out + " 	 " + element(4,3).out + " 	 " + element(4,4).out + " !%N"
		end

feature -- from any
	copy (other: like Current) is
			-- Reinitialize by copying all the items of `other'.
			-- (This is also used by `clone'.)
		do
			area := other.area.twin
		ensure then
			equal_areas: area.is_equal (other.area)
		end
	is_equal (other: like Current): BOOLEAN is
			-- Is array made of the same items as `other'?
		local
			i: INTEGER
		do
			Result := area.same_items (other.area, area.upper )
		end
end