#ifndef Matrix4d4_HH
#define Matrix4d4_HH
//////////////////////////////////////////////////////////////////////////
//! file="amma/Geometry/AnGeo4/Matrix4d.hh"
//! lib=Mag4
//! userlevel=Normal
//! author="Radek Marik"
//! date="26.06.1995"
//! docentry="Geometry.4-D;Basic Types.Numerical.Specialised"
//! rcsid="$Id: Matrix4d.hh,v 1.9 2000/10/04 12:48:59 ees1cg Exp $"

#include "amma/StdType.hh"
#include "amma/Index.hh"
#include "amma/Boolean.hh"

class Point4dC;
class Vector4dC;
class MatrixC;
class ostream;
class istream;

// -------------------------------------------------------------------
// ******** Matrix4d4C ***********************************************
// -------------------------------------------------------------------

//: Matrix 4 x 4
// The class Matrix4d4C represents a matrix 4x4 of real numbers. The indexes
// start from 0.

class Matrix4d4C
{
  public:
    // Type definitions.
    // -----------------

    enum InputDataT {ROWS, COLUMNS};
    // The kind of input data.
    
  public:

    // Constructors, assigment, and destructor.
    // ----------------------------------------

    inline Matrix4d4C();
    // Creates the zero matrix.

    inline 
    Matrix4d4C(const RealT b00, const RealT b01,
                 const RealT b02, const RealT b03,
               const RealT b10, const RealT b11,
                 const RealT b12, const RealT b13,
               const RealT b20, const RealT b21,
                 const RealT b22, const RealT b23,
               const RealT b30, const RealT b31,
                 const RealT b32, const RealT b33);
    // Creates a new matrix and assign the values b[ij].

    Matrix4d4C(const Vector4dC & v0,
               const Vector4dC & v1,
               const Vector4dC & v2,
               const Vector4dC & v3,
               const InputDataT  dType = ROWS);
    // Creates a new matrix and assign the values according to the vectors.
    // If 'dType' is equal to ROWS, the vectors are treated as rows of
    // the matrix. If 'dType' is equal to COLUMNS, the vector are treated as
    // columns of the matrix.

    Matrix4d4C(const Point4dC & p0,
               const Point4dC & p1,
               const Point4dC & p2,
               const Point4dC & p3,
               const InputDataT dType = ROWS);
    // Creates a new matrix and assign the values according to the points
    // coordinates. If 'dType' is equal to ROWS, the points are treated
    // as rows of the matrix. If 'dType' is equal to COLUMNS, the points
    // are treated as columns of the matrix.
    
    Matrix4d4C(BooleanT);
    //: Creates a unit matrix.
    
    inline Matrix4d4C(const Matrix4d4C & mat);
    // Copy constructor.

    Matrix4d4C(const MatrixC & mat);
    // Copy constructor from the N-dimensional matrix.

    Matrix4d4C(istream & inS);
    // Creates a matrix form the input stream.

    inline const Matrix4d4C & operator=(const Matrix4d4C & mat);
    // Assigment.

    // Access to the matrix and its elements.
    // --------------------------------------

    inline IndexT RDim() const;
    // Returns the number of rows.

    inline IndexT CDim() const;
    // Returns the number of columns.

    inline const Matrix4d4C & Matrix() const;
    // Access to the constant matrix.

    inline Matrix4d4C & Matrix();                     
    // Access to the matrix.

    inline const RealT * operator[](IndexT i) const;
    // Returns the pointer to the row 'i' of the constant object.

    inline RealT * operator[](IndexT i);
    // Returns the pointer to the row 'i' of the object.

    inline RealT operator()(IndexT i, IndexT j) const; 
    // Returns the value of the element a[ij].

    inline RealT & operator()(IndexT i, IndexT j); 
    // Access to the value of the element a[ij].

    inline RealT A00() const;
    // Returns the value of the element a[00].

    inline RealT A01() const;
    // Returns the value of the element a[01].

    inline RealT A02() const;
    // Returns the value of the element a[02].

    inline RealT A03() const;
    // Returns the value of the element a[03].

    inline RealT A10() const;
    // Returns the value of the element a[10].

    inline RealT A11() const;
    // Returns the value of the element a[11].

    inline RealT A12() const;
    // Returns the value of the element a[12].

    inline RealT A13() const;
    // Returns the value of the element a[13].

    inline RealT A20() const;
    // Returns the value of the element a[20].

    inline RealT A21() const;
    // Returns the value of the element a[21].

    inline RealT A22() const;
    // Returns the value of the element a[22].

    inline RealT A23() const;
    // Returns the value of the element a[22].

    inline RealT A30() const;
    // Returns the value of the element a[30].

    inline RealT A31() const;
    // Returns the value of the element a[31].

    inline RealT A32() const;
    // Returns the value of the element a[32].

    inline RealT A33() const;
    // Returns the value of the element a[32].

    inline RealT & A00();
    // Access to the value of the element a[00].

    inline RealT & A01();
    // Access to the value of the element a[01].

    inline RealT & A02();
    // Access to the value of the element a[02].

    inline RealT & A03();
    // Access to the value of the element a[03].

    inline RealT & A10();
    // Access to the value of the element a[10].

    inline RealT & A11();
    // Access to the value of the element a[11].

    inline RealT & A12();
    // Access to the value of the element a[12].

    inline RealT & A13();
    // Access to the value of the element a[13].

    inline RealT & A20();
    // Access to the value of the element a[20].

    inline RealT & A21();
    // Access to the value of the element a[21].

    inline RealT & A22();
    // Access to the value of the element a[22].

    inline RealT & A23();
    // Access to the value of the element a[23].

    inline RealT & A30();
    // Access to the value of the element a[30].

    inline RealT & A31();
    // Access to the value of the element a[31].

    inline RealT & A32();
    // Access to the value of the element a[32].

    inline RealT & A33();
    // Access to the value of the element a[33].

    Vector4dC Column(IndexT i) const;
    // Returns the column 'i'.
    
    Vector4dC Row(IndexT i) const;
    // Return the row 'i'.

    Vector4dC Diagonal() const;
    // Returns the vector of main diagonal item.

    // Setting of matrix elements.
    // ---------------------------

    inline void Set(const Matrix4d4C & mat);
    // Sets all elements of the matrix according to the matrix 'mat'.

    inline Matrix4d4C & SetZero();
    // Sets all elements to be 0.0.

    Matrix4d4C & SetRow(IndexT i, const Point4dC & p);
    // Sets the elements of the row 'i' according to the elements of 
    // the point 'p'.

    Matrix4d4C & SetColumn(IndexT i, const Point4dC & p);
    // Sets the elements of the column 'i' according to the elements
    // of the point 'p'.

    // Logical operators.
    // ------------------

    inline BooleanT operator==(const Matrix4d4C & mat) const;
    // Returns TRUE if all elements of the matrix are equal 
    // to the elements of the matrix 'mat'.

    // Arithmetical operations.
    // ------------------------

    inline const Matrix4d4C & operator+=(const Matrix4d4C & mat);
    // Adds the matrix 'mat' to this matrix.

    inline const Matrix4d4C & operator-=(const Matrix4d4C & mat);
    // Subtracts the matrix 'mat' from this matrix.

    inline const Matrix4d4C & operator*=(RealT lambda);
    // Multiplies all elements by the scalar 'lambda'.

    inline const Matrix4d4C & operator/=(RealT lambda);
    // Divides all elements by the scalar 'lambda'.

    inline Matrix4d4C operator+(const Matrix4d4C & mat) const;
    // Sums both matrixes and returns the result.

    inline Matrix4d4C operator-(const Matrix4d4C & mat) const;
    // Subtracts the matrix 'mat' from this matrix and returns the result.

    inline Matrix4d4C operator*(RealT lambda) const;
    // Multiplies all elements of this matrix by the scalar 'lambda' and
    // returns the result.

    inline const Matrix4d4C & AddToDiagonal(RealT a);
    // The value 'a' is added to all diagonal elements.

    inline const Matrix4d4C & SubtractFromDiagonal(RealT a);
    // The value 'a' is subtracted from all diagonal elements.

    // Matrix operators.
    // -----------------

    inline Matrix4d4C T() const;
    // Returns the matrix, which is the transposition of this matrix.

    inline Matrix4d4C I() const;
    // Returns the inversion of this matrix.

    inline RealT Det() const;
    // Returns the determinant of the matrix.

    Point4dC operator*(const Point4dC & point) const;
    // Multiplies this matrix by the point 'p'. The returned point
    // is equal to (*this) * p.

    Vector4dC operator*(const Vector4dC & vector) const;
    // Multiplies the matrix by the vector 'v'. The returned vector
    // is equal to (*this) * v.

    inline Matrix4d4C operator*(const Matrix4d4C & mat) const;
    // Multiplies this matrix by the matrix 'mat'. The result matrix
    // is returned.

    Point4dC LinearSolver(const Point4dC & rightSide) const;
    // Returns the solution of 4 linear equations determined 
    // by the matrix and the right side vector 'rightSide'

    static Matrix4d4C ProductT(const Vector4dC & a);
    // Returns the matrix that is result of the tensor product a*a.T().

    // Error checking.
    // ---------------

    inline void ErrOutOfRange(IndexT i, const char * functionName) const;
    // Triggers the error handling if the index 'i' is smaller than 0
    // or bigger than 3. The name of the calling function 'functionName'
    // is passed to the error hangler.

  private:

    // Special functions.
    // ------------------

    static inline 
    RealT AdjointDet(RealT b00, RealT b01, RealT b02,
                     RealT b10, RealT b11, RealT b12,
                     RealT b20, RealT b21, RealT b22);
    // Computes the determinant of the matrix 3x3 created from b[ij] 
    // elements.
                       
  private:

    // Object representation.
    // ======================
    enum {dim = 4}; // the dimension of matrix

    RealT a[dim][dim]; // The matrix elements.
};    

ostream & operator<<(ostream & outS, const Matrix4d4C & m);
// Saves the matrix 'm' into the output stream 'outS'.

istream & operator>>(istream & inS, Matrix4d4C & matrix);
// Sets the matrix 'm' elements according to the values read from
// input stream 'inS'.

#include "amma/Error.hh"

// -----------------------------------------------------------------
// ******** Matrix4d4C *********************************************
// -----------------------------------------------------------------

inline 
void 
Matrix4d4C::ErrOutOfRange(IndexT i,
                          const char * functionName) const
{
  if (i < 0 || i > (dim-1))
  {
    errAMMA << "the index out of range <0," << (dim-1) << ">: " << i;
    errAMMA.Function(functionName).Exit();
  }
}

inline 
IndexT 
Matrix4d4C::RDim() const
{ return (IndexT) dim; }

inline 
IndexT 
Matrix4d4C::CDim() const
{ return (IndexT) dim; }

inline 
const RealT * 
Matrix4d4C::operator[](IndexT i) const
{
#ifdef AMMA_CHECK
  ErrOutOfRange(i, "Matrix4d4C::operator[](IndexT i) const");
#endif
  return a[i.V()];
}

inline 
RealT * 
Matrix4d4C::operator[](IndexT i)
{
#ifdef AMMA_CHECK
  ErrOutOfRange(i, "Matrix4d4C::operator[](IndexT i) const");
#endif
  return a[i.V()];
}

inline 
RealT
Matrix4d4C::operator()(IndexT i, IndexT j) const
{
  #ifdef AMMA_CHECK
    ErrOutOfRange(i, "Matrix4d4C::operator()(IndexT i, IndexT) const");
    ErrOutOfRange(j, "Matrix4d4C::operator()(IndexT, IndexT j) const");
  #endif
  return a[i.V()][j.V()];
}

inline 
RealT &
Matrix4d4C::operator()(IndexT i, IndexT j)
{
#ifdef AMMA_CHECK
  ErrOutOfRange(i, "Matrix4d4C::operator()(IndexT i, IndexT) const");
  ErrOutOfRange(j, "Matrix4d4C::operator()(IndexT, IndexT j) const");
#endif
  return a[i.V()][j.V()];
}

inline 
RealT 
Matrix4d4C::A00() const
{ return a[0][0]; }

inline 
RealT 
Matrix4d4C::A01() const
{ return a[0][1]; }

inline 
RealT 
Matrix4d4C::A02() const
{ return a[0][2]; }

inline 
RealT 
Matrix4d4C::A03() const
{ return a[0][3]; }

inline 
RealT 
Matrix4d4C::A10() const
{ return a[1][0]; }

inline 
RealT 
Matrix4d4C::A11() const
{ return a[1][1]; }

inline 
RealT 
Matrix4d4C::A12() const
{ return a[1][2]; }

inline 
RealT 
Matrix4d4C::A13() const
{ return a[1][3]; }

inline 
RealT 
Matrix4d4C::A20() const
//====================
{
  return a[2][0];
}

inline 
RealT 
Matrix4d4C::A21() const
{ return a[2][1]; }

inline 
RealT 
Matrix4d4C::A22() const
{ return a[2][2]; }

inline 
RealT 
Matrix4d4C::A23() const
{ return a[2][3]; }

inline 
RealT 
Matrix4d4C::A30() const
{ return a[3][0]; }

inline 
RealT 
Matrix4d4C::A31() const
{ return a[3][1]; }

inline 
RealT 
Matrix4d4C::A32() const
{ return a[3][2]; }

inline 
RealT 
Matrix4d4C::A33() const
{ return a[3][3]; }

inline 
RealT & 
Matrix4d4C::A00()
{ return a[0][0]; }

inline 
RealT & 
Matrix4d4C::A01()
{ return a[0][1]; }

inline 
RealT & 
Matrix4d4C::A02()
{ return a[0][2]; }

inline 
RealT & 
Matrix4d4C::A03()
//===============
{
  return a[0][3];
}

inline 
RealT & 
Matrix4d4C::A10()
{ return a[1][0]; }

inline 
RealT & 
Matrix4d4C::A11()
{ return a[1][1]; }

inline 
RealT & 
Matrix4d4C::A12()
{ return a[1][2]; }

inline 
RealT & 
Matrix4d4C::A13()
{ return a[1][3]; }

inline 
RealT & 
Matrix4d4C::A20()
{ return a[2][0]; }

inline 
RealT & 
Matrix4d4C::A21()
{ return a[2][1]; }

inline 
RealT & 
Matrix4d4C::A22()
//===============
{
  return a[2][2];
}

inline 
RealT & 
Matrix4d4C::A23()
{ return a[2][3]; }

inline 
RealT & 
Matrix4d4C::A30()
{ return a[3][0]; }

inline 
RealT & 
Matrix4d4C::A31()
{ return a[3][1]; }

inline 
RealT & 
Matrix4d4C::A32()
{ return a[3][2]; }

inline 
RealT & 
Matrix4d4C::A33()
//===============
{
  return a[3][3];
}


inline
Matrix4d4C::Matrix4d4C(const RealT b00,
                       const RealT b01,
                       const RealT b02,
                       const RealT b03,

                       const RealT b10,
                       const RealT b11,
                       const RealT b12,
                       const RealT b13,

                       const RealT b20,
                       const RealT b21,
                       const RealT b22,
                       const RealT b23,

                       const RealT b30,
                       const RealT b31,
                       const RealT b32,
                       const RealT b33
                      )
{
  A00() = b00;   A01() = b01;   A02() = b02;   A03() = b03;
  A10() = b10;   A11() = b11;   A12() = b12;   A13() = b13;
  A20() = b20;   A21() = b21;   A22() = b22;   A23() = b23;
  A30() = b30;   A31() = b31;   A32() = b32;   A33() = b33;
}

inline
Matrix4d4C::Matrix4d4C(const Matrix4d4C & mat)
{
  A00() = mat.A00(); A01() = mat.A01(); A02() = mat.A02(); A03() = mat.A03();
  A10() = mat.A10(); A11() = mat.A11(); A12() = mat.A12(); A13() = mat.A13();
  A20() = mat.A20(); A21() = mat.A21(); A22() = mat.A22(); A23() = mat.A23();
  A30() = mat.A30(); A31() = mat.A31(); A32() = mat.A32(); A33() = mat.A33();
}

inline
const Matrix4d4C & 
Matrix4d4C::operator=(const Matrix4d4C & mat)
{
  A00() = mat.A00();
  A01() = mat.A01();
  A02() = mat.A02();
  A03() = mat.A03();

  A10() = mat.A10();
  A11() = mat.A11();
  A12() = mat.A12();
  A13() = mat.A13();

  A20() = mat.A20();
  A21() = mat.A21();
  A22() = mat.A22();
  A23() = mat.A23();

  A30() = mat.A30();
  A31() = mat.A31();
  A32() = mat.A32();
  A33() = mat.A33();

  return *this;
}

inline
void 
Matrix4d4C::Set(const Matrix4d4C & mat)
{ *this = mat; }

inline
Matrix4d4C &
Matrix4d4C::SetZero()
{
  A00() = 0.0;
  A01() = 0.0;
  A02() = 0.0;
  A03() = 0.0;

  A10() = 0.0;
  A11() = 0.0;
  A12() = 0.0;
  A13() = 0.0;

  A20() = 0.0;
  A21() = 0.0;
  A22() = 0.0;
  A23() = 0.0;

  A30() = 0.0;
  A31() = 0.0;
  A32() = 0.0;
  A33() = 0.0;

  return *this;
}

inline
Matrix4d4C::Matrix4d4C()
{ SetZero(); }

inline 
const Matrix4d4C &
Matrix4d4C::Matrix() const
{ return *this; }

inline 
Matrix4d4C & 
Matrix4d4C::Matrix()                     
{ return *this; }

inline 
BooleanT 
Matrix4d4C::operator==(const Matrix4d4C & mat) const
{
  return    (A00() == mat.A00()) && (A01() == mat.A01())
         && (A02() == mat.A02()) && (A03() == mat.A03())
 
         && (A10() == mat.A10()) && (A11() == mat.A11())
         && (A12() == mat.A12()) && (A13() == mat.A13())

         && (A20() == mat.A20()) && (A21() == mat.A21())
         && (A22() == mat.A22()) && (A23() == mat.A23())

         && (A30() == mat.A30()) && (A31() == mat.A31())
         && (A32() == mat.A32()) && (A33() == mat.A33());
}
  
inline 
const Matrix4d4C & 
Matrix4d4C::operator+=(const Matrix4d4C & m)
{
  A00() += m.A00(); A01() += m.A01(); A02() += m.A02(); A03() += m.A03();
  A10() += m.A10(); A11() += m.A11(); A12() += m.A12(); A13() += m.A13();
  A20() += m.A20(); A21() += m.A21(); A22() += m.A22(); A23() += m.A23();
  A30() += m.A30(); A31() += m.A31(); A32() += m.A32(); A33() += m.A33();
  return *this;
}

inline 
const Matrix4d4C & 
Matrix4d4C::operator-=(const Matrix4d4C & m)
{
  A00() -= m.A00(); A01() -= m.A01(); A02() -= m.A02(); A03() -= m.A03();
  A10() -= m.A10(); A11() -= m.A11(); A12() -= m.A12(); A13() -= m.A13();
  A20() -= m.A20(); A21() -= m.A21(); A22() -= m.A22(); A23() -= m.A23();
  A30() -= m.A30(); A31() -= m.A31(); A32() -= m.A32(); A33() -= m.A33();
  return *this;
}

inline 
const Matrix4d4C & 
Matrix4d4C::operator*=(RealT lambda)
{
  A00() *= lambda;  A01() *= lambda;  A02() *= lambda;  A03() *= lambda;
  A10() *= lambda;  A11() *= lambda;  A12() *= lambda;  A13() *= lambda;
  A20() *= lambda;  A21() *= lambda;  A22() *= lambda;  A23() *= lambda;
  A30() *= lambda;  A31() *= lambda;  A32() *= lambda;  A33() *= lambda;
  return *this;
}

inline 
const Matrix4d4C & 
Matrix4d4C::operator/=(RealT lambda)
{
  A00() /= lambda;  A01() /= lambda;  A02() /= lambda;  A03() /= lambda;
  A10() /= lambda;  A11() /= lambda;  A12() /= lambda;  A13() /= lambda;
  A20() /= lambda;  A21() /= lambda;  A22() /= lambda;  A23() /= lambda;
  A30() /= lambda;  A31() /= lambda;  A32() /= lambda;  A33() /= lambda;
  return *this;
}

inline 
Matrix4d4C 
Matrix4d4C::operator+(const Matrix4d4C & m) const
{
  return Matrix4d4C(A00()+m.A00(), A01()+m.A01(), A02()+m.A02(), A03()+m.A03(),
                    A10()+m.A10(), A11()+m.A11(), A12()+m.A12(), A13()+m.A13(),
                    A20()+m.A20(), A21()+m.A21(), A22()+m.A22(), A23()+m.A23(),
                    A30()+m.A30(), A31()+m.A31(), A32()+m.A32(), A33()+m.A33()
                   );
}

inline 
Matrix4d4C 
Matrix4d4C::operator-(const Matrix4d4C & m) const
{
  return Matrix4d4C(A00()-m.A00(), A01()-m.A01(), A02()-m.A02(), A03()-m.A03(),
                    A10()-m.A10(), A11()-m.A11(), A12()-m.A12(), A13()-m.A13(),
                    A20()-m.A20(), A21()-m.A21(), A22()-m.A22(), A23()-m.A23(),
                    A30()-m.A30(), A31()-m.A31(), A32()-m.A32(), A33()-m.A33()
                   );
}

inline 
Matrix4d4C 
Matrix4d4C::operator*(RealT lambda) const
{
  return Matrix4d4C(A00()*lambda, A01()*lambda, A02()*lambda, A03()*lambda,
                    A10()*lambda, A11()*lambda, A12()*lambda, A13()*lambda,
                    A20()*lambda, A21()*lambda, A22()*lambda, A23()*lambda,
                    A30()*lambda, A31()*lambda, A32()*lambda, A33()*lambda
                   );
}

inline 
Matrix4d4C 
Matrix4d4C::operator*(const Matrix4d4C & m) const
{
  return Matrix4d4C(A00()*m.A00()+A01()*m.A10()+A02()*m.A20()+A03()*m.A30(),
                    A00()*m.A01()+A01()*m.A11()+A02()*m.A21()+A03()*m.A31(),
                    A00()*m.A02()+A01()*m.A12()+A02()*m.A22()+A03()*m.A32(),
                    A00()*m.A03()+A01()*m.A13()+A02()*m.A23()+A03()*m.A33(),
 
                    A10()*m.A00()+A11()*m.A10()+A12()*m.A20()+A13()*m.A30(),
                    A10()*m.A01()+A11()*m.A11()+A12()*m.A21()+A13()*m.A31(),
                    A10()*m.A02()+A11()*m.A12()+A12()*m.A22()+A13()*m.A32(),
                    A10()*m.A03()+A11()*m.A13()+A12()*m.A23()+A13()*m.A33(),
 
                    A20()*m.A00()+A21()*m.A10()+A22()*m.A20()+A23()*m.A30(),
                    A20()*m.A01()+A21()*m.A11()+A22()*m.A21()+A23()*m.A31(),
                    A20()*m.A02()+A21()*m.A12()+A22()*m.A22()+A23()*m.A32(),
                    A20()*m.A03()+A21()*m.A13()+A22()*m.A23()+A23()*m.A33(),
 
                    A30()*m.A00()+A31()*m.A10()+A32()*m.A20()+A33()*m.A30(),
                    A30()*m.A01()+A31()*m.A11()+A32()*m.A21()+A33()*m.A31(),
                    A30()*m.A02()+A31()*m.A12()+A32()*m.A22()+A33()*m.A32(),
                    A30()*m.A03()+A31()*m.A13()+A32()*m.A23()+A33()*m.A33()
                   );
}

inline 
const Matrix4d4C & 
Matrix4d4C::AddToDiagonal(RealT a)
{
  A00() += a;
  A11() += a;
  A22() += a;
  A33() += a;
  return *this;
}

inline 
const Matrix4d4C &
Matrix4d4C::SubtractFromDiagonal(RealT a)
{
  A00() -= a;
  A11() -= a;
  A22() -= a;
  A33() -= a;
  return *this;
}

inline 
Matrix4d4C 
Matrix4d4C::T() const
{
  return Matrix4d4C(A00(), A10(), A20(), A30(),
                    A01(), A11(), A21(), A31(),
                    A02(), A12(), A22(), A32(),
                    A03(), A13(), A23(), A33()
                   );
}

inline 
RealT 
Matrix4d4C::AdjointDet(
  RealT b00, RealT b01, RealT b02,
  RealT b10, RealT b11, RealT b12,
  RealT b20, RealT b21, RealT b22)
{
  return    b00 * (b11*b22 - b12*b21)
          + b10 * (b21*b02 - b22*b01)
          + b20 * (b01*b12 - b02*b11);
}

inline 
RealT
Matrix4d4C::Det() const
// The expansion is probably the fastest way.
{
  return    A00() * AdjointDet(A11(), A12(), A13(),
                               A21(), A22(), A23(),
                               A31(), A32(), A33())
          - A01() * AdjointDet(A10(), A12(), A13(),
                               A20(), A22(), A23(),
                               A30(), A32(), A33())
          + A02() * AdjointDet(A10(), A11(), A13(),
                               A20(), A21(), A23(),
                               A30(), A31(), A33())
          - A03() * AdjointDet(A10(), A11(), A12(),
                               A20(), A21(), A22(),
                               A30(), A31(), A32());
}

#endif

// IAPS - Image analysis program system. 
// End of include file Matrix4d.hh