#ifndef __GELMATRIX_H
#define __GELMATRIX_H

#include "ge_lvector.h"
#include "ge_matrix.h"

struct LMatrix
{
	LVector off,v1,v2,v3;

	LMatrix(void)
	{
		off=LVector(0.0,0.0,0.0); v1=LVector(1.0,0.0,0.0);
		v2 =LVector(0.0,1.0,0.0); v3=LVector(0.0,0.0,1.0);
	}
	LMatrix(const LVector &off_in, const LVector &v1_in, const LVector &v2_in, const LVector &v3_in)
	{
		off=off_in;v1=v1_in;v2=v2_in;v3=v3_in;
	}
	LMatrix(_DONTCONSTRUCT v): off(v), v1(v), v2(v), v3(v){}

	friend const LMatrix operator * (const LMatrix &m1, const LMatrix &m2)
	{
		LMatrix erg(DC);

		erg.off.x = m1.off.x + m1.v1.x*m2.off.x + m1.v2.x*m2.off.y + m1.v3.x*m2.off.z;
		erg.off.y = m1.off.y + m1.v1.y*m2.off.x + m1.v2.y*m2.off.y + m1.v3.y*m2.off.z;
		erg.off.z = m1.off.z + m1.v1.z*m2.off.x + m1.v2.z*m2.off.y + m1.v3.z*m2.off.z;

		erg.v1.x = m1.v1.x*m2.v1.x + m1.v2.x*m2.v1.y + m1.v3.x*m2.v1.z;
		erg.v1.y = m1.v1.y*m2.v1.x + m1.v2.y*m2.v1.y + m1.v3.y*m2.v1.z;
		erg.v1.z = m1.v1.z*m2.v1.x + m1.v2.z*m2.v1.y + m1.v3.z*m2.v1.z;

		erg.v2.x = m1.v1.x*m2.v2.x + m1.v2.x*m2.v2.y + m1.v3.x*m2.v2.z;
		erg.v2.y = m1.v1.y*m2.v2.x + m1.v2.y*m2.v2.y + m1.v3.y*m2.v2.z;
		erg.v2.z = m1.v1.z*m2.v2.x + m1.v2.z*m2.v2.y + m1.v3.z*m2.v2.z;

		erg.v3.x = m1.v1.x*m2.v3.x + m1.v2.x*m2.v3.y + m1.v3.x*m2.v3.z;
		erg.v3.y = m1.v1.y*m2.v3.x + m1.v2.y*m2.v3.y + m1.v3.y*m2.v3.z;
		erg.v3.z = m1.v1.z*m2.v3.x + m1.v2.z*m2.v3.y + m1.v3.z*m2.v3.z;
		return erg;
	}

	friend const LVector operator * (const LMatrix &m, const LVector &v)
	{
		return LVector(m.off.x+m.v1.x*v.x+m.v2.x*v.y+m.v3.x*v.z,
									 m.off.y+m.v1.y*v.x+m.v2.y*v.y+m.v3.y*v.z,
									 m.off.z+m.v1.z*v.x+m.v2.z*v.y+m.v3.z*v.z);
	}

		friend const LVector operator ^ (const LVector &v, const LMatrix &m)
	{
		return LVector(m.v1.x*v.x+m.v2.x*v.y+m.v3.x*v.z,
									 m.v1.y*v.x+m.v2.y*v.y+m.v3.y*v.z,
									 m.v1.z*v.x+m.v2.z*v.y+m.v3.z*v.z);
	}

	friend const LVector operator * (const LVector &v, const LMatrix &m)
	{
		return LVector(m.off.x+m.v1.x*v.x+m.v2.x*v.y+m.v3.x*v.z,
									 m.off.y+m.v1.y*v.x+m.v2.y*v.y+m.v3.y*v.z,
									 m.off.z+m.v1.z*v.x+m.v2.z*v.y+m.v3.z*v.z);
	}

	friend const LMatrix operator * (const LReal s, const LMatrix &m)
	{
		return LMatrix(s*m.off,s*m.v1,s*m.v2,s*m.v3);
	}

	friend const LMatrix operator / (const LMatrix &m, const LReal s)
	{
		return LMatrix(m.off/s,m.v1/s,m.v2/s,m.v3/s);
	}

	friend const LMatrix operator + (const LMatrix &m1, const LMatrix &m2)
	{
		return LMatrix(m1.off+m2.off,m1.v1+m2.v1,m1.v2+m2.v2,m1.v3+m2.v3);
	}

	friend const LMatrix operator - (const LMatrix &m1, const LMatrix &m2)
	{
		return LMatrix(m1.off-m2.off,m1.v1-m2.v1,m1.v2-m2.v2,m1.v3-m2.v3);
	}

	friend const LMatrix operator ! (const LMatrix &m)
	{
		LMatrix mi(DC);

    LReal det = (m.v1.x*(m.v2.y*m.v3.z - m.v3.y*m.v2.z) +
                 m.v2.x*(m.v3.y*m.v1.z - m.v1.y*m.v3.z) +
                 m.v3.x*(m.v1.y*m.v2.z - m.v2.y*m.v1.z));
		if (det==0.0) return LMatrix();

		det = 1.0/det;
    mi.off.x = (m.v2.x *(m.off.y*m.v3.z  - m.v3.y *m.off.z)+
                m.v3.x *(m.off.z*m.v2.y  - m.off.y*m.v2.z)+
                m.off.x*(m.v3.y *m.v2.z  - m.v2.y *m.v3.z))*det;
    mi.off.y = (m.v3.x *(m.off.y*m.v1.z  - m.v1.y *m.off.z)+
                m.off.x*(m.v1.y *m.v3.z  - m.v3.y *m.v1.z) +
                m.v1.x *(m.v3.y *m.off.z - m.off.y*m.v3.z))*det;
    mi.off.z = (m.off.x*(m.v2.y *m.v1.z  - m.v1.y *m.v2.z) +
                m.v1.x *(m.v2.z *m.off.y - m.v2.y *m.off.z)+
                m.v2.x *(m.off.z*m.v1.y  - m.off.y*m.v1.z))*det;

    mi.v1.x = (m.v2.y*m.v3.z - m.v3.y*m.v2.z)*det;
    mi.v1.y = (m.v3.y*m.v1.z - m.v1.y*m.v3.z)*det;
    mi.v1.z = (m.v1.y*m.v2.z - m.v2.y*m.v1.z)*det;

    mi.v2.x = (m.v2.z*m.v3.x - m.v3.z*m.v2.x)*det;
    mi.v2.y = (m.v3.z*m.v1.x - m.v1.z*m.v3.x)*det;
    mi.v2.z = (m.v1.z*m.v2.x - m.v2.z*m.v1.x)*det;

    mi.v3.x = (m.v2.x*m.v3.y - m.v3.x*m.v2.y)*det;
    mi.v3.y = (m.v3.x*m.v1.y - m.v1.x*m.v3.y)*det;
    mi.v3.z = (m.v1.x*m.v2.y - m.v2.x*m.v1.y)*det;

		return mi;
	}
};// LMatrix

inline const LMatrix LM(const Matrix &m)  {return LMatrix(LV(m.off),LV(m.v1),LV(m.v2),LV(m.v3));}
inline const Matrix  SM(const LMatrix &m) {return Matrix (SV(m.off),SV(m.v1),SV(m.v2),SV(m.v3));}

#endif
