#ifndef M33_H
#define M33_H

#include <math.h>
#include <stdio.h>
#include "vect3d.hpp"

class m33
{
 public:

  float m[3][3];

  m33() { make_identity(); };

  m33(const m33 &m1)
  {
    m[0][0]=m1.m[0][0]; m[0][1]=m1.m[0][1]; m[0][2]=m1.m[0][2];
    m[1][0]=m1.m[1][0]; m[1][1]=m1.m[1][1]; m[1][2]=m1.m[1][2];
    m[2][0]=m1.m[2][0]; m[2][1]=m1.m[2][1]; m[2][2]=m1.m[2][2];
  };

  m33(float m00, float m10, float m20,
      float m01, float m11, float m21,
      float m02, float m12, float m22)
  {
    m[0][0]=m00; m[0][1]=m01; m[0][2]=m02;
    m[1][0]=m10; m[1][1]=m11; m[1][2]=m12;
    m[2][0]=m20; m[2][1]=m21; m[2][2]=m22;
  };

  void make_identity()
  {
    m[0][0]=1.0; m[0][1]=0.0; m[0][2]=0.0;
    m[1][0]=0.0; m[1][1]=1.0; m[1][2]=0.0;
    m[2][0]=0.0; m[2][1]=0.0; m[2][2]=1.0;
  };

  m33 T() const
  {
    m33 m2(m[0][0],m[0][1],m[0][2],
           m[1][0],m[1][1],m[1][2],
           m[2][0],m[2][1],m[2][2]);
    return m2;
  };

  Vect3d col(int c) const
  {
    Vect3d V(m[0][c],m[1][c],m[2][c]);
    return V;
  };

  Vect3d row(int r) const
  {
    Vect3d V(m[r][0],m[r][1],m[r][2]);
    return V;
  };

  void print()
  {
    printf("%f %f %f\n", m[0][0], m[0][1], m[0][2]);
    printf("%f %f %f\n", m[1][0], m[1][1], m[1][2]);
    printf("%f %f %f\n", m[2][0], m[2][1], m[2][2]);
  };
};

m33 operator+(const m33 &m1, const m33 &m2);
m33 operator-(const m33 &m1, const m33 &m2);
m33 operator*(const m33 &m1, const m33 &m2);
m33 operator*(const m33 &m1, const float &s);
m33 operator*(const float &s, const m33 &m1);
m33 T(m33 M);
m33 Identity();
m33 Rx(float t);
m33 Ry(float t);
m33 Rz(float t);
m33 R(float dx, float dy, float dz, float t);

Vect3d solve_system(const m33 &A, const Vect3d &b);
Vect3d operator*(const m33 &ma, const Vect3d &va);
void eigen(m33 *evec, Vect3d *eval, m33 M);
m33 orthonormalize(m33 m);

#endif
/* M33_H */