I remembered that the isometries for the upper half space model of $\mathbf H^3$ are given by the action of $PSL_2\mathbf C$ on the points at infinity (the $xy$-plane as the complex numbers). So, given four geodesic endpoints, move one to $0$ with a pure translation. We now have three complex numbers $A,B,C,$ where $A$ is the other endpoint of the geodesic with an endpoint at $0.$ Define a complex number $\gamma$ that solves
$$ ABC \gamma^2 + 2 BC \gamma + (B+C-A) =0, $$
where you might as well pick $\gamma = 0$ if the constant term $B+C-A=0.$
Next, apply the linear fractional transformation
$$ h(z) = \frac{z}{\gamma z + 1} $$
to the plane.
The result is $$ h(A) = h(B) + h(C) $$
That is, the midpoint of $0$ and $h(A)$ is the same as the midpoint of $h(B)$ and $h(C).$ So the common orthogonal geodesic is the vertical ray through $h(A)/2$ and allowing the third coordinate to vary. Then map everything back to your originals.
Note that there is no answer if two of your original geodesic endpoints coincide. In that case, they lie in a common flat, and lines asymptotic at infinity do not share a perpendicular.
See if I can do the link the right way this time, Milnor's survey on the first 150 years of hyperbolic geometry can be downloaded.
