#ifndef matrix_h
#define matrix_h

#include <iostream.h>
#include "mathex.h"

class TVector;

class TMatrix33
{
	private:
		double _Mx[3][3];

		// Input output
		ostream &write(ostream &out) const;
		istream &read(istream &in);

	public:
		// Constructors
		TMatrix33();
		TMatrix33(double Phi, double Theta, double Psi);
		TMatrix33(double mx00, double mx01, double mx02, double mx10, double mx11, double mx12, double mx20, double mx21, double mx22);

		// Selectors
		double operator()(int Row, int Column) const { return _Mx[Row][Column]; }
		double &operator()(int Row, int Column) { return _Mx[Row][Column]; }

		// Optimised artimetric methods
		static TMatrix33 &add(const TMatrix33 &m1, const TMatrix33 &m2, TMatrix33 &result);
		static TMatrix33 &subtract(const TMatrix33 &m1, const TMatrix33 &m2, TMatrix33 &result);
		static TMatrix33 &multiply(const TMatrix33 &m1, const TMatrix33 &m2, TMatrix33 &result);
		static TMatrix33 &multiply(const TMatrix33 &m1, const double &scale, TMatrix33 &result);
		static TVector &multiply(const TMatrix33 &m1, const TVector &v, TVector &result);

		// Matrix arithmetic
		TMatrix33 &operator+=(const TMatrix33 &m) { return add(*this, m, *this); }
		TMatrix33 &operator-=(const TMatrix33 &m) { return subtract(*this, m, *this); }
		TMatrix33 &operator*=(const TMatrix33 &m) { TMatrix33 tm(*this); return multiply(tm, m, *this); }
		TMatrix33 &operator*=(const double &scale) { return multiply(*this, scale, *this); }
		TMatrix33 operator+(const TMatrix33 &m) const { TMatrix33 tm; return add(*this, m, tm); }
		TMatrix33 operator-(const TMatrix33 &m) const { TMatrix33 tm; return subtract(*this, m, tm); }
		TMatrix33 operator*(const TMatrix33 &m) const { TMatrix33 tm; return multiply(*this, m, tm); }
		TMatrix33 operator*(const double &scale) const { TMatrix33 tm; return multiply(*this, scale, tm); }
		TVector operator*(const TVector &v) const;

		// Determinants
		double determinant() const;
		static double determinant(const TMatrix33 &m) { return m.determinant(); }

		// Transpose
		TMatrix33 &transpose();
		static TMatrix33 &transpose(const TMatrix33 &m, TMatrix33 &result) { result = m; return result.transpose(); }
		static TMatrix33 transpose(const TMatrix33 &m) { return TMatrix33(m).transpose(); }

		// Inverse
		static TMatrix33 &TMatrix33::inverse(const TMatrix33 &m1, TMatrix33 &result);
		static TMatrix33 TMatrix33::inverse(const TMatrix33 &m1) { TMatrix33 tm; return inverse(m1, tm); }
		TMatrix33 inverse() const { TMatrix33 tm; return inverse(*this, tm); }

		// Stream
		friend ostream &operator<<(ostream &out, const TMatrix33 &o) { return o.write(out); }
		friend istream &operator>>(istream &in, TMatrix33 &o) { return o.read(in); }
};
#endif