text:
# We have rotations along one, two or three axes (e.g. RotX, RotXY and RotXYZ).
# They compose together nicely, so the user can make Euler angles by:
#
#     RotX(θx) * RotY(θy) * RotZ(θz) -> RotXYZ(θx, θy, θz)
#
# and never get confused by the order of application, etc.

#########################
# Single axis rotations #
#########################

for axis in [:X, :Y, :Z]
    RotType = Symbol("Rot" * string(axis))
    @eval begin
        struct $RotType{T} <: Rotation{3,T}
            theta::T
            $RotType{T}(theta) where {T} = new{T}(theta)
            $RotType{T}(r::$RotType) where {T} = new{T}(r.theta)
        end

        @inline $RotType(theta::T) where {T} = $RotType{T}(theta)
        @inline $RotType(r::$RotType{T}) where {T} = $RotType{T}(r)

        @inline (::Type{R})(t::NTuple{9}) where {R<:$RotType} = error("Cannot construct a cardinal axis rotation from a matrix")

        @inline Base.:*(r1::$RotType, r2::$RotType) = $RotType(r1.theta + r2.theta)

        @inline Base.inv(r::$RotType) = $RotType(-r.theta)

        # define null rotations for convenience
        @inline Base.one(::Type{$RotType}) = $RotType(0.0)
        @inline Base.one(::Type{$RotType{T}}) where {T} = $RotType{T}(zero(T))
    end
end

function Base.rand(::Type{R}) where R <: Union{RotX,RotY,RotZ}
    T = eltype(R)
    if T == Any
        T = Float64
    end

    return R(2*pi*rand(T))
end

"""
    struct RotX{T} <: Rotation{3,T}
    RotX(theta)

A 3×3 rotation matrix which represents a rotation by `theta` about the X axis.
"""
RotX

@inline function Base.getindex(r::RotX{T}, i::Int) where T
    T2 = Base.promote_op(sin, T)
    if i == 1
        one(T2)
    elseif i < 5
        zero(T2)
    elseif i == 5
        cos(r.theta)
    elseif i == 6
        sin(r.theta)
    elseif i == 7
        zero(T2)
    elseif i == 8
        -sin(r.theta)
    elseif i == 9
        cos(r.theta)
    else
        throw(BoundsError(r,i))
    end
end

@inline function Base.Tuple(r::RotX{T}) where T
    s, c = sincos(r.theta)
    o = one(s)
    z = zero(s)
    (o,  z,  z,   # transposed representation
     z,  c,  s,
     z, -s,  c)
end

@inline function Base.:*(r::RotX, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    st, ct = sincos(r.theta)
    T = Base.promote_op(*, typeof(st), eltype(v))
    return similar_type(v,T)(v[1],
                             v[2] * ct - v[3] * st,
                             v[3] * ct + v[2] * st)
end

"""
    struct RotY{T} <: Rotation{3,T}
    RotY(theta)

A 3×3 rotation matrix which represents a rotation by `theta` about the Y axis.
"""
RotY

@inline function Base.getindex(r::RotY{T}, i::Int) where T
    T2 = Base.promote_op(sin, T)
    if i == 1
        cos(r.theta)
    elseif i == 2
        zero(T2)
    elseif i == 3
        -sin(r.theta)
    elseif i == 4
        zero(T2)
    elseif i == 5
        one(T2)
    elseif i == 6
        zero(T2)
    elseif i == 7
        sin(r.theta)
    elseif i == 8
        zero(T2)
    elseif i == 9
        cos(r.theta)
    else
        throw(BoundsError(r,i))
    end
end

@inline function Base.Tuple(r::RotY{T}) where T
    s, c = sincos(r.theta)
    o = one(s)
    z = zero(s)
    (c,  z, -s,   # transposed representation
     z,  o,  z,
     s,  z,  c)
end

@inline function Base.:*(r::RotY, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    st, ct = sincos(r.theta)
    T = Base.promote_op(*, typeof(st), eltype(v))
    return similar_type(v,T)(v[1] * ct + v[3] * st,
                             v[2],
                             v[3] * ct - v[1] * st)
end

"""
    struct RotZ{T} <: Rotation{3,T}
    RotZ(theta)

A 3×3 rotation matrix which represents a rotation by `theta` about the Z axis.
"""
RotZ

@inline function Base.getindex(r::RotZ{T}, i::Int) where T
    T2 = Base.promote_op(sin, T)
    if i == 1
        cos(r.theta)
    elseif i == 2
        sin(r.theta)
    elseif i == 3
        zero(T2)
    elseif i == 4
        -sin(r.theta)
    elseif i == 5
        cos(r.theta)
    elseif i < 9
        zero(T2)
    elseif i == 9
        one(T2)
    else
        throw(BoundsError(r,i))
    end
end

@inline function Base.Tuple(r::RotZ{T}) where T
    s, c = sincos(r.theta)
    o = one(s)
    z = zero(s)
    ( c, s, z,   # transposed representation
     -s, c, z,
      z, z, o)
end

@inline function Base.:*(r::RotZ, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    st, ct = sincos(r.theta)
    T = Base.promote_op(*, typeof(st), eltype(v))
    return similar_type(v,T)(v[1] * ct - v[2] * st,
                             v[2] * ct + v[1] * st,
                             v[3])
end

################################################################################
################################################################################

######################
# Two axis rotations #
######################

for axis1 in [:X, :Y, :Z]
    Rot1Type = Symbol("Rot" * string(axis1))
    for axis2 in filter(axis -> axis != axis1, [:X, :Y, :Z])
        Rot2Type = Symbol("Rot" * string(axis2))
        RotType = Symbol("Rot" * string(axis1) * string(axis2))
        InvRotType = Symbol("Rot" * string(axis2) * string(axis1))

        @eval begin
            struct $RotType{T} <: Rotation{3,T}
                theta1::T
                theta2::T
                $RotType{T}(theta1, theta2) where {T} = new{T}(theta1, theta2)
                $RotType{T}(r::$RotType) where {T} = new{T}(r.theta1, r.theta2)
            end

            @inline $RotType(theta1::T1, theta2::T2) where {T1, T2} = $RotType{promote_type(T1, T2)}(theta1, theta2)
            @inline $RotType(r::$RotType{T}) where {T} = $RotType{T}(r)

            @inline function Base.getindex(r::$RotType{T}, i::Int) where T
                Tuple(r)[i] # Slow...
            end

            @inline (::Type{R})(t::NTuple{9}) where {R<:$RotType} = error("Cannot construct a two-axis rotation from a matrix")

            # Composing single-axis rotations to obtain a two-axis rotation:
            @inline Base.:*(r1::$Rot1Type, r2::$Rot2Type) = $RotType(r1.theta, r2.theta)

            # Composing single-axis rotations with two-axis rotations:
            @inline Base.:*(r1::$RotType, r2::$Rot2Type) = $RotType(r1.theta1, r1.theta2 + r2.theta)
            @inline Base.:*(r1::$Rot1Type, r2::$RotType) = $RotType(r1.theta + r2.theta1, r2.theta2)

            @inline Base.inv(r::$RotType) = $InvRotType(-r.theta2, -r.theta1)

            # define null rotations for convenience
            @inline Base.one(::Type{$RotType}) = $RotType(0.0, 0.0)
            @inline Base.one(::Type{$RotType{T}}) where {T} = $RotType{T}(zero(T), zero(T))
        end
    end
end

function Base.rand(::Type{R}) where R <: Union{RotXY,RotYZ,RotZX, RotXZ, RotYX, RotZY}
    T = eltype(R)
    if T == Any
        T = Float64
    end

    # Not really sure what this distribution is, but it's also not clear what
    # it should be! rand(RotXY) *is* invariant to pre-rotations by a RotX and
    # post-rotations by a RotY...
    return R(2*pi*rand(T), 2*pi*rand(T))
end

"""
    struct RotXY{T} <: Rotation{3,T}
    RotXY(theta_x, theta_y)

A 3×3 rotation matrix which represents a rotation by `theta_y` about the Y axis,
followed by a rotation by `theta_x` about the X axis.
"""
RotXY

@inline function Base.Tuple(r::RotXY{T}) where T
    sinθ₁, cosθ₁ = sincos(r.theta1)
    sinθ₂, cosθ₂ = sincos(r.theta2)
    z = zero(sinθ₁)

    # transposed representation
    (cosθ₂,  sinθ₁*sinθ₂,    cosθ₁*-sinθ₂,
     z,      cosθ₁,          sinθ₁,
     sinθ₂,  -sinθ₁*cosθ₂,   cosθ₁*cosθ₂)
end

@inline function Base.:*(r::RotXY, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    sinθ₁, cosθ₁ = sincos(r.theta1)
    sinθ₂, cosθ₂ = sincos(r.theta2)

    T = Base.promote_op(*, typeof(sinθ₁), eltype(v))

    return similar_type(v,T)(cosθ₂*v[1] + sinθ₂*v[3],
                             sinθ₁*sinθ₂*v[1] + cosθ₁*v[2] + -sinθ₁*cosθ₂*v[3],
                             cosθ₁*-sinθ₂*v[1] + sinθ₁*v[2] + cosθ₁*cosθ₂*v[3])
end

"""
    struct RotYX{T} <: Rotation{3,T}
    RotYX(theta_y, theta_x)

A 3×3 rotation matrix which represents a rotation by `theta_x` about the X axis,
followed by a rotation by `theta_y` about the Y axis.
"""
RotYX

@inline function Base.Tuple(r::RotYX{T}) where T
    sinθ₁, cosθ₁ = sincos(r.theta1)
    sinθ₂, cosθ₂ = sincos(r.theta2)
    z = zero(sinθ₁)

    # transposed representation
    (cosθ₁,        z,       -sinθ₁,
     sinθ₁*sinθ₂,  cosθ₂,   cosθ₁*sinθ₂,
     sinθ₁*cosθ₂,  -sinθ₂,  cosθ₁*cosθ₂)
end

@inline function Base.:*(r::RotYX, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    sinθ₁, cosθ₁ = sincos(r.theta1)
    sinθ₂, cosθ₂ = sincos(r.theta2)

    T = Base.promote_op(*, typeof(sinθ₁), eltype(v))

    return similar_type(v,T)(cosθ₁*v[1] + sinθ₁*sinθ₂*v[2] + sinθ₁*cosθ₂*v[3],
                             cosθ₂*v[2] + -sinθ₂*v[3],
                             -sinθ₁*v[1] + cosθ₁*sinθ₂*v[2] + cosθ₁*cosθ₂*v[3])
end

"""
    struct RotXZ{T} <: Rotation{3,T}
    RotXZ(theta_x, theta_z)

A 3×3 rotation matrix which represents a rotation by `theta_z` about the Z axis,
followed by a rotation by `theta_x` about the X axis.
"""
RotXZ

@inline function Base.Tuple(r::RotXZ{T}) where T
    sinθ₁, cosθ₁ = sincos(r.theta1)
    sinθ₂, cosθ₂ = sincos(r.theta2)
    z = zero(sinθ₁)

    # transposed representation
    (cosθ₂,   cosθ₁*sinθ₂,  sinθ₁*sinθ₂,
     -sinθ₂,  cosθ₁*cosθ₂,  sinθ₁*cosθ₂,
     z,       -sinθ₁,       cosθ₁)
end

@inline function Base.:*(r::RotXZ, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    sinθ₁, cosθ₁ = sincos(r.theta1)
    sinθ₂, cosθ₂ = sincos(r.theta2)

    T = Base.promote_op(*, typeof(sinθ₁), eltype(v))

    return similar_type(v,T)(cosθ₂*v[1] + -sinθ₂*v[2],
                             cosθ₁*sinθ₂*v[1] + cosθ₁*cosθ₂*v[2] + -sinθ₁*v[3],
                             sinθ₁*sinθ₂*v[1] + sinθ₁*cosθ₂*v[2] + cosθ₁*v[3])
end

"""
    struct RotZX{T} <: Rotation{3,T}
    RotZX(theta_z, theta_x)

A 3×3 rotation matrix which represents a rotation by `theta_x` about the X axis,
followed by a rotation by `theta_z` about the Z axis.
"""
RotZX

@inline function Base.Tuple(r::RotZX{T}) where T
    sinθ₁, cosθ₁ = sincos(r.theta1)
    sinθ₂, cosθ₂ = sincos(r.theta2)
    z = zero(sinθ₁)

    # transposed representation
    ( cosθ₁,         sinθ₁,         z,
     -sinθ₁*cosθ₂,   cosθ₁*cosθ₂,   sinθ₂,
      sinθ₁*sinθ₂,   cosθ₁*-sinθ₂,  cosθ₂)
end

@inline function Base.:*(r::RotZX, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    sinθ₁, cosθ₁ = sincos(r.theta1)
    sinθ₂, cosθ₂ = sincos(r.theta2)

    T = Base.promote_op(*, typeof(sinθ₁), eltype(v))

    return similar_type(v,T)(cosθ₁*v[1] + -sinθ₁*cosθ₂*v[2] + sinθ₁*sinθ₂*v[3],
                             sinθ₁*v[1] + cosθ₁*cosθ₂*v[2] + cosθ₁*-sinθ₂*v[3],
                             sinθ₂*v[2] + cosθ₂*v[3])
end

"""
    struct RotZY{T} <: Rotation{3,T}
    RotZY(theta_z, theta_y)

A 3×3 rotation matrix which represents a rotation by `theta_y` about the Y axis,
followed by a rotation by `theta_z` about the Z axis.
"""
RotZY

@inline function Base.Tuple(r::RotZY{T}) where T
    sinθ₁, cosθ₁ = sincos(r.theta1)
    sinθ₂, cosθ₂ = sincos(r.theta2)
    z = zero(sinθ₁)

    # transposed representation
    ( cosθ₁*cosθ₂,  sinθ₁*cosθ₂, -sinθ₂,
     -sinθ₁,        cosθ₁,        z,
      cosθ₁*sinθ₂,  sinθ₁*sinθ₂,  cosθ₂)
end

@inline function Base.:*(r::RotZY, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    sinθ₁, cosθ₁ = sincos(r.theta1)
    sinθ₂, cosθ₂ = sincos(r.theta2)

    T = Base.promote_op(*, typeof(sinθ₁), eltype(v))

    return similar_type(v,T)(cosθ₁*cosθ₂*v[1] + -sinθ₁*v[2] + cosθ₁*sinθ₂*v[3],
                             sinθ₁*cosθ₂*v[1] + cosθ₁*v[2] + sinθ₁*sinθ₂*v[3],
                             -sinθ₂*v[1] + cosθ₂*v[3])
end

"""
    struct RotYZ{T} <: Rotation{3,T}
    RotYZ(theta_y, theta_z)

A 3×3 rotation matrix which represents a rotation by `theta_z` about the Z axis,
followed by a rotation by `theta_y` about the Y axis.
"""
RotYZ

@inline function Base.Tuple(r::RotYZ{T}) where T
    sinθ₁, cosθ₁ = sincos(r.theta1)
    sinθ₂, cosθ₂ = sincos(r.theta2)
    z = zero(sinθ₁)

    # transposed representation
    (cosθ₁*cosθ₂,   sinθ₂,    -sinθ₁*cosθ₂,
     cosθ₁*-sinθ₂,  cosθ₂,     sinθ₁*sinθ₂,
     sinθ₁,         z,         cosθ₁)
end

@inline function Base.:*(r::RotYZ, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    sinθ₁, cosθ₁ = sincos(r.theta1)
    sinθ₂, cosθ₂ = sincos(r.theta2)

    T = Base.promote_op(*, typeof(sinθ₁), eltype(v))

    return similar_type(v,T)(cosθ₁*cosθ₂*v[1] + cosθ₁*-sinθ₂*v[2] + sinθ₁*v[3],
                             sinθ₂*v[1] + cosθ₂*v[2],
                             -sinθ₁*cosθ₂*v[1] + sinθ₁*sinθ₂*v[2] + cosθ₁*v[3])
end

################################################################################
################################################################################

########################
# Three axis Rotations #
########################

for axis1 in [:X, :Y, :Z]
    Rot1Type = Symbol("Rot" * string(axis1))
    for axis2 in filter(axis -> axis != axis1, [:X, :Y, :Z])
        Rot2Type = Symbol("Rot" * string(axis2))
        Rot12Type = Symbol("Rot" * string(axis1) * string(axis2))
        for axis3 in filter(axis -> axis != axis2, [:X, :Y, :Z])
            Rot3Type = Symbol("Rot" * string(axis3))
            Rot23Type = Symbol("Rot" * string(axis2) * string(axis3))
            RotType = Symbol("Rot" * string(axis1) * string(axis2) * string(axis3))
            InvRotType = Symbol("Rot" * string(axis3) * string(axis2) * string(axis1))

            @eval begin
                struct $RotType{T} <: Rotation{3,T}
                    theta1::T
                    theta2::T
                    theta3::T
                    $RotType{T}(theta1, theta2, theta3) where {T} = new{T}(theta1, theta2, theta3)
                    $RotType{T}(r::$RotType) where {T} = new{T}(r.theta1, r.theta2, r.theta3)
                end

                @inline $RotType(theta1::T1, theta2::T2, theta3::T3) where {T1, T2, T3} = $RotType{promote_type(promote_type(T1, T2), T3)}(theta1, theta2, theta3)
                @inline $RotType(r::$RotType{T}) where {T} = $RotType{T}(r)

                @inline function Base.getindex(r::$RotType{T}, i::Int) where T
                    Tuple(r)[i] # Slow...
                end

                # Composing single-axis rotations with two-axis rotations:
                @inline Base.:*(r1::$Rot1Type, r2::$Rot23Type) = $RotType(r1.theta, r2.theta1, r2.theta2)
                @inline Base.:*(r1::$Rot12Type, r2::$Rot3Type) = $RotType(r1.theta1, r1.theta2, r2.theta)

                # Composing with single-axis rotations:
                @inline Base.:*(r1::$RotType, r2::$Rot3Type) = $RotType(r1.theta1, r1.theta2, r1.theta3 + r2.theta)
                @inline Base.:*(r1::$Rot1Type, r2::$RotType) = $RotType(r1.theta + r2.theta1, r2.theta2, r2.theta3)

                @inline Base.inv(r::$RotType) = $InvRotType(-r.theta3, -r.theta2, -r.theta1)

                # define null rotations for convenience
                @inline Base.one(::Type{$RotType}) = $RotType(0.0, 0.0, 0.0)
                @inline Base.one(::Type{$RotType{T}}) where {T} = $RotType{T}(zero(T), zero(T), zero(T))
            end
        end
    end
end

################################################################################
################################################################################

##########################
# Proper Euler Rotations #
##########################

"""
    struct RotXYX{T} <: Rotation{3,T}
    RotXYX(theta1, theta2, theta3)

A 3×3 rotation matrix parameterized by the "proper" XYX Euler angle convention,
consisting of first a rotation about the X axis by `theta3`, followed by a
rotation about the Y axis by `theta2`, and finally a rotation about the X axis
by `theta1`.
"""
RotXYX

@inline function (::Type{Rot})(t::NTuple{9}) where Rot <: RotXYX
    R = SMatrix{3,3}(t)

    t1 = atan(R[2, 1], (-R[3, 1] + eps(t[1])) - eps(t[1]))  # TODO: handle denormal numbers better, as atan(0,0) != atan(0,-0)
    st1, ct1 = sincos(t1)

    Rot(t1,
        atan((R[1, 2] * R[1, 2] + R[1, 3] * R[1, 3])^(1/2), R[1, 1]),
        atan(- R[2, 3]*ct1 - R[3, 3]*st1, R[2, 2]*ct1 + R[3, 2]*st1))
end

@inline function Base.Tuple(r::RotXYX{T}) where T
    sinθ₁, cosθ₁ = sincos(r.theta1)
    sinθ₂, cosθ₂ = sincos(r.theta2)
    sinθ₃, cosθ₃ = sincos(r.theta3)

    # transposed representation
    (cosθ₂,        sinθ₁*sinθ₂,                        cosθ₁*-sinθ₂,
     sinθ₂*sinθ₃,  cosθ₁*cosθ₃ + -sinθ₁*cosθ₂*sinθ₃,   sinθ₁*cosθ₃ + cosθ₁*cosθ₂*sinθ₃,
     sinθ₂*cosθ₃,  cosθ₁*-sinθ₃ + -sinθ₁*cosθ₂*cosθ₃,  sinθ₁*-sinθ₃ + cosθ₁*cosθ₂*cosθ₃)
end

@inline function Base.:*(r::RotXYX, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    sinθ₁, cosθ₁ = sincos(r.theta1)
    sinθ₂, cosθ₂ = sincos(r.theta2)
    sinθ₃, cosθ₃ = sincos(r.theta3)

    T = Base.promote_op(*, typeof(sinθ₁), eltype(v))

    return similar_type(v,T)(
        cosθ₂*v[1] + sinθ₂*sinθ₃*v[2] + sinθ₂*cosθ₃*v[3],
        sinθ₁*sinθ₂*v[1] + (cosθ₁*cosθ₃ + -sinθ₁*cosθ₂*sinθ₃)*v[2] + (cosθ₁*-sinθ₃ + -sinθ₁*cosθ₂*cosθ₃)*v[3],
        cosθ₁*-sinθ₂*v[1] + (sinθ₁*cosθ₃ + cosθ₁*cosθ₂*sinθ₃)*v[2] + (sinθ₁*-sinθ₃ + cosθ₁*cosθ₂*cosθ₃)*v[3])
end

"""
    struct RotXZX{T} <: Rotation{3,T}
    RotXZX(theta1, theta2, theta3)

A 3×3 rotation matrix parameterized by the "proper" XZX Euler angle convention,
consisting of first a rotation about the X axis by `theta3`, followed by a
rotation about the Z axis by `theta2`, and finally a rotation about the X axis
by `theta1`.
"""
RotXZX

@inline function (::Type{Rot})(t::NTuple{9}) where Rot <: RotXZX
    R = SMatrix{3,3}(t)

    t1 = atan(R[3, 1], R[2, 1])
    st1, ct1 = sincos(t1)

    Rot(t1,
        atan((R[1, 2] * R[1, 2] + R[1, 3] * R[1, 3])^(1/2), R[1, 1]),
        atan(R[3, 2]*ct1 - R[2, 2]*st1, R[3, 3]*ct1 - R[2, 3]*st1))
end

@inline function Base.Tuple(r::RotXZX{T}) where T
    sinθ₁, cosθ₁ = sincos(r.theta1)
    sinθ₂, cosθ₂ = sincos(r.theta2)
    sinθ₃, cosθ₃ = sincos(r.theta3)

    # transposed representation
    (cosθ₂,         cosθ₁*sinθ₂,                        sinθ₁*sinθ₂,
     -sinθ₂*cosθ₃,  cosθ₁*cosθ₂*cosθ₃ + -sinθ₁*sinθ₃,   sinθ₁*cosθ₂*cosθ₃ + cosθ₁*sinθ₃,
     sinθ₂*sinθ₃,   cosθ₁*cosθ₂*-sinθ₃ + -sinθ₁*cosθ₃,  sinθ₁*cosθ₂*-sinθ₃ + cosθ₁*cosθ₃)
end

@inline function Base.:*(r::RotXZX, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    sinθ₁, cosθ₁ = sincos(r.theta1)
    sinθ₂, cosθ₂ = sincos(r.theta2)
    sinθ₃, cosθ₃ = sincos(r.theta3)

    T = Base.promote_op(*, typeof(sinθ₁), eltype(v))

    return similar_type(v,T)(
        cosθ₂*v[1] + -sinθ₂*cosθ₃*v[2] + sinθ₂*sinθ₃*v[3],
        cosθ₁*sinθ₂*v[1] + (cosθ₁*cosθ₂*cosθ₃ + -sinθ₁*sinθ₃)*v[2] + (cosθ₁*cosθ₂*-sinθ₃ + -sinθ₁*cosθ₃)*v[3],
        sinθ₁*sinθ₂*v[1] + (sinθ₁*cosθ₂*cosθ₃ + cosθ₁*sinθ₃)*v[2] + (sinθ₁*cosθ₂*-sinθ₃ + cosθ₁*cosθ₃)*v[3])
end

"""
    struct RotYXY{T} <: Rotation{3,T}
    RotYXY(theta1, theta2, theta3)

A 3×3 rotation matrix parameterized by the "proper" YXY Euler angle convention,
consisting of first a rotation about the Y axis by `theta3`, followed by a
rotation about the X axis by `theta2`, and finally a rotation about the Y axis
by `theta1`.
"""
RotYXY

@inline function (::Type{Rot})(t::NTuple{9}) where Rot <: RotYXY
    R = SMatrix{3,3}(t)

    t1 = atan(R[1, 2], R[3, 2])
    st1, ct1 = sincos(t1)

    Rot(t1,
        atan((R[2, 1] * R[2, 1] + R[2, 3] * R[2, 3])^(1/2), R[2, 2]),
        atan(R[1, 3]*ct1 - R[3, 3]*st1, R[1, 1]*ct1 - R[3, 1]*st1))
end

@inline function Base.Tuple(r::RotYXY{T}) where T
    sinθ₁, cosθ₁ = sincos(r.theta1)
    sinθ₂, cosθ₂ = sincos(r.theta2)
    sinθ₃, cosθ₃ = sincos(r.theta3)

    # transposed representation
    (cosθ₁*cosθ₃ + sinθ₁*cosθ₂*-sinθ₃,  sinθ₂*sinθ₃,   -sinθ₁*cosθ₃ + cosθ₁*cosθ₂*-sinθ₃,
     sinθ₁*sinθ₂,                       cosθ₂,         cosθ₁*sinθ₂,
     cosθ₁*sinθ₃ + sinθ₁*cosθ₂*cosθ₃,   -sinθ₂*cosθ₃,  -sinθ₁*sinθ₃ + cosθ₁*cosθ₂*cosθ₃)
end

@inline function Base.:*(r::RotYXY, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    sinθ₁, cosθ₁ = sincos(r.theta1)
    sinθ₂, cosθ₂ = sincos(r.theta2)
    sinθ₃, cosθ₃ = sincos(r.theta3)

    T = Base.promote_op(*, typeof(sinθ₁), eltype(v))

    return similar_type(v,T)(
        (cosθ₁*cosθ₃ + sinθ₁*cosθ₂*-sinθ₃)*v[1] + sinθ₁*sinθ₂*v[2] + (cosθ₁*sinθ₃ + sinθ₁*cosθ₂*cosθ₃)*v[3],
        sinθ₂*sinθ₃*v[1] + cosθ₂*v[2] + -sinθ₂*cosθ₃*v[3],
        (-sinθ₁*cosθ₃ + cosθ₁*cosθ₂*-sinθ₃)*v[1] + cosθ₁*sinθ₂*v[2] + (-sinθ₁*sinθ₃ + cosθ₁*cosθ₂*cosθ₃)*v[3])
end

"""
    struct RotYZY{T} <: Rotation{3,T}
    RotYZY(theta1, theta2, theta3)

A 3×3 rotation matrix parameterized by the "proper" YXY Euler angle convention,
consisting of first a rotation about the Y axis by `theta3`, followed by a
rotation about the Z axis by `theta2`, and finally a rotation about the Y axis
by `theta1`.
"""
RotYZY

@inline function (::Type{Rot})(t::NTuple{9}) where Rot <: RotYZY
    R = SMatrix{3,3}(t)

    t1 = atan(R[3, 2], -R[1, 2])  # TODO: handle denormal numbers better, as atan(0,0) != atan(0,-0)
    st1, ct1 = sincos(t1)

    Rot(t1,
        atan((R[2, 1] * R[2, 1] + R[2, 3] * R[2, 3])^(1/2), R[2, 2]),
        atan(- R[3, 1]*ct1 - R[1, 1]*st1, R[3, 3]*ct1 + R[1, 3]*st1))
end

@inline function Base.Tuple(r::RotYZY{T}) where T
    sinθ₁, cosθ₁ = sincos(r.theta1)
    sinθ₂, cosθ₂ = sincos(r.theta2)
    sinθ₃, cosθ₃ = sincos(r.theta3)

    # transposed representation
    (cosθ₁*cosθ₂*cosθ₃ + sinθ₁*-sinθ₃,  sinθ₂*cosθ₃,  -sinθ₁*cosθ₂*cosθ₃ + cosθ₁*-sinθ₃,
     cosθ₁*-sinθ₂,                      cosθ₂,        sinθ₁*sinθ₂,
     cosθ₁*cosθ₂*sinθ₃ + sinθ₁*cosθ₃,   sinθ₂*sinθ₃,  -sinθ₁*cosθ₂*sinθ₃ + cosθ₁*cosθ₃)
end

@inline function Base.:*(r::RotYZY, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    sinθ₁, cosθ₁ = sincos(r.theta1)
    sinθ₂, cosθ₂ = sincos(r.theta2)
    sinθ₃, cosθ₃ = sincos(r.theta3)

    T = Base.promote_op(*, typeof(sinθ₁), eltype(v))

    return similar_type(v,T)(
        (cosθ₁*cosθ₂*cosθ₃ + sinθ₁*-sinθ₃)*v[1] + cosθ₁*-sinθ₂*v[2] + (cosθ₁*cosθ₂*sinθ₃ + sinθ₁*cosθ₃)*v[3],
        sinθ₂*cosθ₃*v[1] + cosθ₂*v[2] + sinθ₂*sinθ₃*v[3],
        (-sinθ₁*cosθ₂*cosθ₃ + cosθ₁*-sinθ₃)*v[1] + sinθ₁*sinθ₂*v[2] + (-sinθ₁*cosθ₂*sinθ₃ + cosθ₁*cosθ₃)*v[3])
end

"""
    struct RotZXZ{T} <: Rotation{3,T}
    RotZXZ(theta1, theta2, theta3)

A 3×3 rotation matrix parameterized by the "proper" ZXZ Euler angle convention,
consisting of first a rotation about the Z axis by `theta3`, followed by a
rotation about the X axis by `theta2`, and finally a rotation about the Z axis
by `theta1`.
"""
RotZXZ

@inline function (::Type{Rot})(t::NTuple{9}) where Rot <: RotZXZ
    R = SMatrix{3,3}(t)

    t1 = atan(R[1, 3], (-R[2, 3] + eps()) - eps())  # TODO: handle denormal numbers better, as atan(0,0) != atan(0,-0)
    st1, ct1 = sincos(t1)

    Rot(t1,
        atan((R[3, 1] * R[3, 1] + R[3, 2] * R[3, 2])^(1/2), R[3, 3]),
        atan(- R[1, 2]*ct1 - R[2, 2]*st1, R[1, 1]*ct1 + R[2, 1]*st1))
end

@inline function Base.Tuple(r::RotZXZ{T}) where T
    sinθ₁, cosθ₁ = sincos(r.theta1)
    sinθ₂, cosθ₂ = sincos(r.theta2)
    sinθ₃, cosθ₃ = sincos(r.theta3)

    # transposed representation
    (cosθ₁*cosθ₃ + -sinθ₁*cosθ₂*sinθ₃,   sinθ₁*cosθ₃ + cosθ₁*cosθ₂*sinθ₃,   sinθ₂*sinθ₃,
     cosθ₁*-sinθ₃ + -sinθ₁*cosθ₂*cosθ₃,  sinθ₁*-sinθ₃ + cosθ₁*cosθ₂*cosθ₃,  sinθ₂*cosθ₃,
     sinθ₁*sinθ₂,                        cosθ₁*-sinθ₂,                      cosθ₂)
end

@inline function Base.:*(r::RotZXZ, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    sinθ₁, cosθ₁ = sincos(r.theta1)
    sinθ₂, cosθ₂ = sincos(r.theta2)
    sinθ₃, cosθ₃ = sincos(r.theta3)

    T = Base.promote_op(*, typeof(sinθ₁), eltype(v))

    return similar_type(v,T)(
         (cosθ₁*cosθ₃ + -sinθ₁*cosθ₂*sinθ₃)*v[1] + (cosθ₁*-sinθ₃ + -sinθ₁*cosθ₂*cosθ₃)*v[2] + sinθ₁*sinθ₂*v[3],
         (sinθ₁*cosθ₃ + cosθ₁*cosθ₂*sinθ₃)*v[1] + (sinθ₁*-sinθ₃ + cosθ₁*cosθ₂*cosθ₃)*v[2] + cosθ₁*-sinθ₂*v[3],
         sinθ₂*sinθ₃*v[1] + sinθ₂*cosθ₃*v[2] + cosθ₂*v[3])
end

"""
    struct RotZYZ{T} <: Rotation{3,T}
    RotZYZ(theta1, theta2, theta3)

A 3×3 rotation matrix parameterized by the "proper" ZXZ Euler angle convention,
consisting of first a rotation about the Z axis by `theta3`, followed by a
rotation about the Y axis by `theta2`, and finally a rotation about the Z axis
by `theta1`.
"""
RotZYZ

@inline function (::Type{Rot})(t::NTuple{9}) where Rot <: RotZYZ
    R = SMatrix{3,3}(t)

    t1 = atan(R[2, 3], R[1, 3])
    st1, ct1 = sincos(t1)

    Rot(t1,
        atan((R[3, 1] * R[3, 1] + R[3, 2] * R[3, 2])^(1/2), R[3, 3]),
        atan(R[2, 1]*ct1 - R[1, 1]*st1, R[2, 2]*ct1 - R[1, 2]*st1))
end

@inline function Base.Tuple(r::RotZYZ{T}) where T
    sinθ₁, cosθ₁ = sincos(r.theta1)
    sinθ₂, cosθ₂ = sincos(r.theta2)
    sinθ₃, cosθ₃ = sincos(r.theta3)

    # transposed representation
    (cosθ₁*cosθ₂*cosθ₃ + -sinθ₁*sinθ₃,   sinθ₁*cosθ₂*cosθ₃ + cosθ₁*sinθ₃,   -sinθ₂*cosθ₃,
     cosθ₁*cosθ₂*-sinθ₃ + -sinθ₁*cosθ₃,  sinθ₁*cosθ₂*-sinθ₃ + cosθ₁*cosθ₃,  sinθ₂*sinθ₃,
     cosθ₁*sinθ₂,                        sinθ₁*sinθ₂,                       cosθ₂)
end

@inline function Base.:*(r::RotZYZ, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    sinθ₁, cosθ₁ = sincos(r.theta1)
    sinθ₂, cosθ₂ = sincos(r.theta2)
    sinθ₃, cosθ₃ = sincos(r.theta3)

    T = Base.promote_op(*, typeof(sinθ₁), eltype(v))

    return similar_type(v,T)(
        (cosθ₁*cosθ₂*cosθ₃ + -sinθ₁*sinθ₃)*v[1] + (cosθ₁*cosθ₂*-sinθ₃ + -sinθ₁*cosθ₃)*v[2] + cosθ₁*sinθ₂*v[3],
        (sinθ₁*cosθ₂*cosθ₃ + cosθ₁*sinθ₃)*v[1] + (sinθ₁*cosθ₂*-sinθ₃ + cosθ₁*cosθ₃)*v[2] + sinθ₁*sinθ₂*v[3],
        -sinθ₂*cosθ₃*v[1] + sinθ₂*sinθ₃*v[2] + cosθ₂*v[3])
end

###############################
# Tait-Bryant Euler Rotations #
###############################

"""
    struct RotXYZ{T} <: Rotation{3,T}
    RotXYZ(theta1, theta2, theta3)
    RotXYZ(roll=r, pitch=p, yaw=y)

A 3×3 rotation matrix parameterized by the "Tait-Bryant" XYZ Euler angle
convention, consisting of first a rotation about the Z axis by `theta3`,
followed by a rotation about the Y axis by `theta2`, and finally a rotation
about the X axis by `theta1`.

The keyword argument form applies roll, pitch and yaw to the X, Y and Z axes
respectively, in XYZ order. (Because it is a right-handed coordinate system,
note that positive pitch is heading in the negative Z axis).
"""
RotXYZ

@inline (::Type{Rot})(; roll=0, pitch=0, yaw=0) where {Rot<:RotXYZ} = Rot(roll, pitch, yaw)

@inline function (::Type{Rot})(t::NTuple{9}) where Rot <: RotXYZ
    R = SMatrix{3,3}(t)

    t1 = atan(-R[2, 3], R[3, 3])
    st1, ct1 = sincos(t1)

    Rot(t1,
        atan(R[1, 3], (R[1, 1] * R[1, 1] + R[1, 2] * R[1, 2])^(1/2)),
        atan(R[2, 1]*ct1 + R[3, 1]*st1, R[2, 2]*ct1 + R[3, 2]*st1))
end

@inline function Base.Tuple(r::RotXYZ{T}) where T
    sinθ₁, cosθ₁ = sincos(r.theta1)
    sinθ₂, cosθ₂ = sincos(r.theta2)
    sinθ₃, cosθ₃ = sincos(r.theta3)

    # transposed representation
    (cosθ₂*cosθ₃,   sinθ₁*sinθ₂*cosθ₃ + cosθ₁*sinθ₃,   cosθ₁*-sinθ₂*cosθ₃ + sinθ₁*sinθ₃,
     cosθ₂*-sinθ₃,  sinθ₁*sinθ₂*-sinθ₃ + cosθ₁*cosθ₃,  cosθ₁*-sinθ₂*-sinθ₃ + sinθ₁*cosθ₃,
     sinθ₂,         -sinθ₁*cosθ₂,                      cosθ₁*cosθ₂)
end

@inline function Base.:*(r::RotXYZ, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    sinθ₁, cosθ₁ = sincos(r.theta1)
    sinθ₂, cosθ₂ = sincos(r.theta2)
    sinθ₃, cosθ₃ = sincos(r.theta3)

    T = Base.promote_op(*, typeof(sinθ₁), eltype(v))

    return similar_type(v,T)(
        cosθ₂*cosθ₃*v[1] + cosθ₂*-sinθ₃*v[2] + sinθ₂*v[3],
        (sinθ₁*sinθ₂*cosθ₃ + cosθ₁*sinθ₃)*v[1] + (sinθ₁*sinθ₂*-sinθ₃ + cosθ₁*cosθ₃)*v[2] + -sinθ₁*cosθ₂*v[3],
        (cosθ₁*-sinθ₂*cosθ₃ + sinθ₁*sinθ₃)*v[1] + (cosθ₁*sinθ₂*sinθ₃ + sinθ₁*cosθ₃)*v[2] + cosθ₁*cosθ₂*v[3])
end

"""
    struct RotZYX{T} <: Rotation{3,T}
    RotZYX(theta1, theta2, theta3)
    RotZYX(roll=r, pitch=p, yaw=y)

A 3×3 rotation matrix parameterized by the "Tait-Bryant" ZYX Euler angle
convention, consisting of first a rotation about the X axis by `theta3`,
followed by a rotation about the Y axis by `theta2`, and finally a rotation
about the Z axis by `theta1`.

The keyword argument form applies roll, pitch and yaw to the X, Y and Z axes
respectively, in ZYX order. (Because it is a right-handed coordinate system,
note that positive pitch is heading in the negative Z axis).
"""
RotZYX

@inline (::Type{Rot})(; roll=0, pitch=0, yaw=0) where {Rot<:RotZYX} = Rot(yaw, pitch, roll)

@inline function (::Type{Rot})(t::NTuple{9}) where Rot <: RotZYX
    R = SMatrix{3,3}(t)

    t1 = atan(R[2, 1], R[1, 1])
    st1, ct1 = sincos(t1)

    Rot(t1,
        atan(-R[3, 1], (R[3, 2] * R[3, 2] + R[3, 3] * R[3, 3])^(1/2)),
        atan(R[1, 3]*st1 - R[2, 3]*ct1, R[2, 2]*ct1 - R[1, 2]*st1))
end

@inline function Base.Tuple(r::RotZYX{T}) where T
    sinθ₁, cosθ₁ = sincos(r.theta1)
    sinθ₂, cosθ₂ = sincos(r.theta2)
    sinθ₃, cosθ₃ = sincos(r.theta3)

    # transposed representation
    ( cosθ₁*cosθ₂,                       sinθ₁*cosθ₂,                      -sinθ₂,
     -sinθ₁*cosθ₃ + cosθ₁*sinθ₂*sinθ₃,   cosθ₁*cosθ₃ + sinθ₁*sinθ₂*sinθ₃,   cosθ₂*sinθ₃,
      sinθ₁*sinθ₃ + cosθ₁*sinθ₂*cosθ₃,   cosθ₁*-sinθ₃ + sinθ₁*sinθ₂*cosθ₃,  cosθ₂*cosθ₃)
end

@inline function Base.:*(r::RotZYX, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    sinθ₁, cosθ₁ = sincos(r.theta1)
    sinθ₂, cosθ₂ = sincos(r.theta2)
    sinθ₃, cosθ₃ = sincos(r.theta3)

    T = Base.promote_op(*, typeof(sinθ₁), eltype(v))

    return similar_type(v,T)(
        cosθ₁*cosθ₂*v[1] + (-sinθ₁*cosθ₃ + cosθ₁*sinθ₂*sinθ₃)*v[2] + (sinθ₁*sinθ₃ + cosθ₁*sinθ₂*cosθ₃)*v[3],
        sinθ₁*cosθ₂*v[1] + (cosθ₁*cosθ₃ + sinθ₁*sinθ₂*sinθ₃)*v[2] + (cosθ₁*-sinθ₃ + sinθ₁*sinθ₂*cosθ₃)*v[3],
        -sinθ₂*v[1] + cosθ₂*sinθ₃*v[2] + cosθ₂*cosθ₃*v[3])
end

"""
    struct RotXZY{T} <: Rotation{3,T}
    RotXZY(theta1, theta2, theta3)
    RotXZY(roll=r, pitch=p, yaw=y)

A 3×3 rotation matrix parameterized by the "Tait-Bryant" XZY Euler angle
convention, consisting of first a rotation about the Y axis by `theta3`,
followed by a rotation about the Z axis by `theta2`, and finally a rotation
about the X axis by `theta1`.

The keyword argument form applies roll, pitch and yaw to the X, Y and Z axes
respectively, in XZY order. (Because it is a right-handed coordinate system,
note that positive pitch is heading in the negative Z axis).
"""
RotXZY

@inline (::Type{Rot})(; roll=0, pitch=0, yaw=0) where {Rot<:RotXZY} = Rot(roll, yaw, pitch)

@inline function (::Type{Rot})(t::NTuple{9}) where Rot <: RotXZY
    R = SMatrix{3,3}(t)

    t1 = atan(R[3, 2], R[2, 2])
    st1, ct1 = sincos(t1)

    Rot(t1,
        atan(-R[1, 2], (R[1, 1] * R[1, 1] + R[1, 3] * R[1, 3])^(1/2)),
        atan(R[2, 1]*st1 - R[3, 1]*ct1, R[3, 3]*ct1 - R[2, 3]*st1))
end

@inline function Base.Tuple(r::RotXZY{T}) where T
    sinθ₁, cosθ₁ = sincos(r.theta1)
    sinθ₂, cosθ₂ = sincos(r.theta2)
    sinθ₃, cosθ₃ = sincos(r.theta3)

    # transposed representation
    ( cosθ₂*cosθ₃,  cosθ₁*sinθ₂*cosθ₃ + sinθ₁*sinθ₃,   sinθ₁*sinθ₂*cosθ₃ + cosθ₁*-sinθ₃,
     -sinθ₂,        cosθ₁*cosθ₂,                       sinθ₁*cosθ₂,
      cosθ₂*sinθ₃,  cosθ₁*sinθ₂*sinθ₃ + -sinθ₁*cosθ₃,  sinθ₁*sinθ₂*sinθ₃ + cosθ₁*cosθ₃)
end

@inline function Base.:*(r::RotXZY, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    sinθ₁, cosθ₁ = sincos(r.theta1)
    sinθ₂, cosθ₂ = sincos(r.theta2)
    sinθ₃, cosθ₃ = sincos(r.theta3)

    T = Base.promote_op(*, typeof(sinθ₁), eltype(v))

    return similar_type(v,T)(
        cosθ₂*cosθ₃*v[1] + -sinθ₂*v[2] + cosθ₂*sinθ₃*v[3],
        (cosθ₁*sinθ₂*cosθ₃ + sinθ₁*sinθ₃)*v[1] + cosθ₁*cosθ₂*v[2] + (cosθ₁*sinθ₂*sinθ₃ + -sinθ₁*cosθ₃)*v[3],
        (sinθ₁*sinθ₂*cosθ₃ + cosθ₁*-sinθ₃)*v[1] + sinθ₁*cosθ₂*v[2] + (sinθ₁*sinθ₂*sinθ₃ + cosθ₁*cosθ₃)*v[3])
end

"""
    struct RotYZX{T} <: Rotation{3,T}
    RotYZX(theta1, theta2, theta3)
    RotYZX(roll=r, pitch=p, yaw=y)

A 3×3 rotation matrix parameterized by the "Tait-Bryant" YZX Euler angle
convention, consisting of first a rotation about the X axis by `theta3`,
followed by a rotation about the Z axis by `theta2`, and finally a rotation
about the Y axis by `theta1`.

The keyword argument form applies roll, pitch and yaw to the X, Y and Z axes
respectively, in YZX order. (Because it is a right-handed coordinate system,
note that positive pitch is heading in the negative Z axis).
"""
RotYZX

@inline (::Type{Rot})(; roll=0, pitch=0, yaw=0) where {Rot<:RotYZX} = Rot(pitch, yaw, roll)

@inline function (::Type{Rot})(t::NTuple{9}) where Rot <: RotYZX
    R = SMatrix{3,3}(t)

    t1 = atan(-R[3, 1], R[1, 1])
    st1, ct1 = sincos(t1)

    Rot(t1,
        atan(R[2, 1], (R[2, 2] * R[2, 2] + R[2, 3] * R[2, 3])^(1/2)),
        atan(R[3, 2]*ct1 + R[1, 2]*st1, R[3, 3]*ct1 + R[1, 3]*st1))
end

@inline function Base.Tuple(r::RotYZX{T}) where T
    sinθ₁, cosθ₁ = sincos(r.theta1)
    sinθ₂, cosθ₂ = sincos(r.theta2)
    sinθ₃, cosθ₃ = sincos(r.theta3)

    # transposed representation
    (cosθ₁*cosθ₂,                        sinθ₂,         -sinθ₁*cosθ₂,
     cosθ₁*-sinθ₂*cosθ₃ + sinθ₁*sinθ₃,   cosθ₂*cosθ₃,   sinθ₁*sinθ₂*cosθ₃ + cosθ₁*sinθ₃,
     cosθ₁*sinθ₂*sinθ₃ + sinθ₁*cosθ₃,    cosθ₂*-sinθ₃,  sinθ₁*sinθ₂*-sinθ₃ + cosθ₁*cosθ₃)
end

@inline function Base.:*(r::RotYZX, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    sinθ₁, cosθ₁ = sincos(r.theta1)
    sinθ₂, cosθ₂ = sincos(r.theta2)
    sinθ₃, cosθ₃ = sincos(r.theta3)

    T = Base.promote_op(*, typeof(sinθ₁), eltype(v))

    return similar_type(v,T)(
        cosθ₁*cosθ₂*v[1] + (cosθ₁*-sinθ₂*cosθ₃ + sinθ₁*sinθ₃)*v[2] + (cosθ₁*-sinθ₂*-sinθ₃ + sinθ₁*cosθ₃)*v[3],
        sinθ₂*v[1] + cosθ₂*cosθ₃*v[2] + cosθ₂*-sinθ₃*v[3],
        -sinθ₁*cosθ₂*v[1] + (sinθ₁*sinθ₂*cosθ₃ + cosθ₁*sinθ₃)*v[2] + (sinθ₁*sinθ₂*-sinθ₃ + cosθ₁*cosθ₃)*v[3])
end

"""
    struct RotYXZ{T} <: Rotation{3,T}
    RotYXZ(theta1, theta2, theta3)
    RotYXZ(roll=r, pitch=p, yaw=y)

A 3×3 rotation matrix parameterized by the "Tait-Bryant" YXZ Euler angle
convention, consisting of first a rotation about the Z axis by `theta3`,
followed by a rotation about the X axis by `theta2`, and finally a rotation
about the Y axis by `theta1`.

The keyword argument form applies roll, pitch and yaw to the X, Y and Z axes
respectively, in YXZ order. (Because it is a right-handed coordinate system,
note that positive pitch is heading in the negative Z axis).
"""
RotYXZ

@inline (::Type{Rot})(; roll=0, pitch=0, yaw=0) where {Rot<:RotYXZ} = Rot(pitch, roll, yaw)

@inline function (::Type{Rot})(t::NTuple{9}) where Rot <: RotYXZ
    R = SMatrix{3,3}(t)

    t1 = atan(R[1, 3], R[3, 3])
    st1, ct1 = sincos(t1)

    Rot(t1,
        atan(-R[2, 3], (R[2, 1] * R[2, 1] + R[2, 2] * R[2, 2])^(1/2)),
        atan(R[3, 2]*st1 - R[1, 2]*ct1, R[1, 1]*ct1 - R[3, 1]*st1))
end

@inline function Base.Tuple(r::RotYXZ{T}) where T
    sinθ₁, cosθ₁ = sincos(r.theta1)
    sinθ₂, cosθ₂ = sincos(r.theta2)
    sinθ₃, cosθ₃ = sincos(r.theta3)

    # transposed representation
    (cosθ₁*cosθ₃ + sinθ₁*sinθ₂*sinθ₃,   cosθ₂*sinθ₃,  -sinθ₁*cosθ₃ + cosθ₁*sinθ₂*sinθ₃,
     cosθ₁*-sinθ₃ + sinθ₁*sinθ₂*cosθ₃,  cosθ₂*cosθ₃,  sinθ₁*sinθ₃ + cosθ₁*sinθ₂*cosθ₃,
     sinθ₁*cosθ₂,                       -sinθ₂,       cosθ₁*cosθ₂)
end

@inline function Base.:*(r::RotYXZ, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    sinθ₁, cosθ₁ = sincos(r.theta1)
    sinθ₂, cosθ₂ = sincos(r.theta2)
    sinθ₃, cosθ₃ = sincos(r.theta3)

    T = Base.promote_op(*, typeof(sinθ₁), eltype(v))

    return similar_type(v,T)(
        (cosθ₁*cosθ₃ + sinθ₁*sinθ₂*sinθ₃)*v[1] + (cosθ₁*-sinθ₃ + sinθ₁*sinθ₂*cosθ₃)*v[2] + sinθ₁*cosθ₂*v[3],
        cosθ₂*sinθ₃*v[1] + cosθ₂*cosθ₃*v[2] + -sinθ₂*v[3],
        (-sinθ₁*cosθ₃ + cosθ₁*sinθ₂*sinθ₃)*v[1] + (sinθ₁*sinθ₃ + cosθ₁*sinθ₂*cosθ₃)*v[2] + cosθ₁*cosθ₂*v[3])
end

"""
    struct RotZXY{T} <: Rotation{3,T}
    RotZXY(theta1, theta2, theta3)
    RotZXY(roll=r, pitch=p, yaw=y)

A 3×3 rotation matrix parameterized by the "Tait-Bryant" ZXY Euler angle
convention, consisting of first a rotation about the Y axis by `theta3`,
followed by a rotation about the X axis by `theta2`, and finally a rotation
about the Z axis by `theta1`.

The keyword argument form applies roll, pitch and yaw to the X, Y and Z axes
respectively, in ZXY order. (Because it is a right-handed coordinate system,
note that positive pitch is heading in the negative Z axis).
"""
RotZXY

@inline (::Type{Rot})(; roll=0, pitch=0, yaw=0) where {Rot<:RotZXY} = Rot(yaw, roll, pitch)

@inline function (::Type{Rot})(t::NTuple{9}) where Rot <: RotZXY
    R = SMatrix{3,3}(t)

    t1 = atan(-R[1, 2], R[2, 2])
    st1, ct1 = sincos(t1)

    Rot(t1,
        atan(R[3, 2], (R[3, 1] * R[3, 1] + R[3, 3] * R[3, 3])^(1/2)),
        atan(R[1, 3]*ct1 + R[2, 3]*st1, R[1, 1]*ct1 + R[2, 1]*st1))
end

@inline function Base.Tuple(r::RotZXY{T}) where T
    sinθ₁, cosθ₁ = sincos(r.theta1)
    sinθ₂, cosθ₂ = sincos(r.theta2)
    sinθ₃, cosθ₃ = sincos(r.theta3)

    # transposed representation
    ( cosθ₁*cosθ₃ + sinθ₁*sinθ₂*-sinθ₃,  sinθ₁*cosθ₃ + cosθ₁*-sinθ₂*-sinθ₃,  cosθ₂*-sinθ₃,
     -sinθ₁*cosθ₂,                       cosθ₁*cosθ₂,                        sinθ₂,
      cosθ₁*sinθ₃ + sinθ₁*sinθ₂*cosθ₃,   sinθ₁*sinθ₃ + cosθ₁*-sinθ₂*cosθ₃,   cosθ₂*cosθ₃)
end

@inline function Base.:*(r::RotZXY, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    sinθ₁, cosθ₁ = sincos(r.theta1)
    sinθ₂, cosθ₂ = sincos(r.theta2)
    sinθ₃, cosθ₃ = sincos(r.theta3)

    T = Base.promote_op(*, typeof(sinθ₁), eltype(v))

    return similar_type(v,T)(
        (cosθ₁*cosθ₃ + sinθ₁*sinθ₂*-sinθ₃)*v[1] + -sinθ₁*cosθ₂*v[2] + (cosθ₁*sinθ₃ + sinθ₁*sinθ₂*cosθ₃)*v[3],
        (sinθ₁*cosθ₃ + cosθ₁*-sinθ₂*-sinθ₃)*v[1] + cosθ₁*cosθ₂*v[2] + (sinθ₁*sinθ₃ + cosθ₁*-sinθ₂*cosθ₃)*v[3],
         cosθ₂*-sinθ₃*v[1] + sinθ₂*v[2] + cosθ₂*cosθ₃*v[3])
end


text:
using Base.Test
using StochasticProcesses

let p = BrownianBridge(), t=linspace(0, 1, 10)
  cumsim(p, t, 3)
end

ext: jl


text:
"""
Generators struct

"""

using Parameters

include("constants.jl")

mutable struct SysGenerator
	unit_id::Array{String}
	status_init::Array{Int64,1}
	power_init::Array{Float64,1}
	uptime_init::Array{Int64,1}
	downtime_init::Array{Int64,1}
	hr_b1::Array{Float64,1}
	hr_b2::Array{Float64,1}
	hr_b3::Array{Float64,1}
	hr_b4::Array{Float64,1}
	hr_b5::Array{Float64,1}
	hr_b6::Array{Float64,1}
	hr_b7::Array{Float64,1}
	hr_b1_qty::Array{Float64,1}
	hr_b2_qty::Array{Float64,1}
	hr_b3_qty::Array{Float64,1}
	hr_b4_qty::Array{Float64,1}
	hr_b5_qty::Array{Float64,1}
	hr_b6_qty::Array{Float64,1}
	hr_b7_qty::Array{Float64,1}
	noload_hr::Array{Float64,1}
	shutdown_cost::Array{Int64,1}
	hot_startup_fixcost::Array{Float64,1}
	hot_startup_hr::Array{Float64,1}
	down_time_min::Array{Int64,1}
	up_time_min::Array{Int64,1}
	power_maxout::Array{Int64,1}
	power_minout::Array{Int64,1}
	spin_res_limit::Array{Int64,1}
	ramp_up_limit::Array{Int64,1}
	ramp_shutdn_limit::Array{Int64,1}
	isnuclear::Array{Int64,1}
	iscogen::Array{Int64,1}
	isNG_fired::Array{Int64,1}
end

mutable struct Storage
	name::Array{String}
	power::Array{Int64,1}
	pwr_ratio_min::Array{Float64,1}
	selfdischarge::Array{Int64,1}
	trip_effic_up::Array{Float64,1}
	trip_effic_dn::Array{Float64,1}
	pwr_toenergy_ratio::Array{Float64,1}
	soc_init::Array{Float64,1}
end

#=
@with_kw struct A
           a::Int = 6
           b::Float64 = -1.1
           c::UInt8
       end
=#
"""

"""
# change power_out by genout
@with_kw mutable struct GensResults
	onoff_init::Array{Int64,1} = zeros(GENS)
	power_out_init::Array{Float64,1} = zeros(GENS)
	uptime_init::Array{Float64,1} = zeros(GENS)
	dntime_init::Array{Float64,1} = zeros(GENS)
	onoff::Array{Int64,2} = zeros(GENS, INIT_HR_SUCR+INIT_HR_FUCR)
	power_out::Array{Float64,2} = zeros(GENS, INIT_HR_SUCR+INIT_HR_FUCR)
	startup::Array{Int64,2} = zeros(GENS, INIT_HR_SUCR+INIT_HR_FUCR)
	shutdown::Array{Int64,2} = zeros(GENS, INIT_HR_SUCR+INIT_HR_FUCR)
	genout_block::Array{Float64,3} = zeros(GENS, BLOCKS, INIT_HR_SUCR+INIT_HR_FUCR)

end

# in BUCRs
#TODO: Include these variables in UC_Results structures

@with_kw mutable struct PeakersResults
	onoff_init::Array{Int64,1} = zeros(PEAKERS)
	power_out_init::Array{Float64,1} = zeros(PEAKERS)
	uptime_init::Array{Float64,1} = zeros(PEAKERS)
	dntime_init::Array{Float64,1} = zeros(PEAKERS)
	onoff::Array{Int64,2} = zeros(PEAKERS, INIT_HR_SUCR+INIT_HR_FUCR)
	power_out::Array{Float64,2} = zeros(PEAKERS, INIT_HR_SUCR+INIT_HR_FUCR)
	startup::Array{Int64,2} = zeros(PEAKERS, INIT_HR_SUCR+INIT_HR_FUCR)
	shutdown::Array{Int64,2} = zeros(PEAKERS, INIT_HR_SUCR+INIT_HR_FUCR)
	genout_block::Array{Float64,3} = zeros(PEAKERS,BLOCKS,INIT_HR_SUCR-INIT_HR_FUCR)
end

#=
BUCR1_commit.gens_lb

BUCR1_commit.gens_lb::Array{Int64,1}   <- BUCR1_Commit_LB
BUCR1_commit.gens_ub::Array{Int64,1}   <- BUCR1_Commit_UB

Peakers
BUCR1_commit.peakers_lb::Array{Int64,1}  <-  BUCR1_Commit_Peaker_LB
BUCR1_commit.peakers_ub::Array{Int64,1}  <-  BUCR1_Commit_Peaker_UB
=#

@with_kw mutable struct StorageResults
	soc_init::Array{Float64,1} = zeros(STORG_UNITS)
	chrg::Array{Int64,2} = zeros(STORG_UNITS, 24-INIT_HR_SUCR+INIT_HR_FUCR)
	disc::Array{Int64,2} = zeros(STORG_UNITS, 24-INIT_HR_SUCR+INIT_HR_FUCR)
	idle::Array{Int64,2} = zeros(STORG_UNITS, 24-INIT_HR_SUCR+INIT_HR_FUCR)
	chrgpwr::Array{Float64,2} = zeros(STORG_UNITS, 24-INIT_HR_SUCR+INIT_HR_FUCR)
	discpwr::Array{Float64,2} = zeros(STORG_UNITS, 24-INIT_HR_SUCR+INIT_HR_FUCR)
	soc::Array{Float64,2} = zeros(STORG_UNITS, 24-INIT_HR_SUCR+INIT_HR_FUCR)
end;

#=
mutable struct UnitsResults
	onoff_init::Array{Int64,1}
	power_out_init::Array{Float64,1}
	uptime_init::Array{Float64,1}
	dntime_init::Array{Float64,1}
end
=#
#=
struct UC_Results
onoff_init::Array{Int64,1} 			<- FUCR_Init_genOnOff
power_out_init::Array{Float64,1}    <- FUCR_Init_genOut
uptime_init::Array{Float64,1}       <- FUCR_Init_UpTime
downtime_init::Array{Float64,1}     <- FUCR_Init_DownTime

startup  					<- FUCRtoBUCR1_genStartUp::Array{Int64,1}
shutdown 					<- FUCRtoBUCR1_genShutDown::Array{Int64,1}
onoff::Array{Int64,2} 		<- FUCRtoBUCR1_genOnOff[g,h]
power_out::Array{Float64,2} <- FUCRtoBUCR1_genOut[g,h]
genout_block 				<- FUCRtoBUCR1_genOut_Block[g,b,h]

StorageResults
soc_init::Array{Float64,1}          <- FUCR_Init_storgSOC

Still defining
for p=1:STORG_UNITS
	global FUCRtoBUCR1_storgChrg[p,h]=JuMP.value.(FUCR_storgChrg[p,h]);
	global FUCRtoBUCR1_storgDisc[p,h]=JuMP.value.(FUCR_storgDisc[p,h]);
	global FUCRtoBUCR1_storgIdle[p,h]=JuMP.value.(FUCR_storgIdle[p,h]);
	global FUCRtoBUCR1_storgChrgPwr[p,h]=JuMP.value.(FUCR_storgChrgPwr[p,h]);
	global FUCRtoBUCR1_storgDiscPwr[p,h]=JuMP.value.(FUCR_storgDiscPwr[p,h]);
	global FUCRtoBUCR1_storgSOC[p,h]=JuMP.value.(FUCR_storgSOC[p,h]);

global FUCRtoBUCR1_genOnOff[g,h]= JuMP.value.(FUCR_genOnOff[g,h]);
global FUCRtoBUCR1_genOut[g,h]= JuMP.value.(FUCR_genOut[g,h]);
global FUCRtoBUCR1_genStartUp[g,h]= JuMP.value.(FUCR_genStartUp[g,h]);
global FUCRtoBUCR1_genShutDown[g,h]=JuMP.value.(FUCR_genShutDown[g,h]);
for b=1:BLOCKS
	FUCRtoBUCR1_genOut_Block[g,b,h]=JuMP.value.(FUCR_genOut_Block[g,b,h]);
=#

mutable struct UC_Results
	gens::GensResults
	peakers::PeakersResults
	storg::StorageResults
end;

#= UC_Results
FUCRtoBUCR1.gens.onoff 			<- FUCRtoBUCR1_genOnOff
FUCRtoBUCR1.gens.power_out 		<- FUCRtoBUCR1_genOut
FUCRtoBUCR1.gens.startup 		<- FUCRtoBUCR1_genStartUp
FUCRtoBUCR1.gens.shutdown		<- FUCRtoBUCR1_genShutDown
FUCRtoBUCR1.gens.genout_block   <-  FUCRtoBUCR1_genOut_Block

FUCRtoBUCR1.peakers.onoff 			<- FUCRtoBUCR1_peakerOnOff
FUCRtoBUCR1.peakers.power_out 		<- FUCRtoBUCR1_peakerOut
FUCRtoBUCR1.peakers.startup 		<- FUCRtoBUCR1_peakerStartUp
FUCRtoBUCR1.peakers.shutdown		<- FUCRtoBUCR1_peakerShutDown
FUCRtoBUCR1.peakers.genout_block   <-  FUCRtoBUCR1_peakerOut_Block

FUCRtoBUCR1.storg.chrg				<- FUCRtoBUCR1_storgChrg
FUCRtoBUCR1.storg.disc				<- FUCRtoBUCR1_storgDisc
FUCRtoBUCR1.storg.idle				<- FUCRtoBUCR1_storgIdle
FUCRtoBUCR1.storg.chrgpwr			<- FUCRtoBUCR1_storgChrgPwr
FUCRtoBUCR1.storg.discpwr			<- FUCRtoBUCR1_storgDiscPwr
FUCRtoBUCR1.storg.soc				<- FUCRtoBUCR1_storgSOC
=#

"""
Demand Data Pre-Processing for FUCR and SUCR models
 wk_ahead: week-ahead demand data for the first UC run at 6 am
 solarg_wa: week-ahead SolarG data for the first UC run at 6 am
 windg_wa: week-ahead WindG data for the first UC run at 6 am
 hydrog_wa: week-ahead HydroG data for the first UC run at 6 am
 nuclearg_wa: week-ahead WindG data for the first UC run at 6 am
 cogeng_wa: week-ahead HydroG data for the first UC run at 6 am
"""
#=
#TODO add these variables to the DemandWAPreprocGens structure:
lb_MUT::Array{Int64,2} = zeros(Int64, GENS, HRS_FUCR);
lb_MDT::Array{Int64,2} = zeros(Int64, GENS, HRS_FUCR);
lb_MUT_Peaker::Array{Int64,2} =zeros(Int64, PEAKERS, HRS_FUCR)
lb_MDT_Peaker::Array{Int64,2} =zeros(Int64, PEAKERS, HRS_FUCR)
=#
mutable struct DemandWAPreprocGens
	wk_ahead::Array{Int64,2}
	solar_wa::Array{Int64,2}
	wind_wa::Array{Int64,2}
	hydro_wa::Array{Int64,2}
	nuclear_wa::Array{Int64,2}
	cogen_wa::Array{Int64,2}
end;

"""
Demand Data Pre-Processing for BUCR models
"""
#=
#TODO add these variables to the DemandHrPreprocGens structure:
	gens_lb::Array{Int64,1} = zeros(GENS)
	gens_ub::Array{Int64,1} = zeros(GENS)
	peakers_lb::Array{Int64,1} = zeros(PEAKERS)
	peakers_ub::Array{Int64,1} = zeros(PEAKERS)
=#
mutable struct DemandHrPreprocGens
	demand::Array{Int64,1}
	solar::Array{Int64,1}
	wind::Array{Int64,1}
	hydro::Array{Int64,1}
	nuclear::Array{Int64,1}
	cogen::Array{Int64,1}
end;

#=
bucr1_prepoc.demand::Array{Int64,2}  <- BUCR1_Hr_Demand
bucr1_prepoc.solar::Array{Int64,2}   <- BUCR1_Hr_SolarG
bucr1_prepoc.wind::Array{Int64,2}    <- BUCR1_Hr_WindG
bucr1_prepoc.hydro::Array{Int64,2}   <- BUCR1_Hr_HydroG
bucr1_prepoc.nuclear::Array{Int64,2} <- BUCR1_Hr_NuclearG
bucr1_prepoc.cogen::Array{Int64,2}   <- BUCR1_Hr_CogenG
=#

@with_kw mutable struct StorageResults
	soc_init::Array{Float64,1} = zeros(STORG_UNITS)
	chrg::Array{Int64,2} = zeros(STORG_UNITS, 24-INIT_HR_SUCR+INIT_HR_FUCR)
	disc::Array{Int64,2} = zeros(STORG_UNITS, 24-INIT_HR_SUCR+INIT_HR_FUCR)
	idle::Array{Int64,2} = zeros(STORG_UNITS, 24-INIT_HR_SUCR+INIT_HR_FUCR)
	chrgpwr::Array{Float64,2} = zeros(STORG_UNITS, 24-INIT_HR_SUCR+INIT_HR_FUCR)
	discpwr::Array{Float64,2} = zeros(STORG_UNITS, 24-INIT_HR_SUCR+INIT_HR_FUCR)
	soc::Array{Float64,2} = zeros(STORG_UNITS, 24-INIT_HR_SUCR+INIT_HR_FUCR)
end;

#= STRUCT DemandPreprocGens, Also for SUCR model
wk_ahead::Array{Int64,1} 		<- FUCR_WA_Demand
solar_wa::Array{Int64,1}		<- FUCR_WA_SolarG
wind_wa::Array{Int64,1}      <- FUCR_WA_WindG
hydro_wa::Array{Int64,1}     <- FUCR_WA_HydroG
nuclear_wa::Array{Int64,1}   <- FUCR_WA_NuclearG
cogen_wa::Array{Int64,1}     <- FUCR_WA_CogenG
=#

#=
struct  			<- DF_Generators
unit_id				<-UNIT_ID::Array{String}
status_init 		<- StatusInit
power_init 			<-	PowerInit
uptime_init 		<-	UpTimeInit
downtime_init		<-	DownTimeInit
hr_b1 				<- IHRC_B1_HR
...
hr_b7 				<- IHRC_B8_HR
hr_b1_qty 			<- IHRC_B1_Q
...
hr_b7_qty 			<- IHRC_B8_Q
noload_hr 			<- NoLoadHR::Array{Float64,1}
shutdown_cost 		<- ShutdownCost::Array{Int64,1}
hot_startup_fixcost <- HotStartU_FixedCost::Array{Float64,1}
hot_startup_hr 		<- HotStartU_HeatRate::Array{Float64,1}
down_time_min		<- MinDownTime::Array{Int64,1}
up_time_min			<- MinUpTime::Array{Int64,1}
power_maxout		<- MaxPowerOut::Array{Int64,1}
power_minout		<- MinPowerOut::Array{Int64,1}
spin_res_limit		<- SpinningRes_Limit::Array{Int64,1}
ramp_up_limit		<- RampUpLimit::Array{Int64,1}
ramp_shutdn_limit	<- RampShutDownLimit::Array{Int64,1}
isnuclear 			<- Nuclear::Array{Int64,1}
iscogen				<- Cogen::Array{Int64,1}
isNG_fired			<- NaturalGasFired::Array{Int64,1}

------------
Storage, data_storage.csv
name::Array{String}                  <- DF_Storage.Name::Array{String}
power::Array{Int64,1}                <- DF_Storage.Power::Array{Int64,1}
pwr_ratio_min::Array{Float64,1}      <- DF_Storage.MinPwrRatio::Array{Float64,1}
selfdischarge::Array{Int64,1}        <- DF_Storage.SelfDischarge::Array{Int64,1}
trip_effic_up::Array{Float64,1}      <- DF_Storage.TripEfficUp::Array{Float64,1}
trip_effic_dn::Array{Float64,1}      <- DF_Storage.TripEfficDown::Array{Float64,1}
pwr_toenergy_ratio::Array{Float64,1} <- DF_Storage.PowerToEnergRatio::Array{Float64,1}
soc_init::Array{Float64,1}           <- DF_Storage.SOCInit::Array{Float64,1}
----------
=#

#=
peakers = SysGenerator(DF_Generators.UNIT_ID,
            DF_Peakers.StatusInit,
            DF_Peakers.PowerInit,
            DF_Peakers.UpTimeInit,
            DF_Peakers.DownTimeInit,
            DF_Peakers.IHRC_B1_HR,
            DF_Peakers.IHRC_B2_HR,
            DF_Peakers.IHRC_B3_HR,
            DF_Peakers.IHRC_B4_HR,
            DF_Peakers.IHRC_B5_HR,
            DF_Peakers.IHRC_B6_HR,
            DF_Peakers.IHRC_B7_HR,
            DF_Peakers.IHRC_B1_Q,
            DF_Peakers.IHRC_B2_Q,
            DF_Peakers.IHRC_B3_Q,
            DF_Peakers.IHRC_B4_Q,
            DF_Peakers.IHRC_B5_Q,
            DF_Peakers.IHRC_B6_Q,
            DF_Peakers.IHRC_B7_Q,
            DF_Peakers.NoLoadHR,
            DF_Peakers.ShutdownCost,
            DF_Peakers.HotStartU_FixedCost,
            DF_Peakers.HotStartU_HeatRate ,
            DF_Peakers.MinDownTime,
            DF_Peakers.MinUpTime,
            DF_Peakers.MaxPowerOut,
            DF_Peakers.MinPowerOut,
            DF_Peakers.SpinningRes_Limit,
            DF_Peakers.RampUpLimit,
            DF_Peakers.RampShutDownLimit,
            DF_Peakers.Nuclear,
            DF_Peakers.Cogen,
            DF_Peakers.NaturalGasFired,
)

storage = Storage(
	DF_Storage.Name,
	DF_Storage.Power,
	DF_Storage.MinPwrRatio,
	DF_Storage.SelfDischarge,
	DF_Storage.TripEfficUp,
	DF_Storage.TripEfficDown,
	DF_Storage.PowerToEnergRatio,
	DF_Storage.SOCInit
)

=#

max_stars_repo_name: power-systems-op/grace-arpae

max_issues_repo_name: power-systems-op/grace-arpae

max_forks_repo_name: power-systems-op/grace-arpae


text:
# to be able to run this standalone via include()
using PyCall
import SDtoolbox: zndsolve, cvsolve, CJspeed, PostShock_fr, PostShock_eq, ct

include("correlations.jl")

function cell_size(T₁,P₁,X₁::Union{AbstractString,AbstractDict},mech;kwargs...)

    # Find CJ speed
    try # first try CEA method (fastest)
        global cj_speed = CJspeed(P₁,T₁,X₁,mech)::Float64
    catch PyError # if this fails, use min velocity method
        global cj_speed = CJspeed(P₁,T₁,X₁,mech,method="umin")::Float64
    end
    println("CJ velocity: $(round(cj_speed,digits=3)) m/s")

    # Set up gas object
    gas₁ = ct.Solution(mech)
    gas₁.TPX = T₁,P₁,X₁

    # Find equilibrium post shock state for given speed
    gas = PostShock_eq(cj_speed, P₁, T₁, X₁, mech)
    u_cj = cj_speed*gas₁.density/gas.density
    println("u_CJ: $(round(u_cj,digits=3)) m/s")

    # Find frozen post shock state for given speed
    gas = PostShock_fr(cj_speed, P₁, T₁, X₁, mech)
    Tₛ,Pₛ = gas.TP::Tuple{Float64, Float64}
    println("Tₛ: $(round(Tₛ,digits=3)), Pₛ: $(round(Pₛ,digits=3))")

    # Solve ZND ODEs
    out = zndsolve(gas,gas₁,cj_speed,advanced_output=true;kwargs...)

    # Find CV parameters including effective activation energy
    T₊ = Tₛ*1.02
    gas.TPX = T₊,Pₛ,X₁
    CV_out₊ = cvsolve(gas;kwargs...)
    T₋ = Tₛ*0.98
    gas.TPX = T₋,Pₛ,X₁
    CV_out₋ = cvsolve(gas;kwargs...)
    # Approximate effective activation energy for CV explosion
    τ₊ = CV_out₊["ind_time"]
    τ₋ = CV_out₋["ind_time"]
    if τ₊ == zero(τ₊) && τ₋ == zero(τ₋)
        θₑ_CV = 0.0
    else
        θₑ_CV = 1.0/Tₛ*((log(τ₊)-log(τ₋))/((1.0/T₊)-(1.0/T₋)))
    end

    Δᵣ = u_cj / out["max_thermicity_ZND"]
    Δᵢ = out["ind_len_ZND"]
    χ = θₑ_CV * Δᵢ / Δᵣ

    # obtain cell size through correlation
    λ = ng(Δᵢ,χ)::Float64
end


text:
abstract type ModelType end
abstract type GradientRegression <: ModelType end
abstract type L1Regression <: ModelType end
abstract type QuantileRegression <: ModelType end
abstract type MultiClassRegression <: ModelType end
abstract type GaussianRegression <: ModelType end
struct Linear <: GradientRegression end
struct Poisson <: GradientRegression end
struct Logistic <: GradientRegression end
struct L1 <: L1Regression end
struct Quantile <: QuantileRegression end
struct Softmax <: MultiClassRegression end
struct Gaussian <: GaussianRegression end

mutable struct EvoTreeRegressor{T<:AbstractFloat, U<:ModelType, S<:Int} <: MLJBase.Deterministic
    loss::U
    nrounds::S
    λ::T
    γ::T
    η::T
    max_depth::S
    min_weight::T # real minimum number of observations, different from xgboost (but same for linear)
    rowsample::T # subsample
    colsample::T
    nbins::S
    α::T
    metric::Symbol
    seed::S
end

function EvoTreeRegressor(;
    loss=:linear,
    nrounds=10,
    λ=0.0, #
    γ=0.0, # gamma: min gain to split
    η=0.1, # eta: learning rate
    max_depth=5,
    min_weight=1.0, # minimal weight, different from xgboost (but same for linear)
    rowsample=1.0,
    colsample=1.0,
    nbins=64,
    α=0.5,
    metric=:mse,
    seed=444)

    if loss == :linear model_type = Linear()
    elseif loss == :logistic model_type = Logistic()
    elseif loss == :L1 model_type = L1()
    elseif loss == :quantile model_type = Quantile()
    end

    model = EvoTreeRegressor(model_type, nrounds, λ, γ, η, max_depth, min_weight, rowsample, colsample, nbins, α, metric, seed)

    return model
end

mutable struct EvoTreeCount{T<:AbstractFloat, U<:ModelType, S<:Int} <: MLJBase.Probabilistic
    loss::U
    nrounds::S
    λ::T
    γ::T
    η::T
    max_depth::S
    min_weight::T # real minimum number of observations, different from xgboost (but same for linear)
    rowsample::T # subsample
    colsample::T
    nbins::S
    α::T
    metric::Symbol
    seed::S
end

function EvoTreeCount(;
    loss=:poisson,
    nrounds=10,
    λ=0.0, #
    γ=0.0, # gamma: min gain to split
    η=0.1, # eta: learning rate
    max_depth=5,
    min_weight=1.0, # minimal weight, different from xgboost (but same for linear)
    rowsample=1.0,
    colsample=1.0,
    nbins=64,
    α=0.5,
    metric=:poisson,
    seed=444)

    if loss == :poisson model_type = Poisson() end
    model = EvoTreeCount(model_type, nrounds, λ, γ, η, max_depth, min_weight, rowsample, colsample, nbins, α, metric, seed)

    return model
end

mutable struct EvoTreeClassifier{T<:AbstractFloat, U<:ModelType, S<:Int} <: MLJBase.Probabilistic
    loss::U
    nrounds::S
    λ::T
    γ::T
    η::T
    max_depth::S
    min_weight::T # real minimum number of observations, different from xgboost (but same for linear)
    rowsample::T # subsample
    colsample::T
    nbins::S
    α::T
    metric::Symbol
    seed::S
end

function EvoTreeClassifier(;
    loss=:softmax,
    nrounds=10,
    λ=0.0, #
    γ=0.0, # gamma: min gain to split
    η=0.1, # eta: learning rate
    max_depth=5,
    min_weight=1.0, # minimal weight, different from xgboost (but same for linear)
    rowsample=1.0,
    colsample=1.0,
    nbins=64,
    α=0.5,
    metric=:mlogloss,
    seed=444)

    if loss == :softmax model_type = Softmax() end
    model = EvoTreeClassifier(model_type, nrounds, λ, γ, η, max_depth, min_weight, rowsample, colsample, nbins, α, metric, seed)

    return model
end

mutable struct EvoTreeGaussian{T<:AbstractFloat, U<:ModelType, S<:Int} <: MLJBase.Probabilistic
    loss::U
    nrounds::S
    λ::T
    γ::T
    η::T
    max_depth::S
    min_weight::T # real minimum number of observations, different from xgboost (but same for linear)
    rowsample::T # subsample
    colsample::T
    nbins::S
    α::T
    metric::Symbol
    seed::S
end

function EvoTreeGaussian(;
    loss=:gaussian,
    nrounds=10,
    λ=0.0, #
    γ=0.0, # gamma: min gain to split
    η=0.1, # eta: learning rate
    max_depth=5,
    min_weight=1.0, # minimal weight, different from xgboost (but same for linear)
    rowsample=1.0,
    colsample=1.0,
    nbins=64,
    α=0.5,
    metric=:gaussian,
    seed=444)

    if loss == :gaussian model_type = Gaussian() end
    model = EvoTreeGaussian(model_type, nrounds, λ, γ, η, max_depth, min_weight, rowsample, colsample, nbins, α, metric, seed)

    return model
end

# For R-package
function EvoTreeRModels(
    loss,
    nrounds,
    λ,
    γ,
    η,
    max_depth,
    min_weight,
    rowsample,
    colsample,
    nbins,
    α,
    metric,
    seed)

    if loss ∈ [:linear, :L1, :logistic, :quantile]
        if loss == :linear model_type = Linear()
        elseif loss == :logistic model_type = Logistic()
        elseif loss == :quantile model_type = Quantile()
        elseif loss == :L1 model_type = L1()
        end
        model = EvoTreeRegressor(model_type, nrounds, λ, γ, η, max_depth, min_weight, rowsample, colsample, nbins, α, metric, seed)
    elseif loss == :poisson
        model = EvoTreeCount(Poisson(), nrounds, λ, γ, η, max_depth, min_weight, rowsample, colsample, nbins, α, metric, seed)
    elseif loss == :softmax
        model = EvoTreeClassifier(Softmax(), nrounds, λ, γ, η, max_depth, min_weight, rowsample, colsample, nbins, α, metric, seed)
    elseif loss == :gaussian
        model = EvoTreeGaussian(Gaussian(), nrounds, λ, γ, η, max_depth, min_weight, rowsample, colsample, nbins, α, metric, seed)
    else
        throw("invalid loss")
    end

    return model
end

const EvoTypes = Union{EvoTreeRegressor,EvoTreeCount,EvoTreeClassifier,EvoTreeGaussian}