rotmat Function

    function rotmat(axis, angle)
    !! Rotate around around an axis by a given angle taken from [here](http://www.neilmendoza.com/glsl-rotation-about-an-arbitrary-axis/)
        use utils, only : deg2rad
        !> Axis to rotate around
        type(vector),  intent(in) :: axis
        !> Angle to rotate by in degrees
        real(kind=wp), intent(in) :: angle

        type(vector) :: axist

        real(kind=wp) :: rotmat(4, 4), s, c, oc, a

        axist = axis%magnitude()
        a = deg2rad(angle)

        s = sin(a)
        c = cos(a)
        oc = 1._wp - c

    rotmat(:, 1) = [oc * axist%x * axist%x + c, oc * axist%x * axist%y - axist%z * s,&
                    oc * axist%z * axist%x + axist%y * s, 0.0_wp]
    rotmat(:, 2) = [oc * axist%x * axist%y + axist%z * s, oc * axist%y * axist%y + c,&
                    oc * axist%y * axist%z - axist%x * s, 0.0_wp]
    rotmat(:, 3) = [oc * axist%z * axist%x - axist%y * s, oc * axist%y * axist%z + axist%x * s,&
                    oc * axist%z * axist%z + c, 0.0_wp]
    rotmat(:, 4) = [0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp]

    end function rotmat