Matrix3.hh

Go to the documentation of this file.

#pragma once

#include "boost/range/numeric.hpp"

#include "TQuaternion.h"

#include "boca/math/GradedVector3.hh"
#include "boca/math/Matrix2.hh"

namespace boca
{

template<typename Value_>
using Array3 = std::array<Value_, 3>;

template <typename Value_>
class Matrix3 : public boost::totally_ordered<Matrix3<Value_>>
            , boost::additive<Matrix3<Value_>>
{

    template<typename Value_2_>
    using ValueProduct = ValueProduct<Value_, Value_2_>;
    using ValueSquare = boca::ValueSquare<Value_>;
    using ValueCubed = boca::ValueCubed<Value_>;
    template<typename Value_2>
    using ValueQuotient = ValueQuotient<Value_, Value_2>;
    using ValueInverse = boost::units::divide_typeof_helper<double, Value_>;
    template<typename Value_2>
    using OnlyIfNotOrSameQuantity = typename std::enable_if < !IsQuantity<Value_2>::value || std::is_same<Value_, Value_2>::value >::type;
    template<typename Value_2>
    using OnlyIfNotQuantity = typename std::enable_if < !IsQuantity<Value_2>::value >::type;

public:

    constexpr Matrix3() {}

    Matrix3(Value_ scalar, Matrix matrix = Matrix::diagonal)
    {
        switch (matrix) {
        case Matrix::diagonal :
            SetDiagonal(scalar);
            break;
        case Matrix::uniform :
            SetUniform(scalar);
            break;
        default :
            std::cout << "Maformed matrix constructor: " << Name(matrix) << '\n';
        }
    }

    Matrix3(Vector3<Value_> const &x, Vector3<Value_> const &y, Vector3<Value_> const &z, Matrix matrix = Matrix::row)
    {
        switch (matrix) {
        case Matrix::row:
            SetRows(x, y, z);
            break;
        case Matrix::column :
            SetColumns(x, y, z);
            break;
        default :
            std::cout << "Maformed matrix constructor: " << Name(matrix) << '\n';
        }
    }

    Matrix3(Vector3<Value_> const &vector_1, Vector3<Value_> const &vector_2, Matrix matrix = Matrix::symmetric)
    {
        switch (matrix) {
        case Matrix::symmetric:
            *this = MatrixProduct(vector_1, vector_2);
            break;
        default :
            std::cout << "Maformed matrix constructor: " << Name(matrix) << '\n';
        }
    }

    Matrix3(Vector3<Value_> const &vector, Matrix matrix = Matrix::antisymmetric)
    {
        switch (matrix) {
        case Matrix::antisymmetric:
            SetAntisymmetric(vector);
            break;
        case Matrix::diagonal:
            SetDiagonal(vector);
            break;
        default :
            std::cout << "Maformed matrix constructor: " << Name(matrix) << '\n';
        }
    }

    Matrix3(Vector3<Value_> const &vector, Dim3 dim, Matrix matrix = Matrix::row)
    {
        switch (matrix) {
        case Matrix::row:
            SetRow(vector, dim);
            break;
        case Matrix::column:
            SetColumn(vector, dim);
            break;
        default :
            std::cout << "Maformed matrix constructor: " << Name(matrix) << '\n';
        }
    }



    Matrix3(TQuaternion const &quaternion)
    {
        auto const mag2 = quaternion.QMag2();
        if (mag2 <= 0) *this = Matrix3(1.);
        else {
            *this = 2 * (Matrix3(sqr(quaternion.fRealPart)) + Matrix3(Vector3<double>(quaternion.fVectorPart), Vector3<double>(quaternion.fVectorPart)) - Matrix3(Vector3<double>(quaternion.fVectorPart * quaternion.fRealPart)));
            // protect agains non-unit quaternion
            if (std::abs(mag2 - 1) > 1e-10) *this /= mag2;
            // diago : remove identity
            *this -= Matrix3(1);
        }
    }

    template<typename Value_2>
    constexpr Matrix3(Matrix3<Value_2> const &matrix) :
        x_(matrix.X()),
        y_(matrix.Y()),
        z_(matrix.Z())
    {}


    Matrix3 &SetDiagonal(Value_ value)
    {
        x_ = {value, Dim3::x};
        y_ = {value, Dim3::y};
        z_ = {value, Dim3::z};
        return *this;
    }

    Matrix3 &SetDiagonal(Vector3<Value_> vector)
    {
        x_ = {vector.X(), Dim3::x};
        y_ = {vector.Y(), Dim3::y};
        z_ = {vector.Z(), Dim3::z};
        return *this;
    }

    Matrix3 &SetUpperTriangle(Vector3<Value_> vector)
    {
        x_.Y() = vector.X();
        x_.Z() = vector.Y();
        y_.Z() = vector.Z();
        return *this;
    }

    Matrix3 &SetLowerTriangle(Vector3<Value_> vector)
    {
        y_.X() = vector.X();
        z_.X() = vector.Y();
        z_.Y() = vector.Z();
        return *this;
    }

    Matrix3 &SetToIdentity()
    {
        SetUpperTriangle({0});
        SetLowerTriangle({0});
        return SetDiagonal(1);
    }

    Matrix3 &SetUniform(Value_ value)
    {
        x_ = {value, value, value};
        y_ = {value, value, value};
        z_ = {value, value, value};
        return *this;
    }

    void SetRows(Vector3<Value_> const &x, Vector3<Value_> const &y, Vector3<Value_> const &z)
    {
        x_ = x;
        y_ = y;
        z_ = z;
    }

    void SetRow(Vector3<Value_> const &vector, Dim3 dim)
    {
        switch (dim) {
        case Dim3::x :
            x_ = vector;
            break;
        case Dim3::y :
            y_ = vector;
            break;
        case Dim3::z :
            z_ = vector;
            break;
        default :
            ;
        }
    }

    void SetColumns(Vector3<Value_> const &x, Vector3<Value_> const &y, Vector3<Value_> const &z)
    {
        x_ = {x.X(), y.X(), z.X()};
        y_ = {x.Y(), y.Y(), z.Y()};
        z_ = {x.Z(), y.Z(), z.Z()};
    }

    void SetColumn(Vector3<Value_> const &vector, Dim3 dim)
    {
        x_ = {vector.X(), dim};
        y_ = {vector.Y(), dim};
        z_ = {vector.Z(), dim};
    }

    void SetAntisymmetric(Vector3<Value_> const &vector)
    {
        x_ = {Value_(0), vector.Z(), -vector.Y()};
        y_ = { -vector.Z(), Value_(0), vector.X()};
        z_ = {vector.Y(), -vector.X(), Value_(0)};
    }

    Matrix3 &SetEulerAngles(Angle const &phi, Angle const &theta, Angle const &psi, Dim2 dim = Dim2::x)
    {
        SetToIdentity();
        Rotate(phi, Dim3::z);
        switch (dim) {
        case Dim2::x :
            Rotate(theta, Dim3::x);
            break;
        case Dim2::y :
            Rotate(theta, Dim3::y);
            break;
        default :
            ;
        }
        Rotate(psi, Dim3::z);
        return *this;
    }

    void SetPhi(Angle const &phi, Dim2 dim = Dim2::x)
    {
        SetEulerAngles(phi, Theta(dim), Psi(dim), dim);
    }

    void SetTheta(Angle const &theta, Dim2 dim = Dim2::x)
    {
        SetEulerAngles(Phi(dim), theta, Psi(dim), dim);
    }

    void SetPsi(Angle const &psi, Dim2 dim = Dim2::x)
    {
        SetEulerAngles(Phi(), Theta(dim), psi, dim);
    }

    Matrix3 &SetAxis(Vector3<Value_> const &axis, Dim3 dim)
    {
        return SetAxis(axis, {axis.Mag(), Next(dim)});
    }

    Matrix3 &SetAxis(Vector3<Value_> const &axis, Vector3<Value_> const &plane, Dim3 dim)
    {
        *this = MakeBasis(plane, axis).MoveTo(dim, Dim3::x).Transpose();
        return *this;
    }


    constexpr Vector3<Value_> const &X() const
    {
        return x_;
    }

    Vector3<Value_> &X()
    {
        return x_;
    }

    constexpr Vector3<Value_> const &Y() const
    {
        return y_;
    }

    Vector3<Value_> &Y()
    {
        return y_;
    }

    Vector3<Value_> &Z()
    {
        return z_;
    }

    constexpr Vector3<Value_> const &Z() const
    {
        return z_;
    }


    constexpr Matrix Symmetry()
    {
        if (UpperTriangle() == LowerTriangle()) return Matrix::symmetric;
        if (UpperTriangle() == -LowerTriangle() && Diagonal() == Vector3<Value_> {0, 0, 0}) return Matrix::antisymmetric;
        return Matrix::none;
    }

    constexpr Vector3<Value_> ColumnX() const
    {
        return {x_.X(), y_.X(), z_.X()};
    }

    constexpr Vector3<Value_> ColumnY() const
    {
        return {x_.Y(), y_.Y(), z_.Y()};
    }

    constexpr Vector3<Value_> ColumnZ() const
    {
        return {x_.Z(), y_.Z(), z_.Z()};
    }

    constexpr Vector3<Value_> Diagonal()
    {
        return {x_.X(), y_.Y(), z_.Z()};
    }

    constexpr Vector3<Value_> UpperTriangle()
    {
        return {x_.Y(), x_.Z(), y_.Z()};
    }

    constexpr Vector3<Value_> LowerTriangle()
    {
        return {y_.X(), z_.X(), z_.Y()};
    }


    constexpr auto Sign(Dim3 i, Dim3 j) const
    {
        return (static_cast<int>(i) + static_cast<int>(j)) % 2 ? -1 : 1;
    }

    constexpr auto Trace()const
    {
        return x_.X() + y_.Y() + z_.Z();
    }

    template<typename Value_2_>
    constexpr auto ProductTrace(Matrix3<Value_2_> const &matrix) const
    {
        return x_ * matrix.ColumnX() + y_ * matrix.ColumnY() + z_ * matrix.ColumnZ();
    }

    constexpr auto Determinant()const
    {
        return  boost::accumulate(Dimensions3(), ValueCubed(0), [&](ValueCubed & sum, Dim3 dim) {
            return sum + Laplace(Dim3::x, dim);
        });
    }

    constexpr auto Laplace(Dim3 dim_1, Dim3 dim_2) const
    {
        return (*this)[dim_1][dim_2] * Cofactor(dim_1, dim_2);
    }

    constexpr auto Cofactor(Dim3 dim_1, Dim3 dim_2) const
    {
        return static_cast<double>(Sign(dim_1, dim_2)) * Minor(dim_1, dim_2);
    }

    constexpr auto Minor(Dim3 delete_1, Dim3 delete_2) const
    {
        return SubMatrix(delete_1, delete_2).Determinant();
    }

    constexpr auto ReducedDeterminant(Dim3 dim_1, Dim3 dim_2) const
    {
        return Determinant() - Laplace(dim_1, dim_2);
    }

    Angle Phi(Dim3 dim) const
    {
        return Vector2<Value_>(x_(dim), y_(dim)).Phi();
    }

    Angle Theta(Dim3 dim) const
    {
        return acos(z_(dim));
    }

    Angle Phi(Dim2 dim) const
    {
        switch (dim) {
        case Dim2::x :
            return XPhi();
        case Dim2::y :
            return XPhi() + PiRad() / 2.;
        default :
            return XPhi();
        }
    }

    Angle Theta(Dim2 = Dim2::x) const
    {
        return Theta(Dim3::z);
    }

    Angle Psi(Dim2 dim) const
    {
        switch (dim) {
        case Dim2::x :
            return XPsi();
        case Dim2::y :
            return XPsi() - PiRad() / 2.;
        default :
            return XPhi();
        }
    }



    constexpr Matrix3 Transposed() const
    {
        return {ColumnX(), ColumnY(), ColumnZ()};
    }

    constexpr Matrix3 &Transpose()
    {
        return *this = Transposed();
    }

    constexpr Matrix3<ValueInverse> Inverse()
    {
        auto const det = Determinant();
        if (det == ValueCubed(0)) return {};
        auto const x = ColumnX();
        auto const y = ColumnY();
        auto const z = ColumnZ();
        return Matrix3<ValueSquare>(y ^ z, z ^ x, x ^ y) / det;
    }

    auto SubMatrix(Dim3 delete_1, Dim3 delete_2) const
    {
        auto dim2_1 = EnumIterator<Dim2> {Dim2::x};
        auto dim2_2 = EnumIterator<Dim2> {Dim2::x};
        auto matrix = Matrix2<Value_> {};
        for (auto dim3_1 : Dimensions3()) {
            if (dim3_1 == delete_1) continue;
            for (auto dim3_2 : Dimensions3()) {
                if (dim3_2 == delete_2) continue;
                matrix[*dim2_1][*dim2_2] = (*this)[dim3_1][dim3_2];
                ++dim2_2;
            }
            ++dim2_1;
            dim2_2.Set(Dim2::x);
        }
        return matrix;
    }

    constexpr Matrix3 &MoveTo(Dim3 dim_1, Dim3 dim_2)
    {
        std::swap((*this)(dim_1), (*this)(dim_2));
        std::swap((*this)(dim_1), (*this)(Third(dim_1, dim_2)));
        return *this;
    }


    Matrix3 &Rotate(Angle const &phi, Dim3 dim_1,  Dim3 dim_2)
    {
        if (phi == 0_rad)  return *this;
        auto const cos = boost::units::cos(phi);
        auto const sin = boost::units::sin(phi);
        auto const row = (*this)(dim_1);
        (*this)(dim_1) = cos * row - sin * (*this)(dim_2);
        (*this)(dim_2) = sin * row + cos * (*this)(dim_2);
        return *this;
    }

    constexpr Matrix3 &Rotate(Angle const &phi, Dim3 dim)
    {
        switch (dim) {
        case Dim3::x :
            return Rotate(phi, Dim3::y, Dim3::z);
        case Dim3::y :
            return Rotate(phi, Dim3::z, Dim3::x);
        case Dim3::z :
            return Rotate(phi, Dim3::x, Dim3::y);
        default :
            ;
        }
    }

    Matrix3 &Rotate(Angle const &phi, Vector3<Value_> const &axis)
    {
        if (phi == 0_rad)  return *this;
        auto const sin = boost::units::sin(phi);
        auto const cos = boost::units::cos(phi);
        auto const unit = axis.Unit();
        return Transform(Matrix3(cos) + Matrix3(unit, unit) * (1 - cos) - Matrix3(unit * sin));
    }

    constexpr Matrix3 &RotateAxes(Vector3<Value_> const &vector_x, Vector3<Value_> const &vector_y, Vector3<Value_> const &vector_z)
    {
        auto const epsilon = 0.001;
        auto const cross = vector_x.Cross(vector_y);
        if (std::abs(vector_z.X() - cross.X()) > epsilon || std::abs(vector_z.Y() - cross.Y()) > epsilon || std::abs(vector_z.Z() - cross.Z()) > epsilon || std::abs(vector_x.Mag2() - 1) > epsilon || std::abs(vector_y.Mag2() - 1) > epsilon || std::abs(vector_z.Mag2() - 1) > epsilon || std::abs(vector_x.Dot(vector_y)) > epsilon || std::abs(vector_y.Dot(vector_z)) > epsilon || std::abs(vector_z.Dot(vector_x)) > epsilon) {
            std::cout << "RotateAxes bad axis vectors" << '\n';
            return *this;
        }
        return Transform(Matrix3(vector_x.X(), vector_y.X(), vector_z.X(), vector_x.Y(), vector_y.Y(), vector_z.Y(), vector_x.Z(), vector_y.Z(), vector_z.Z()));
    }


    constexpr Matrix3 &RotateEulerAngles(Angle const &phi, Angle const &theta, Angle const &psi, Dim2 dim = Dim2::x)
    {
        auto euler = Matrix3 {};
        euler.SetEulerAngles(phi, theta, psi, dim);
        return Transform(euler);
    }

    std::pair<Vector3<Value_>, Angle> AngleAxis() const
    {
        auto const cosa  = 0.5 * (x_.X() + y_.Y() + z_.Z() - 1);
        auto const cosa1 = 1 - cosa;
        if (cosa1 <= Value_(0)) return std::make_pair<Vector3<Value_>, Angle>({0, 0, 1},  0_rad);
        auto x = x_.X() > cosa ? sqrt((x_.X() - cosa) / cosa1) : 0;
        auto y = y_.Y() > cosa ? sqrt((y_.Y() - cosa) / cosa1) : 0;
        auto z = z_.Z() > cosa ? sqrt((z_.Z() - cosa) / cosa1) : 0;
        if (z_.Y() < y_.Z()) x = -x;
        if (x_.Z() < z_.X()) y = -y;
        if (y_.X() < x_.Y()) z = -z;
        return std::make_pair<Vector3<Value_>, Angle>({x, y, z},  acos(cosa));
    }


    template<typename Value_2_>
    constexpr Matrix3<ValueProduct<Value_2_>> Multiply(Matrix3<Value_2_> const &matrix) const
    {
        return {{x_ * matrix.ColumnX(), x_ * matrix.ColumnY(), x_ * matrix.ColumnZ()}, {y_ * matrix.ColumnX(), y_ * matrix.ColumnY(), y_ * matrix.ColumnZ()}, {z_ * matrix.ColumnX(), z_ * matrix.ColumnY(), z_ * matrix.ColumnZ()}};
    }

    template<typename Value_2_>
    constexpr Vector3<ValueProduct<Value_2_>> Multiply(Vector3<Value_2_> const &vector) const
    {
        return {x_ * vector, y_ * vector, z_ * vector};
    }

    template<typename Value_2_>
    constexpr Matrix3<ValueProduct<Value_2_>> Scale(Value_2_ scalar) const
    {
        return {x_ * scalar, y_ * scalar, z_ * scalar};
    }

    template<typename Value_2_, typename = OnlyIfNotOrSameQuantity<Value_2_>>
    constexpr Matrix3<ValueProduct<Value_2_>> &Transform(Matrix3<Value_2_> const &matrix)
    {
        return *this = matrix * (*this);
    }

    template<typename Value_2_>
    constexpr Matrix3<ValueProduct<Value_2_>> &Transformed(Matrix3<Value_2_> const &matrix)
    {
        return matrix * (*this);
    }



    constexpr bool operator<(Matrix3 const &matrix) const
    {
        return abs(Determinant()) < abs(matrix.Determinant());
    }

    constexpr bool operator==(Matrix3 const &matrix) const
    {
        return x_ == matrix.x_ && y_ == matrix.y_ && z_ == matrix.z_;
    }

    template <typename Value_2_, typename = OnlyIfNotOrSameQuantity<Value_2_>>
    Matrix3 &operator+=(Matrix3<Value_2_> const &matrix)
    {
        x_ += matrix.x_;
        y_ += matrix.y_;
        z_ += matrix.z_;
        return *this;
    }

    template <typename Value_2_, typename = OnlyIfNotOrSameQuantity<Value_2_>>
    Matrix3 &operator-=(Matrix3<Value_2_> const &matrix)
    {
        x_ -= matrix.x_;
        y_ -= matrix.y_;
        z_ -= matrix.z_;
        return *this;
    }

    template<typename Value_2_, typename = OnlyIfNotQuantity<Value_2_>>
    constexpr Matrix3 &operator*=(Matrix3<Value_2_> const &matrix)
    {
        return *this = Multiply(matrix);
    }

    template<typename Value_2_>
    constexpr Matrix3<ValueProduct<Value_2_>> operator*(Matrix3<Value_2_> const &matrix) const
    {
        return Multiply(matrix);
    }

    template<typename Value_2_>
    constexpr Vector3<ValueProduct<Value_2_>> operator*(Vector3<Value_2_> const &vector) const
    {
        return Multiply(vector);
    }

    template<typename Value_2_, typename = OnlyIfNotOrSameQuantity<Value_2_>>
    Matrix3<ValueQuotient<Value_2_>> operator/=(Value_2_ scalar)
    {
        x_ /= scalar;
        y_ /= scalar;
        z_ /= scalar;
        return *this;
    }

    template<typename Value_2_>
    constexpr Matrix3<ValueQuotient<Value_2_>> operator/(Value_2_ scalar) const
    {
        return Scale(1. / scalar);
    }

    constexpr Vector3<Value_> const &operator()(Dim3 dim) const
    {
        switch (dim) {
        case Dim3::x :
            return x_;
        case Dim3::y :
            return y_;
        case Dim3::z :
            return z_;
        default :
            Default("Matrix3", to_int(dim));
            return x_;
        }
    }

    Vector3<Value_> &operator()(Dim3 dim)
    {
        switch (dim) {
        case Dim3::x :
            return x_;
        case Dim3::y :
            return y_;
        case Dim3::z :
            return z_;
        default :
            Default("Matrix3", to_int(dim));
            return x_;
        }
    }

    constexpr Value_ operator()(Dim3 i, Dim3 j) const
    {
        return operator()(i)(j);
    }

    Value_ &operator()(Dim3 i, Dim3 j)
    {
        return operator()(i)(j);
    }

    constexpr Vector3<Value_> const& operator[](Dim3 dim_3) const
    {
        switch (dim_3) {
        case Dim3::x :
            return x_;
        case Dim3::y :
            return y_;
        case Dim3::z :
            return z_;
        default :
            Default("Matrix3", to_int(dim_3));
            return x_;
        }
    }

    Vector3<Value_> &operator[](Dim3 dim_3)
    {
        return const_cast<Vector3<Value_> &>(static_cast<Matrix3<Value_> const &>(*this)[dim_3]);
    }

    friend auto &operator<<(std::ostream &stream, Matrix3<Value_> const &matrix)
    {
        for(auto const& vector : matrix) stream << vector;
        return stream;
    }


    using Dimension = Dim3;

    constexpr ConstSubIterator<boca::Matrix3, Vector3, Value_> begin() const
    {
        return {this, Dim3::x};
    }

    constexpr ConstSubIterator<boca::Matrix3, Vector3, Value_> end() const
    {
        return {this, Dim3::last};
    }

    SubIterator<boca::Matrix3, Vector3, Value_> begin()
    {
        return {this, Dim3::x};
    }

    SubIterator<boca::Matrix3, Vector3, Value_> end()
    {
        return {this, Dim3::last};
    }


    constexpr Array3<Value_> EigenValues() const
    {
        return Eigen().Values();
    }

    constexpr Array3<Vector3<Value_>> EigenVectors() const
    {
        return Eigen().Vectors();
    }

    constexpr Array3<GradedVector3<Value_>> EigenSystem() const
    {
        return Eigen().System();
    }


private:

    class Characteristic_
    {

    public:

        Characteristic_(Matrix3 const &matrix)
        {
            qudratic_ = - matrix.Trace();
            linear_ = - (matrix.ProductTrace(matrix) - sqr(qudratic_)) / 2;
            constant_ = - matrix.Determinant();
        }

        constexpr Value_ Qudratic() const
        {
            return qudratic_;
        }

        constexpr ValueSquare Linear() const
        {
            return linear_;
        }

        constexpr ValueCubed Constant() const
        {
            return constant_;
        }

    private:

        Value_ qudratic_;

        ValueSquare linear_;

        ValueCubed constant_;

    };

    class Depressed_
    {

    public:

        Depressed_(Characteristic_ const &c)
        {
            auto const quadratic = c.Qudratic() / 3;
            linear_ = c.Linear() - sqr(quadratic) * 3;
            constant_ = 2. * cube(quadratic) - quadratic * c.Linear() + c.Constant();
        }

        constexpr ValueSquare Linear() const
        {
            return linear_;
        }

        constexpr ValueCubed Constant() const
        {
            return constant_;
        }

    private:

        ValueSquare linear_;

        ValueCubed constant_;

    };

    class Eigen_
    {

    public:

        constexpr Eigen_() {}

        Eigen_(Depressed_ const &d, Matrix3<Value_> const &matrix)
        {
            factor_ = 2. * sqrt(- d.Linear() / 3.);
            angle_ = acos(3. * d.Constant() / d.Linear() / factor_);
            complex_ = 4 * cube(d.Linear()) > - 27 * sqr(d.Constant());
            matrix_ = &matrix;
        }

        constexpr Array3<Value_> Values() const
        {
            return values_.Get([&]() {
                auto values = Array3<Value_> {};
                if (complex_) {
                    for (auto &value : values) value = -1;
                    std::cerr << "Eigensystem has no real Eigenvalues!\n";
                    return values;
                }
                for (auto const index : IntegerRange(values.size())) values.at(index) = Value(index);
                return boost::range::sort(values, [](Value_ i, Value_ j) {
                    return abs(i) > abs(j);
                });
            });
        }

        constexpr Array3<Vector3<Value_>> Vectors() const
        {
            return vectors_.Get([&]() {
                auto vectors = Array3<Vector3<Value_>> {};
                for (auto const index : IntegerRange(vectors.size())) vectors.at(index) = Vector(index);
                return vectors;
            });
        }

        constexpr Array3<GradedVector3<Value_>> System() const
        {
            auto system = Array3<GradedVector3<Value_>> {};
            for (auto const index : IntegerRange(system.size())) system.at(index) = {Vectors().at(index), Values().at(index)};
            return system;
        }

    private:

        constexpr Value_ Factor() const
        {
            return factor_;
        }

        boca::Angle Angle() const
        {
            return angle_;
        }

        constexpr Value_ Value(int index) const
        {
            return Factor() * cos((Angle() - TwoPiRad() * static_cast<double>(index)) / 3.) + matrix_->Trace() / 3;
        }

        constexpr Vector3<Value_> Vector(int index) const
        {
            auto const matrix = *matrix_ - Matrix3<Value_>(Values().at(index));
            return Vector3<Value_>(matrix.Cofactor(Dim3::x, Dim3::x), matrix.Cofactor(Dim3::x, Dim3::y), matrix.Cofactor(Dim3::x, Dim3::z)).Unit();
        }

        Value_ factor_;

        boca::Angle angle_;

        bool complex_;

        Mutable<Array3<Value_>> values_;

        Mutable<Array3<Vector3<Value_>>> vectors_;

        Matrix3<Value_> const *matrix_;

    };

    constexpr Eigen_ const &Eigen() const
    {
        return eigen_.Get([&]() {
            return Eigen_(Depressed_(Characteristic_(*this)), *this);
        });
    }

    Mutable<Eigen_> eigen_;

    Vector3<Value_> x_;

    Vector3<Value_> y_;

    Vector3<Value_> z_;

    Angle XPhi() const
    {
        auto s2 =  ValueSquare(1) - z_.Z() * z_.Z();
        if (s2 < ValueSquare(0)) {
            std::cout << "Phi() |z_.Z()| > 1 " << '\n';
            s2 = 0;
        }
        auto const sinTheta = std::sqrt(s2);
        if (sinTheta != Value_(0)) {
            auto const cscTheta = 1 / sinTheta;
            auto const cosAbsPhi =  z_.Y() * cscTheta;
            if (std::abs(cosAbsPhi) > 1) {         // NaN-proofing
                std::cout << "Phi() finds | cos phi | > 1" << '\n';
                cosAbsPhi = 1;
            }
            auto const absPhi = acos(cosAbsPhi);
            if (z_.X() > 0) return absPhi;
            if (z_.X() < 0) return -absPhi;
            if (z_.Y() > 0) return 0_rad;
            return PiRad();
        } else {              // sinTheta == Value(0) so |Fzz| = 1
            auto const absPhi = .5 * acos(x_.X());
            if (x_.Y() > 0) return -absPhi;
            if (x_.Y() < 0) return absPhi;
            if (x_.X() > 0) return 0_rad;
            return z_.Z() * PiRad() / 2;
        }
    }

// Return the euler angles of the rotation.  See SetYEulerAngles for a
// note about conventions.
    Angle XPsi() const
    {
        // psi angle
        auto s2 =  1 - sqr(z_.Z());
        if (s2 < 0) {
            std::cout << "Psi() |z_.Z()| > 1 " << '\n';
            s2 = 0;
        }
        auto const sinTheta = std::sqrt(s2);
        if (sinTheta != 0) {
            auto const cscTheta = 1 / sinTheta;
            auto const cosAbsPsi =  - y_.Z() * cscTheta;
            if (std::abs(cosAbsPsi) > 1) {         // NaN-proofing
                std::cout << "Psi() | cos psi | > 1 " << '\n';
                cosAbsPsi = 1;
            }
            auto const absPsi = boca::acos(cosAbsPsi);
            if (x_.Z() > 0) return absPsi;
            if (x_.Z() < 0) return -absPsi;
            return (y_.Z() < 0) ? 0_rad : PiRad();
        } else {              // sinTheta == Value(0) so |Fzz| = 1
            auto absPsi = x_.X();
            // NaN-proofing
            if (std::abs(x_.X()) > 1) {
                std::cout << "Psi() | x_.X() | > 1 " << '\n';
                absPsi = 1;
            }
            auto const absPsi2 = .5 * acos(absPsi);
            if (y_.X() > 0) return absPsi2;
            if (y_.X() < 0) return -absPsi2;
            return (x_.X() > 0) ? 0_rad : PiRad() / 2.;
        }
    }

    Matrix3<double> MakeBasis(Vector3<Value_> const &plane, Vector3<Value_> const &axis) const
    {
        auto const mag_z = axis.Mag();
        if (mag_z == Value_(0)) {
            if (plane.Mag() != Value_(0)) return MakeBasis(axis, plane);
            return {};
        }
        auto matrix = Matrix3<double> {axis / mag_z, Dim3::z};

        auto mag_x = plane.Mag();
        if (mag_x == Value_(0)) {
            matrix.X() = matrix.Z().Orthogonal();
            mag_x = 1;
        } else matrix.X() = plane;

        matrix.Y() = matrix.Z().Cross(matrix.X()) / mag_x;
        auto const mag_y = matrix.Y().Mag();
        matrix.Y() = mag_y == Value_(0) ? matrix.Z().Orthogonal() : matrix.Y() / mag_y;
        matrix.X() = matrix.Y().Cross(matrix.Z());
        return matrix;
    }

};

template <typename>
struct IsMatrix3 : std::false_type { };

template <typename Value>
struct IsMatrix3<Matrix3<Value>> : std::true_type { };

template<typename Value>
using OnlyIfNotMatrix3 = typename std::enable_if < !IsMatrix3<Value>::value >::type;

template < class Value_, class Value_2_, typename = OnlyIfNotMatrix3<Value_2_> >
constexpr auto operator*(Matrix3<Value_> const &matrix, Value_2_ scalar)
{
    return matrix.Scale(scalar);
}

template < class Value_, class Value_2_, typename = OnlyIfNotMatrix3<Value_2_> >
constexpr auto operator*(Value_2_ scalar, Matrix3<Value_> const &matrix)
{
    return matrix.Scale(scalar);
}

template < class Value_, class Value_2_>
constexpr auto operator*(Matrix3<Value_> const &matrix, Vector3<Value_2_> const &vector)
{
    return matrix.Multiply(vector);
}

template < class Value_, class Value_2_>
constexpr Matrix3<ValueProduct<Value_,  Value_2_>> operator*(Vector3<Value_2_> const &vector, Matrix3<Value_> const &matrix)
{
    return {vector.X() *matrix.ColumnX(), vector.Y() *matrix.ColumnY(), vector.Z() *matrix.ColumnZ()};
}

template<typename Value_1_, typename Value_2_>
constexpr Matrix3<ValueProduct<Value_1_,  Value_2_>> MatrixProduct(Vector3<Value_1_> const &vector_1, Vector3<Value_2_> const &vector_2)
{
    return {vector_1.X() *vector_2, vector_1.Y() *vector_2, vector_1.Z() *vector_2};
}


template<typename Value_1_>
template<typename Value_2_>
auto &Vector3<Value_1_>::Rotate(boca::Angle angle, Vector3<Value_2_> const &axis)
{
    Matrix3<double> matrix;
    matrix.Rotate(angle, axis);
    operator*=(matrix);
}

template<typename Value_1_>
template<typename Value_2_>
auto &Vector3<Value_1_>::Transform(Matrix3<Value_2_> const &matrix)
{
    return *this = matrix * (*this);
}

template<typename Value_1_>
template<typename Value_2_>
auto &Vector3<Value_1_>::operator*=(Matrix3<Value_2_> const &matrix)
{
    return *this = matrix * (*this);
}

}

namespace boost
{

template<typename Value_>
struct range_const_iterator< boca::Matrix3<Value_> > {
    typedef typename boca::ConstSubIterator<boca::Matrix3, boca::Vector3, Value_> type;
};

template<typename Value_>
struct range_mutable_iterator< boca::Matrix3<Value_> > {
    typedef typename boca::SubIterator<boca::Matrix3, boca::Vector3, Value_> type;
};

}