#ifndef __JMATRIX3D__
#define __JMATRIX3D__

#include <cmath>

#include "JIO/JSerialisable.hh"
#include "JGeometry2D/JMatrix2D.hh"


namespace JGEOMETRY3D {

  namespace {
    using JIO::JReader;
    using JIO::JWriter;
    using JGEOMETRY2D::JMatrix2D;
  }


  /**
   * Data structure for an arbitrary 3x3 matrix.
   */
  class JMatrix3D :
    public JMatrix2D
  {
  public:
    /**
     * Default constructor (= identity matrix).
     */
    JMatrix3D() :
      JMatrix2D()
    {
      a00_ =  1;  a01_ = 0;  a02_ = 0;
      a10_ =  0;  a11_ = 1;  a12_ = 0;
      a20_ =  0;  a21_ = 0;  a22_ = 1;
    }
    
    double a02() const { return a02_; }
    double a12() const { return a12_; }
    double a20() const { return a20_; }
    double a21() const { return a21_; }
    double a22() const { return a22_; }

    
    /**
     * Transpose.
     */
    JMatrix3D& transpose()
    {
      std::swap(a10_,a01_);
      std::swap(a20_,a02_);
      std::swap(a21_,a12_);
      
      return *this;
    }


    /**
     * Check whether this is the identity matrix.
     *
     * \param  eps        precision
     */
    bool isIdentity(const double eps = 1.0e-10) const
    {
      return (fabs(1.0 - a00()) <= eps &&
	      fabs(1.0 - a11()) <= eps &&
	      fabs(1.0 - a22()) <= eps &&
	      fabs(a01())       <= eps &&
	      fabs(a02())       <= eps &&
	      fabs(a10())       <= eps &&
	      fabs(a12())       <= eps &&
	      fabs(a21())       <= eps &&
	      fabs(a20())       <= eps);
    }


    /**
     * Matrix multiplication.
     *
     * \param  B          JMatrix3D
     * \return            this matrix
     */
    JMatrix3D& multiply(const JMatrix3D& B)
    {
      const double c00 = a00() * B.a00()  +  a01() * B.a10()  +  a02() * B.a20();
      const double c01 = a00() * B.a01()  +  a01() * B.a11()  +  a02() * B.a21();
      const double c02 = a00() * B.a02()  +  a01() * B.a12()  +  a02() * B.a22();

      const double c10 = a10() * B.a00()  +  a11() * B.a10()  +  a12() * B.a20();
      const double c11 = a10() * B.a01()  +  a11() * B.a11()  +  a12() * B.a21();
      const double c12 = a10() * B.a02()  +  a11() * B.a12()  +  a12() * B.a22();

      const double c20 = a20() * B.a00()  +  a21() * B.a10()  +  a22() * B.a20();
      const double c21 = a20() * B.a01()  +  a21() * B.a11()  +  a22() * B.a21();
      const double c22 = a20() * B.a02()  +  a21() * B.a12()  +  a22() * B.a22();

      a00_ = c00;  a01_ = c01;  a02_ = c02; 
      a10_ = c10;  a11_ = c11;  a12_ = c12; 
      a20_ = c20;  a21_ = c21;  a22_ = c22; 

      return *this;
    }


    /**
     * Transform.
     *
     * \param  __x      x value
     * \param  __y      y value
     * \param  __z      z value
     */
    void transform(double& __x, double& __y, double& __z) const
    {
      const double x = a00() * __x  +  a01() * __y  +  a02() * __z;
      const double y = a10() * __x  +  a11() * __y  +  a12() * __z;
      const double z = a20() * __x  +  a21() * __y  +  a22() * __z;

      __x = x;
      __y = y;
      __z = z;
    }


    /**
     * Read JMatrix3D from input.
     *
     * \param  in         JReader
     * \param  matrix     JMatrix3D
     * \return            JReader
     */
    friend inline JReader& operator>>(JReader& in, JMatrix3D& matrix)
    {
      in >> matrix.a00_;  in >> matrix.a01_;  in >> matrix.a02_;
      in >> matrix.a10_;  in >> matrix.a11_;  in >> matrix.a12_;
      in >> matrix.a20_;  in >> matrix.a21_;  in >> matrix.a22_;

      return in;
    }


    /**
     * Write JMatrix3D to output.
     *
     * \param  out        JWriter
     * \param  matrix     JMatrix3D
     * \return            JWriter
     */
    friend inline JWriter& operator<<(JWriter& out, const JMatrix3D& matrix)
    {
      out << matrix.a00_;  out << matrix.a01_;  out << matrix.a02_;
      out << matrix.a10_;  out << matrix.a11_;  out << matrix.a12_;
      out << matrix.a20_;  out << matrix.a21_;  out << matrix.a22_;

      return out;
    }

  protected:
    double             a02_;
    double             a12_;
    double a20_, a21_, a22_;
  };
}

#endif