//******************************************************************************
//** SCATMECH: Polarized Light Scattering C++ Class Library
//** 
//** File: matrix3d.h
//**
//** Thomas A. Germer
//** Optical Technology Division, National Institute of Standards and Technology
//** 100 Bureau Dr. Stop 8443; Gaithersburg, MD 20899-8443
//** Phone: (301) 975-2876; FAX: (301) 975-6991
//** Email: thomas.germer@nist.gov
//**
//** Version: 6.00 (February 2008)
//**
//******************************************************************************
#ifndef SCATMECH_MATRIX3D_H
#define SCATMECH_MATRIX3D_H

#include "vector3d.h"


namespace SCATMECH {


    // 
    // template class Matrix3D_Base is used by the class Matrix3D
    //
    template <class TYPE>
    class Matrix3D_Base {
        public:
            //The elements of the matrix:
            TYPE x[3][3];

            // Default constructor:
            Matrix3D_Base() {}; // <--This doesn't initialize array

            // Copy constructor:
            Matrix3D_Base(const Matrix3D_Base<TYPE>& a)
            {
                for (int i=0;i<3;++i)
                    for (int j=0;j<3;++j)
                        x[i][j]=a.x[i][j];
            }

            // Element by element constructor:
            Matrix3D_Base(const TYPE& a,const TYPE& b,const TYPE& c,
                     const TYPE& d,const TYPE& e,const TYPE& f,
                     const TYPE& g,const TYPE& h,const TYPE& i)
            {
                x[0][0]=a;
                x[0][1]=b;
                x[0][2]=c;
                x[1][0]=d;
                x[1][1]=e;
                x[1][2]=f;
                x[2][0]=g;
                x[2][1]=h;
                x[2][2]=i;
            }

            // Constructor which specifies three vectors:
            Matrix3D_Base(const Vector3D<TYPE>& a,const Vector3D<TYPE>& b,const Vector3D<TYPE>& c)
            {
                x[0][0]=a.x;
                x[0][1]=a.y;
                x[0][2]=a.z;
                x[1][0]=b.x;
                x[1][1]=b.y;
                x[1][2]=b.z;
                x[2][0]=c.x;
                x[2][1]=c.y;
                x[2][2]=c.z;  
            };

            // Assignment operator:
            Matrix3D_Base<TYPE>& operator=(const Matrix3D_Base<TYPE>& a)
            {
                for (int i=0;i<3;++i)
                    for (int j=0;j<3;++j)
                        x[i][j]=a.x[i][j];
                return *this;
            };
    
            // Addition of two matrices:
            Matrix3D_Base<TYPE> operator+(const Matrix3D_Base<TYPE>& a) const
            {
                Matrix3D_Base<TYPE> result;
                for (int i=0;i<3;++i)
                    for (int j=0;j<3;++j)
                        result.x[i][j]=x[i][j]+a.x[i][j];
                return result;
            };

            // Subtraction between two matrices:
            Matrix3D_Base<TYPE> operator-(const Matrix3D_Base<TYPE>& a) const
            {
                Matrix3D_Base<TYPE> result;
                for (int i=0;i<3;++i)
                    for (int j=0;j<3;++j)
                        result.x[i][j]=x[i][j]-a.x[i][j];
                return result;
            };

            // Unary minus sign:
            Matrix3D_Base<TYPE> operator-() const
            {
                Matrix3D_Base<TYPE> result;
                for (int i=0;i<3;++i)
                    for (int j=0;j<3;++j)
                        result.x[i][j]=-x[i][j];
                return result;
            };

            // Scalar multiplication of two matrices:
            Matrix3D_Base<TYPE> operator*(const Matrix3D_Base<TYPE>& a) const
            {
                Matrix3D_Base<TYPE> result;
                int i,j;
                for (i=0;i<3;++i)
                    for (j=0;j<3;++j)
                        result.x[i][j] = 0.;
                for (i=0;i<3;++i)
                    for (j=0;j<3;++j)
                        for (int k=0;k<3;++k)
                            result.x[i][k]+=x[i][j]*a.x[j][k];
                return result;
            }
        
            // Left multiplication of a vector by a matrix:
            Vector3D<TYPE> operator*(const Vector3D<TYPE>& a) const
            {
                return Vector3D<TYPE>(x[0][0]*a.x+x[0][1]*a.y+x[0][2]*a.z,
                                      x[1][0]*a.x+x[1][1]*a.y+x[1][2]*a.z,
                                      x[2][0]*a.x+x[2][1]*a.y+x[2][2]*a.z);
            };
        
            // Right multiplication of a vector by a matrix:
            friend Vector3D<TYPE> operator*(const Vector3D<TYPE>& b,const Matrix3D_Base<TYPE>& a)
            {
                return Vector3D<TYPE>(a.x[0][0]*b.x + a.x[1][0]*b.y + a.x[2][0]*b.z,
                                      a.x[0][1]*b.x + a.x[1][1]*b.y + a.x[2][1]*b.z,
                                      a.x[0][2]*b.x + a.x[1][2]*b.y + a.x[2][2]*b.z);   
            };

            // Multiplication of a matrix by a scalar:
            Matrix3D_Base<TYPE> operator*(TYPE b) const
            {
                Matrix3D_Base<TYPE> result;
                for (int i=0;i<3;++i)
                    for (int j=0;j<3;++j)
                        result.x[i][j]=x[i][j]*b;
                return result;
            };
            friend Matrix3D_Base<TYPE> operator*(TYPE b,const Matrix3D_Base<TYPE>& a)
                { return a*b; };
        
            // Division of matrix by a vector:
            Matrix3D_Base<TYPE> operator/(TYPE b) const { return (*this)*(1./b);}; 
        
            // The outer product of two vectors is a matrix:
            friend Matrix3D_Base<TYPE> outer(const Vector3D<TYPE>& a,const Vector3D<TYPE>& b)
            {
                Matrix3D_Base<TYPE> result;
                result.x[0][0]= a.x*b.x;
                result.x[0][1]= a.x*b.y;
                result.x[0][2]= a.x*b.z;
                result.x[1][0]= a.y*b.x;
                result.x[1][1]= a.y*b.y;
                result.x[1][2]= a.y*b.z;
                result.x[2][0]= a.z*b.x;
                result.x[2][1]= a.z*b.y;
                result.x[2][2]= a.z*b.z;
                return result;
            };

        protected:

            // The type-converting copy function:
            template <class TYPE2>
            void copy(const Matrix3D_Base<TYPE2>& a) {
                     for (int i=0;i<3;++i)
                        for (int j=0;j<3;++j)
                            x[i][j]=a.x[i][j];
            }

    };

    //
    // Define Matrix3D<TYPE>
    //
    template <class TYPE> 
    class Matrix3D : public Matrix3D_Base<TYPE> {};

    // 
    // Define Matrix3D<double>
    //
    template <>
    class Matrix3D<double> : public Matrix3D_Base<double> {
    public:
        typedef double TT;
        // Constructors:
        Matrix3D<TT>() {};
        Matrix3D<TT>(const Matrix3D_Base<TT>& a) : Matrix3D_Base<TT>(a) {};
        Matrix3D<TT>(TT a,TT b,TT c,TT d,TT e,TT f,TT g,TT h,TT i) 
            : Matrix3D_Base<TT>(a,b,c,d,e,f,g,h,i) {}
        Matrix3D<TT>(const Vector3D<TT>& a,const Vector3D<TT>& b,const Vector3D<TT>& c)
            : Matrix3D_Base<TT>(a,b,c) {}  
    };

    // 
    // Define Matrix3D<COMPLEX >
    //
    template <>
    class Matrix3D<COMPLEX > : public Matrix3D_Base<COMPLEX > {
    public:
        typedef COMPLEX TT;
        // Constructors:
        Matrix3D<TT>() {};
        Matrix3D<TT>(const Matrix3D_Base<TT>& a) : Matrix3D_Base<TT>(a) {};
        Matrix3D<TT>(const TT& a,const TT& b,const TT& c,
                   const TT& d,const TT& e,const TT& f,
                   const TT& g,const TT& h,const TT& i)
                   :Matrix3D_Base<TT>(a,b,c,d,e,f,g,h,i) {}
        Matrix3D<TT>(const Vector3D<TT>& a,const Vector3D<TT>& b,const Vector3D<TT>& c)
            : Matrix3D_Base<TT>(a,b,c) {}

        // Constructor which converts a real matrix to a complex matrix:
        Matrix3D<TT>(const Matrix3D_Base<double>& a) {copy(a);};

    };

    // The following typedefs are used by the SCATMECH library:
    typedef Matrix3D<double> Matrix;
    typedef Matrix3D< COMPLEX > CMatrix;



} // namespace SCATMECH 


#endif // MATRIX3D_H