generating two orthogonal vectors that are orthogonal to a particular direction

Sorry, I can't put it as a comment because of a lack of reputation.

Regarding @behzad.nouri's answer, note that if k is not a unit vector the code will not give an orthogonal vector anymore!

The correct and general way to do so is to subtract the longitudinal part of the random vector. The general formula for this is here

So you simply have to replace this in the original code:

>>> x -= x.dot(k) * k / np.linalg.norm(k)**2

Assume the vector that supports the orthogonal basis is u.

b1 = np.cross(u, [1, 0, 0])   # [1, 0, 0] can be replaced by other vectors, just get a vector orthogonal to u
b2 = np.cross(u, b1)
b1, b2 = b1 / np.linalg.norm(b1), b2 / np.linalg.norm(b2)

A shorter answer if you like.

Get a transformation matrix

B = np.array([b1, b2])
TransB = np.dot(B.T, B)
u2b = TransB.dot(u) # should be like [0, 0, 0]

This will do:

>>> k  # an arbitrary unit vector k is not array. k is must be numpy class. np.array
np.array([ 0.59500984,  0.09655469, -0.79789754])

To obtain the 1st one:

>>> x = np.random.randn(3)  # take a random vector
>>> x -= x.dot(k) * k       # make it orthogonal to k
>>> x /= np.linalg.norm(x)  # normalize it

To obtain the 2nd one:

>>> y = np.cross(k, x)      # cross product with k

and to verify:

>>> np.linalg.norm(x), np.linalg.norm(y)
(1.0, 1.0)
>>> np.cross(x, y)          # same as k
array([ 0.59500984,  0.09655469, -0.79789754])
>>> np.dot(x, y)            # and they are orthogonal
-1.3877787807814457e-17
>>> np.dot(x, k)
-1.1102230246251565e-16
>>> np.dot(y, k)
0.0