rotationAlign Function

    function rotationAlign(a, b) result(res)
        !! Calculate the rotation matrix to rotate vector a onto b
        !! [ref1](https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula)
        !! [ref2](https://math.stackexchange.com/questions/180418/calculate-rotation-matrix-to-align-vector-a-to-vector-b-in-3d)
        
        !> Vector to rotate. Unit vector
        type(vector), intent(in) :: a
        !> Vector to be rotated onto. Unit vector
        type(vector), intent(in) :: b
        
        type(vector)  :: v
        real(kind=wp) :: c, k, res(4, 4), v_x(4, 4), v_x2(4, 4)

        v = a .cross. b
        c = a .dot. b
        k = 1._wp / (1._wp + c)
        
        !skew-symmetric matrix
        v_x(:, 1) = [0._wp     , -1._wp*v%z, v%y       , 0._wp]
        v_x(:, 2) = [v%z       , 0._wp     , -1._wp*v%x, 0._wp]
        v_x(:, 3) = [-1._wp*v%y, v%x       , 0._wp     , 0._wp]
        v_x(:, 4) = [0._wp, 0._wp, 0._wp, 0._wp]
        
        v_x2 = matmul(v_x, v_x)
        res = identity() + v_x + v_x2*k

    end function rotationAlign