The two given cones share the same vertex which is the origin. The axis of the first is the unit vector $u_1 = (a,b,c)$ while the second's axis is the unit vector $u_2 = (a',b',c')$. The semi-apical angle of the first is $\theta_1$ and the second's is $\theta_2$.
The two lines of the intersection, are those vectors $r = (x,y,z)$ that make an angle of $\theta_1$ with $u_1$ and $ \theta_2$ with $u_2$.
Imposing $r$ to have a unit length results in the following three equations
$ a x + b y + c z = \cos \theta_1 $
$ a' x + b' y + c' z = \cos \theta_2 $
$ x^2 + y^2 + z^2 = 1 $
The solution of this system is relatively easy because the first two equations are linear in $x,y,z$. Their solution is the line
$ L(\lambda) = V_0 + \lambda V_1 $
where $ V_1 = (a,b,c) \times (a',b',c') $ is the direction vector, and
$ V_0 $ is any point that solves both equations. Such a point can be found by setting $z = 0 $ and solving the remaing $2 \times 2$ system. Doing so, gives
$ x_1 = \dfrac{ b' \cos \theta_1 - b \cos \theta_2 }{ a b' - a' b } $
$ y_1 = \dfrac{ a \cos \theta_2 - a' \cos \theta_1 }{a b' - a' b } $
So that
$ V_0 = (x_1, y_1, 0) $
Now plugging $L(\lambda)$ into the third equation, gives
$ (V_0 + \lambda V_1) \cdot (V_0 + \lambda V_1) = 1 $
which is
$ \lambda^2 \ (V_1 \cdot V_1) + 2 \lambda (V_0 \cdot V_1) + (V_0 \cdot V_0 - 1) = 0 $
This will have two possible solutions, which will give the two unit vectors $r$.
As an explicit example, suppose the unit vectors $(a,b,c)$ and $(a',b',c') $ are given by
$ (a,b,c) = \dfrac{1}{3}(1,2,2) $
$ (a',b',c') = \dfrac{1}{9}(1, -4, 8) $
And let
$ \theta_1 = \dfrac{\pi}{3} $
$\theta_2 = \dfrac{\pi}{4}$
Applying the above formulas,
$V_1 = (a,b,c) \times (a',b',c') = \dfrac{1}{9} (8, -2, -2) $
And since this a direction vector, we can scale it, so we'll take
$ V_1 = (4, -1, -1) $
As $V_0$, we have
$x_1 = 1 +\dfrac{3}{2} \sqrt{2} $
$ y_1 = -\dfrac{3}{ 4} \sqrt{2} + \dfrac{1}{4} $
so that $V_0 = (x_1, y_1, 0) $
Using these two vector into the third (quadratic) equation, gives two roots.
$ \lambda_{1,2} = -0.89167 , -0.58565 $
And these two values of $\lambda$ generate the following unit vectors solutions to the problem.
$ u_1 = \begin{pmatrix} -0.44537 \\ 0.081013 \\ 0.891673 \end{pmatrix} \hspace{25pt} u_2 = \begin{pmatrix} 0.778705 \\ -0.22501 \\ 0.585654 \end{pmatrix} $