Rotate a unit sphere such as to align it two orthogonal unit vectors

Let $q = a \times b$, let $w = c \times d$.

Then the matrix $M$ whose columns are $a, b, q$, sends the standard basis vectors $e_1 = \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}$, $e_2 = \begin{bmatrix} 0 \\ 1 \\ 0 \end{bmatrix}$, and $e_3 = \begin{bmatrix} 0 \\ 0 \\ 1 \end{bmatrix}$ to $a, b,$ and $q$ respectively. Since each of these triples is a positively oriented orthonormal frame, the map "multiplication by $M$" must be a rotation.

Similarly, the matrix $K$, whose columns are $c, d,$ and $w$, sends the standard basis vectors to those three vectors.

Therefore

$$ v \mapsto KM^t v $$ sends $a$ to $c$ (because $M^t$ sends $a$ to $e_1$, and $K$ sends $e_1$ to $c$), $b$ to $d$ and $q$ to $w$ as needed.

And because both $K$ and $M$ are rotations, the product matrix $KM^t$ is also a rotation.

If you actually need to know the axis and angle of the rotation rather than just a matrix for the rotation, let me know and I can expand on this a little bit.

Addition after comments OP Asks for a way to find the axis and angle from the rotation matrix, so here we go. If $M$ is a rotation matrix that represents rotation by angle $\theta$ about axis $v$, then $$ \theta = \cos^{-1} \left(\frac{tr(M) -1}{2} \right). $$ where $tr(M)$ denotes the sum of the diagonal entries. You'll notice that if $tr(M) = 2$ (i.e., $M$ is the identity), then you get $\theta = 0$.

At this point, there are three cases:

1. If $\theta = 0$, then any unit vector serves as an axis for the rotation.

2. If $\theta = \pi$, then look at $M + I$, for which at least one column is nonzero. Take such a nonzero column and normalize it (i.e., divide by its length) to get a unit vector $v$ which is the axis of the rotation. (You might want to try this for the diagonal matrix with diagonal entries $-1, -1, 1$ just to see what happens).

3. For all other $\theta$, compute $Q = \frac{M - M^t}{2 \sin \theta}$; this will be skew-symmetric matrix (i.e., its transpose is its negative). Letting $x = -q_{23}, y = q_{13}, z = -q_{12}$, the vector $$ v = \begin{bmatrix} x\\y\\z\end{bmatrix} $$ is the axis of rotation.

Code for this, in pidgin C#, is below. (Note that C# uses 0-based indexing, while my math above uses 1-based indexing). The output vector is called "omega" instead of "v".

void RotationToAxisAngle(
    Mat33 m,
    out Vector3D omega,
    out double theta)
{
    // convert a 3x3 rotation matrix m to an axis-angle representation
    theta = Math.acos( (m.trace()-1)/2);
    if (theta is near zero)
    {
        omega = Vector3D(1,0,0); // any vector works
        return;
    }
    if (theta is near pi)
    {
        int col = column with largest entry of m in absolute value;
        omega = Vector3D(m[0, col], m[1, col], m[2, col]);
        return;
    }
    else
    {
        mat 33 s = m - m.transpose();
        double x = -s[1,2], y = s[0,2]; z = -s[0,1];
        double t = 2 * Math.Sin(theta);
        omega = Vector3D(x/t, y/t, z/t);
        return;
    }
}

Of course, if you can do it with two or more rotations, you can do it with one: the composition of rotations is a rotation.

Conceptually the strategy should be clear: $a$ and $c$ lie in a plane, and you can find a rotation $R_1$ that turns that plane carrying $a$ onto $c$. After that, $R_1(b)$ and $d$ lie in the plane perpendicular to $R_1(a)=c$, and you just need to rotate around $c$ through some angle with $R_2$ to get $b$ to lie on $d$. Then the rotation $R_2R_1$ carries the first orthonormal pair onto the second pair.

It would be great exercise for you to work the details of this strategy out with matrices. Just start with a rotation around $a\times c$ and follow up with a rotation around $c$.

I'm not sure if you're familiar with doing such things in terms of quaternions, but that would also be a great exercise.