Determine the rotation necessary to transform one point on a sphere to another. Given two points on the surface of a sphere in 3D space, is it possible to/how does one systematically determine the set of rotations around the x, y, and z-axes in order to transform one point to the other?
For example, on the unit sphere centred at the origin with a radius of 1, the points $(1, 0, 0)$ and $(0, \frac{1}{\sqrt{2}}, \frac{1}{\sqrt{2}})$ are both (I believe) on the sphere's surface. One can translate from $(1, 0, 0)$ to $(0, \frac{1}{\sqrt{2}}, \frac{1}{\sqrt{2}})$ by rotating $90^o$ degrees (counter-clockwise?) around the z-axis, then rotating $45^0$ degrees (clockwise?) around the x-axis.
I feel like the answer to this should be simple, but I have been struggling to answer it for a few days now. 
EDIT: Just to add my attempt to follow Robert Israel's approach below with my real world example; perhaps someone can spot my mistake.
My input point is $(0.248729, -0.617981, 0.07822)$, the target is $(0.31209, -0.524457, 0.07822)$ - note that these both have a length (distance from $(0, 0, 0)$) of $0.670735...$, so are on the same sphere (give or take floating point errors.)
$a = 0.248739,\ b = -0.617981,\ c =0.07822,\ d = 0.31209,\ e = -0.524457,\ f =0.07822$
I compute $r = \sqrt{a^2 + b^2} = 0.6615811...$.
Therefore, $\theta = \operatorname{sgn}(b)\arccos(a/r) = -1 * 1.1881482... = -1.1881482...$
$|d|$ is $0.31209$, which is less than $r$, so $\alpha$ can just be computed as $\alpha = \arccos(d/r) - \theta = 2.271361...$
Rotating the original point by this, to $(d, g, h)$, gives $(0.31209, 0.588529, 0.078229)$ - this is good, since $d$ is correct and the length of this is also $0.670735...$, so it's on the sphere. Additionally, $g^2 + h^2 = 0.3524... = e^2 + f^2$.
However, this is where I flounder. I attempt to calculate $\phi$:
$g = \sqrt{g^2 + h^2} \,\cos(\phi)$
$\phi = \arccos\left(\sqrt{g^2 + h^2}/g\right)$
However, $\sqrt{g^2 + h^2}/g = 1.0087...$, so it doesn't work to take the $\arccos$ of it... and now I don't know where to go.
 A: In general, suppose you want to start with $(a,b,c)$ and transform it to $(d,e,f)$ (where $a^2 + b^2 + c^2 = d^2 + e^2 + f^2 = 1$) using a rotation around the $z$ axis followed by a rotation around the $x$ axis.
Now rotation by angle $\alpha$ around the $z$ axis takes $(a,b,c)$ to $(a \cos(\alpha) - b \sin(\alpha), a \sin(\alpha) + b \cos(\alpha), c)$.  Since a rotation around the $x$ axis preserves the $x$ coordinate, you need $a \cos(\alpha) - b \sin(\alpha) = d$.  Write $a = r \cos(\theta)$, $b = r \sin(\theta)$ where $r = \sqrt{a^2 + b^2}$ and $\theta = \text{signum}(b)  \arccos(a/r)$, and this equation becomes $r \cos(\alpha + \theta) = d$.  For this to be possible you need $|d|\le r$, and then you can take $\alpha = \arccos(d/r) - \theta$.  If the result of this rotation is $(d, g, h)$, we will have $g^2 + h^2 = 1 - d^2 = e^2 + f^2$, and rotation by angle $\beta$ around the $x$ axis will take this to
$(d, g \cos(\beta) - h \sin(\beta), g \sin(\beta) + h \cos(\beta))$.  Write $g = \sqrt{g^2+h^2} \cos(\phi)$, $h = \sqrt{g^2 + h^2} \sin(\phi)$, $e = \sqrt{g^2+h^2} \cos(\tau)$ and $f = \sqrt{g^2+h^2} \sin(\tau)$, and you'll see that $\beta = \tau - \phi$ will work.     
Of course you can always throw in a rotation around the $y$ axis too (and you'll need one if $|d| > r$).
A: As is mentioned here, another approach to the problem is to use quaternions. The basic wiki page and the rotation wiki page are very good resources for understanding how these generalizations of complex numbers work. The procedure is as follows:


*

*Compute the cross product of the vectors going to the two points: this will be the axis of rotation, and gets inserted into the $\mathbf{i}, \mathbf{j},$ and $\mathbf{k}$ components of the rotation quaternion. 

*Compute the dot product of the vectors going to the two points: this allows you to find the angle of rotation needed, which goes into the real part of the quaternion, as well as multiplying the complex parts. The angle is the $\arccos$ of the dot product (recall that $\mathbf{x}\cdot\mathbf{y}=\|x\|\,\|y\|\,\cos(\theta)$, so if $\|x\|=\|y\|=1$, we have $\theta=\arccos(\mathbf{x}\cdot\mathbf{y})$).


So, for your example: you want to find the rotation needed to get from point $\mathbf{v}_1=(1,0,0)$ to point $\mathbf{v}_2=\left(0,\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}}\right)$. We compute: 
\begin{align*}
\mathbf{v}_1\times\mathbf{v}_2&=\left(0,-\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}}\right)\quad\text{and}\\
\mathbf{v}_1\cdot\mathbf{v}_2&=0.
\end{align*}
Hence, the angle of rotation is $\theta=\arccos(0)=\dfrac{\pi}{2}.$ From the rotation quaternion wiki page, we see that the necessary rotation quaternion is
\begin{align*}q&=\cos\left(\frac{\theta}{2}\right)+(u_x\mathbf{i}+u_y\mathbf{j}+u_z\mathbf{k})\,\sin\left(\frac{\theta}{2}\right)\\
&=\frac{1}{\sqrt{2}}+\left(0\mathbf{i}-\frac{\mathbf{j}}{\sqrt{2}}+\frac{\mathbf{k}}{\sqrt{2}}\right)\frac{1}{\sqrt{2}}\\
&=\frac{1}{\sqrt{2}}+\left(0\mathbf{i}-\frac{\mathbf{j}}{2}+\frac{\mathbf{k}}{2}\right).\end{align*}
Now we double-check that the rotation represented by this quaternion actually works, and here's where the fun really begins: quaternion multiplication! The rotation from $\mathbf{v}_1$ to $\mathbf{v}_2$ ought to work as follows:
$$\mathbf{v}_2=q\mathbf{v}_1q^{-1}.$$
Here I'm going to refer to the Wolfram Language (Mathematica), as implemented on the Wolfram Dev Platform. We use the following code:
<<Quaternions`
q=Quaternion[1/Sqrt[2],0,-1/2,1/2]
v1=Quaternion[0,1,0,0]
NonCommutativeMultiply[q,NonCommutativeMultiply[v1,Conjugate[q]]]

Mathematica produces the answer
$$\operatorname{Quaternion}\left[0,0,\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}}\right]$$
Notice here we must use the NonCommutativeMultiply function to get the correct quaternion multiplication. Also notice that for a unit quaternion (which rotation quaternions are), $q^{-1}=\overline{q}$. Finally, notice that we must "quaternion-ify" the vector $\mathbf{v}_1$ as the quaternion $(0,v_{1x},v_{1y},v_{1z})$, and the result is a quaternion which we can "vector-ify" by extracting the last three components. The result, as you can see, is correct.
If you prefer to have the rotation matrix, you can get such an $\mathbf{R}$ as explained on the rotation quaternion wiki page under the section "Quaternion-derived rotation matrix". 
A: I interpret the problem statement to mean that we are allowed to specify rotations only around the three orthogonal coordinate axes, and may not rotate around any other axis.
I will also assume we always want to be able to move from one point to the other in no more than two rotations. The problem is almost trivially easy
if we are allowed three rotations.
Let the starting point be $(x_1,y_1,z_1)$ and the desired ending point be 
$(x_2,y_2,z_2).$
We consider the six planes through those two points parallel to the coordinate planes,
that is, the six planes whose equations are respectively 
$x = x_1,$ $y = y_1,$ $z = z_1,$ $x = x_2,$ $y = y_2,$ and $z = z_2.$
These planes partition three-dimensional space into $27$ compartments,
including a parallelepiped (rectangular box) of which
$(x_1,y_1,z_1)$ and $(x_2,y_2,z_2)$ are diagonally opposite vertices.
A suitable rotation around the $x$ axis can put the image of $(x_1,y_1,z_1)$ anywhere on the circle in the plane $x = x_1$ that is centered on the $x$-axis and passes through $(x_1,y_1,z_1).$
Suppose that circle intersects the circle in the plane $y = y_2$ that is centered on the $y$-axis and passes through $(x_2,y_2,z_2),$
and suppose $(x_3,y_3,z_3)$ is a point of intersection of those two circles.
Observe that since one circle lies entirely in the plane $x = x_1$ and the other lies entirely in the plane $y = y_2,$
it must be that $x_3 = x_1$ and $y_3 = y_2.$
Then we can rotate around the $x$-axis to put the image of $(x_1,y_1,z_1)$
at $(x_1,y_2,z_3),$ after which we can rotate around the $y$-axis so that the image of $(x_1,y_2,z_3)$ is $(x_2,y_2,z_2).$
If there is a second distinct point of intersection of the two circles,
$(x_1,y_2,z_4),$ then there is another sequence of rotations about the $x$- and $y$-axes that take $(x_1,y_1,z_1)$ to $(x_1,y_2,z_4),$ and then to $(x_2,y_2,z_2).$
In order for us to be able to do the two particular rotations described in the previous paragraph, it is sufficient that the sphere intersects the line given by the pair of equations $x = x_1, y = y_2.$
The number of intersections with that line gives the number of intersections of the two circles.
There are a total of six lines (call them $L_1,$ $L_2,$ $L_3,$ $L_4,$ $L_5,$ and $L_6$) given by pairs of equations where the first equation is one of the equations $x = x_1,$ $y = y_1,$ or $z = z_1$
and the second equation is one of the equations 
$x = x_2,$ $y = y_2,$ or $z = z_2$; provided that the sphere intersects at least one of those lines, some variation of the procedure given in the previous paragraph will take $(x_1,y_1,z_1)$ to $(x_2,y_2,z_2)$
in just two rotations.
To show that the sphere must intersect at least one of the six lines
$L_i, 1\leq i\leq 6,$ consider the plane through the four points $(x_1,y_1,z_1),$ $(x_2,y_2,z_2),$ $(x_1,y_1,z_2),$ and $(x_2,y_2,z_1).$
The sphere intersects that plane in a circle that passes through the points
$(x_1,y_1,z_1)$ and $(x_2,y_2,z_2).$ 
At least one of the points $(x_1,y_1,z_2)$ and $(x_2,y_2,z_1)$ lies inside or on that circle, and therefore is inside or on the sphere.
The sphere therefore intersects 
the lines $x = x_1, z = z_2$ and $y = y_1, z = z_2,$
or the lines $x = x_2, z = z_1$ and $y = y_2, z = z_1,$ or all four lines.
The number of distinct ways that we can transform
$(x_1,y_1,z_1)$ to $(x_2,y_2,z_2)$ by two rotations around the coordinate axes is equal to the number of distinct intersections of the sphere with the six lines $L_i, 1\leq i\leq 6.$
(A point at the intersection of two of these lines and the sphere counts as one intersection.)
The solution is never unique, because (as already seen) at least two of the lines $L_i, 1\leq i\leq 6$ intersect the sphere, and it is not possible for the sphere to be tangent to both lines at one point;
hence the two lines give at least two points of intersection, each of which gives a different pair of rotations about the coordinate axes.
A: I came up with this as an answer for a different question, but I think it's more appropriate to post it here. What follows is a comprehensive explanation of Robert's approach, as I understand it, using notation that I find more straightforward.
The plan is this: we will move point 1 to point 2 using a rotation about the z-axis, followed by a rotation about the $y$-axis if needed, followed by a rotation about the x-axis. This is possible as long as both points have the same "length" (distance from the origin), so for simplicity we'll scale things such that both points have length $1$. Write the coordinates as follows:
$$
P_1: (x_1,y_1,z_1), \quad P_2:(x_2,y_2,z_2)
$$
First, we'll try to do without the $y$-axis rotation. The idea is as follows: ultimately, we're trying to move point 1 to point 2, following some path along the sphere. There is an "orbit" of points that can be reached from or brought to point 1 (shown in blue in the figure below) by rotating about the $z$-axis (the vertical axis) and an orbit of points that can be reached from or brought to point 2 (show in orange in the figure below) by rotating about the $x$-axis (the axis pointing to the "right"). All points on the blue circle have the same $z$-coordinate as point 1, and all points on the orange circle have the same $x$-coordinate as point 2.

Our first goal is to rotate point 1 so that lies in the orbit of point 2. That is, we want to rotate point 1 so that its $x$-coordinate matches that of point 2. Let $\theta$ denote the azimuthal angle of $P_1$ (that is, $\theta = \operatorname{arctan2}(y_1,x_1)$). Let $\alpha_z$ denote the angle of the rotation (about the $z$-axis) that we're looking for. That is, the angle required to bring $P_1$ into the orbit of $P_2$. The equation to be solved is
$$
\sqrt{x_1^2 + y_1^2} \cos(\alpha_z + \theta) = x_2 \implies \alpha_z = \arccos \left( \frac{x_2}{\sqrt{x_1^2 + y_1^2}}\right) - \theta.
$$
Note that $\arccos$ is only defined for arguments between $1$ and $-1$. Consequently, we can only reach the orbit of $P_2$ if $\left|x_2\big/\sqrt{x_1^2 + y_1^2} \right| < 1$, which is to say that
$$
x_2^2 \leq x_1^2 + y_1^2. 
$$
We'll assume for now that this holds and that we get a solution $\alpha_z$. Let $(x_3,y_3,z_3)$ denote the coordinates of the points that we reach with this first rotation (note that $z_3 = z_1$ and $x_3 = x_2$). The angle of the second rotation (about the x-axis) that we'll need is
$$
\alpha_x = \operatorname{arctan2}(z_2,y_2) - \operatorname{arctan2}(z_3,y_3).
$$
In the case that we can't reach the orbit of $P_2$ about the $x$-axis, we can "bridge the gap" using a rotation about the $y$-axis; see the image below.

First, rotate $P_1$ so that it lies within the $xz$-plane. The required rotation here is by the angle $\alpha_z = -\theta$, with $\theta$ as above. That is,
$$
\alpha_x = -\operatorname{arctan2}(y_1,x_1).
$$
Next, we rotate the resulting point onto the orbit of $P_2$. The required angle can simply be computed as
$$
\alpha_y = \arcsin(x_2) - \arccos(z_1).
$$
Finally, we rotate from this point to reach $P_2$. The required angle is
$$
\alpha_x = \operatorname{arctan2}(z_2,y_2).
$$
In fact, we could also apply this technique in the case that the orbits intersect, as in the figure below. The formulas for the angles are the same; notably, $\alpha_y$ is now negative. The only "caveat" here, I suppose, is that our rotation around the $y$-axis was not strictly necessary.

A: While considering another problem, I ended up with another approach that leads to the same rotation as that obtained in G Cab's and Adrian's answers.
Let $v_1$ denote the first point and $v_2$ the second, and suppose that $\|v_1\| = \|v_2\| = 1$. Let $u = (x_u,y_u,z_u)$ denote the unit vector parallel to $v_1 \times v_2$, $u = \frac{v_1\times v_2}{\|v_1 \times v_2\|}$. With Rodrigues' rotation formula, we find that this rotation has the matrix $R = I + \sin\theta K + (1-\cos \theta) K^2$, where $\theta = \arccos^{-1}v_1 \cdot v_2$ $K$ denotes the "cross-product matrix" for the vector $u$,
$$
K = \pmatrix{
0 & -z_u & y_u\\
z_u&0&x_u\\
-y_u&x_u&0}.
$$
From there, the $z$-$x$-$Z$ extrinsic Euler angles for the rotation can be obtained using the formulas given here:
\begin{align} 
 \phi   &=  \operatorname{arctan2}\left(R_{31}, R_{32}\right),\\
 \theta &=  \arccos\left(R_{33}\right),\\
 \psi   &= -\operatorname{arctan2}\left(R_{13}, R_{23}\right).
\end{align}
As it turns out, we can neatly (without any trigonometric functions!) rewrite the formula for $R$ using the fact that $u = \frac 1{\sin\theta}v_1 \times v_2$ and $\cos \theta = v_1^\top v_2$. I claim (without proof for now) that
$$
R = (v_1^\top v_2) I + (v_2v_1^\top - v_1v_2^\top) + (1 + v_1^\top v_2)(v_1 \times v_2)(v_1 \times v_2)^\top.
$$
Note that $v_1,v_2,v_1 \times v_2$ are used as column-vectors in the above. Another nice reformulation:
$$
M = \pmatrix{v_1 & v_2 & v_1\times v_2}, \\
R = (v_1^\top v_2) I + M\pmatrix{0&-1&0\\1&0&0\\0&0&1 + v_1^\top v_2}M^\top.
$$
