rotation_matrix_dot Function

    pure function rotation_matrix_dot(axis,angle,angledot) result(rotmatdot)

    implicit none

    real(wp),dimension(3,3) :: rotmatdot   !! the rotation matrix derivative \( d \mathbf{C} / d t \)
    integer,intent(in)      :: axis        !! x_axis, y_axis, or z_axis
    real(wp),intent(in)     :: angle       !! angle in radians
    real(wp),intent(in)     :: angledot    !! time derivative of angle in radians/sec

    real(wp) :: c,s

    !precompute these:
    c = cos(angle)
    s = sin(angle)

    !first compute d[C]/da (time derivate w.r.t. the angle):
    select case (axis)
    case(x_axis); rotmatdot = reshape([zero, zero, zero, zero, -s, -c, zero, c, -s],[3,3])
    case(y_axis); rotmatdot = reshape([-s, zero, c, zero, zero, zero, -c, zero, -s],[3,3])
    case(z_axis); rotmatdot = reshape([-s, -c, zero, c, -s, zero, zero, zero, zero],[3,3])
    case default
        rotmatdot = zero
        return
    end select

    rotmatdot = rotmatdot * angledot    ! d[C]/dt = d[C]/da * da/dt

    end function rotation_matrix_dot