#include "machdefs.h"
#include "pexceptions.h"

/*
#include <math.h>
#include <iostream.h>
#include <fstream.h>

#include "tarrinit.h"
#include "array.h"
#include "ctimer.h"
#include "timing.h"
*/



// -------------------------------------------------------
// Classe simple de matrice 3x3 - sans allocation memoire
// -------------------------------------------------------

class Matrix3x3 {
public:
  inline        Matrix3x3(bool raz=true) 
  { 
    if (raz) 
      for(int i=0; i<9; i++) a_[i]=0.; 
  } 
  inline	Matrix3x3(const Matrix3x3 & m)
  { 
    for(int i=0; i<9; i++)  a_[i] = m.a_[i]; 
  }

  inline        ~Matrix3x3() { }

  inline Matrix3x3& Set(const Matrix3x3 & m)
                   { for(int i=0; i<9; i++) a_[i] = m.a_[i]; return (*this); }
  inline Matrix3x3& operator = (const Matrix3x3 & a)
                   { return Set(a); }

  inline double operator()(int r, int c) const 
                { return a_[r*3+c]; }
  inline double& operator()(int r, int c) 
                { return a_[r*3+c]; }

  inline int NRows() const {return 3; }
  inline int NCols() const {return 3; }

  inline Matrix3x3& AddElt(const Matrix3x3 & m)
                   { for(int i=0; i<9; i++) a_[i] += m.a_[i]; return (*this); }
  inline Matrix3x3& MulElt(const Matrix3x3 & m)
                   { for(int i=0; i<9; i++) a_[i] *= m.a_[i]; return (*this); }
  
  // Multiplication de deux Matrix3x3 entre elles  
/*    inline Matrix3x3& MulMat(const Matrix3x3 & m) */
/*      { */
/*        Matrix3x3 c; */
/*        // c.Set(0.); */
/*        for(int r=0; r<=2; r++) */
/*  	for(int c=0; c<=2; c++) */
/*  	  c.a_[r*3+c] = 0.; */
      
/*        for(int r=0; r<=2; r++) */
/*  	for(int c=0; c<=2; c++) */
/*  	  for(int k=0; k<=2; k++) */
/*  	    c.a_[3*r+c] += a_[3*r+k] * m.a_[3*k+c]; */
      
/*        a_ = c; */
/*        return(*this); */
/*      } */

  inline Matrix3x3& Transpose()
    {
      double c;
      c = a_[1]; a_[1] = a_[3]; a_[3] = c;
      c = a_[2]; a_[2] = a_[6]; a_[6] = c;
      c = a_[5]; a_[5] = a_[7]; a_[7] = c;
      return (*this);
    }

  inline Matrix3x3& Inverse()
    {
      // Calcul du determinant global
      double determ;
      determ =  a_[0]*a_[4]*a_[8] + a_[3]*a_[7]*a_[2] + a_[6]*a_[1]*a_[5];
      determ -= a_[2]*a_[4]*a_[6] + a_[5]*a_[7]*a_[0] + a_[8]*a_[1]*a_[3];
      // prevoir un break si determ <= threshold

      double am1_[9];
      // Calcul des cofacteurs
      am1_[0] =  (a_[4]*a_[8] - a_[5]*a_[7])/determ ; 
      am1_[1] = -(a_[3]*a_[8] - a_[5]*a_[6])/determ ; 
      am1_[2] =  (a_[3]*a_[7] - a_[4]*a_[6])/determ ; 
      am1_[3] = -(a_[1]*a_[8] - a_[7]*a_[2])/determ ; 
      am1_[4] =  (a_[0]*a_[8] - a_[6]*a_[2])/determ ; 
      am1_[5] = -(a_[0]*a_[7] - a_[6]*a_[1])/determ ; 
      am1_[6] =  (a_[1]*a_[5] - a_[4]*a_[2])/determ ; 
      am1_[7] = -(a_[0]*a_[5] - a_[3]*a_[2])/determ ; 
      am1_[8] =  (a_[0]*a_[4] - a_[3]*a_[1])/determ ; 
      /****************************************/
      /* c'est la version qui marche, mais la convention semble 
         differente de celle de [3*r + c]...? */


/*        // Transposition */
/*        am1_ = am1_.Transpose() ; */
/*        a_ = am1_; // marche pas bien : "request for member `Transpose' in `Matrix3x3::am1_', */
/*                   //which is of non-aggregate type 'double[9]' */

      // Calcul de l'inverse avec transpostion incluse
/*        am1_[0] =  (a_[4]*a_[8] - a_[5]*a_[7])/determ ; */
/*        am1_[3] = -(a_[3]*a_[8] - a_[5]*a_[6])/determ ; */
/*        am1_[6] =  (a_[3]*a_[7] - a_[4]*a_[6])/determ ; */
/*        am1_[1] = -(a_[1]*a_[8] - a_[7]*a_[2])/determ ; */
/*        am1_[4] =  (a_[0]*a_[8] - a_[6]*a_[2])/determ ; */
/*        am1_[7] = -(a_[0]*a_[7] - a_[6]*a_[1])/determ ; */
/*        am1_[2] =  (a_[1]*a_[5] - a_[4]*a_[2])/determ ; */
/*        am1_[5] = -(a_[0]*a_[5] - a_[3]*a_[2])/determ ; */
/*        am1_[8] =  (a_[0]*a_[4] - a_[3]*a_[1])/determ ; */
      for(int k=0; k<9; k++) a_[k] = am1_[k];
      return (*this);}


protected:
  double a_[9];
  
};


class RingOfMatrix3x3 {
public:
                RingOfMatrix3x3(int_4 npix);
  virtual      ~RingOfMatrix3x3();

  inline int_4  Size() { return size_; }
  inline Matrix3x3& GetElement(int_4 n) 
  {
    if ((n<0) || (n>=size_))
      throw RangeCheckError("RingOfMatrix3x3::GetElement() - Out of range index");
    return mtxtab_[n];
  }
  inline Matrix3x3& operator()(int_4 n) { return GetElement(n); }

protected:
  Matrix3x3 *mtxtab_;
  int_4 size_;  
};