“aim at” in 3D space algorithm

I'd to calculate the cartesian rotate values of "aiming" one point at another in 3D space.

For example, assuming rotation order of XYZ:
Point A is at [0, 0, 0] Point B is at [2, 2, 0]

If point A was "aiming" at point B with primary (aim) axis X, secondary (up) Y; point A's rotate values would be [0, 0, 45]. These values are what I'd like to calculate.

• Using a rotation order is often asking for trouble. If your system uses transformation matrices for each object, then to have object A aim at object B along A's X-axis, you can just normalize (POSB-POSA) and create a orthonormal bases using that as X, and by taking cross product of (0,0,1) with X to create Y. Normalize that Y, and do X cross Y to get your Z. Now you have 3 orthonormal axes that form your transformation matrix for object A. Just plug in the translation of A, and you are done. No angles or rotations required. No chance of messing up rotation orders. – Bram Jul 2 at 22:57

It's assumed that point A is at the origin. Let d be the unit vector in the direction point A is currently looking at, and p be the unit vector in the direction of B - A.

Notice that d and p necessarily define a plane in 3D-space. If you can imagine this, then it becomes clear that a single rotation of d by an amount equal to the angle between these two vectors about the normal to this plane will orient point A towards point B. The goal is to find out what the angle between these two vectors and the normal to the plane they define are. (Don't worry, I will get to Euler Angles soon!)

With basic linear algebra we have (I am assuming a familiarity with these operations):

acos(dot(d,p)) = the angle between d and p, I'll call it theta.

cross(d, p)/length(cross(d,p)) = unit normal to the plane defined by d and p, I'll call it n.

With an angle to rotate by and an axis of rotation, we can define a quaternion to represent this rotation. (If you aren't familiar with quaternions at all, it's fine; they're only used as an intermediate step here to get to the Euler angles you want.)

This quaternion will be equal to:

q = cos(theta/2) + (nxi + nyj + nzk) * sin(theta/2), where nx is the x-component of n, and so on.

Finally, there is the formula for converting from quaternion to Euler angle representation of a rotation.

$$\begin{bmatrix} \phi \\ \theta \\ \psi \end{bmatrix} = \begin{bmatrix} atan2(2(q_0q_1+q_2q_3),1-2(q^2_1+q^2_2)) \\ arcsin(2(q_0q_2-q_3q_1)) \\ atan2(2(q_0q_3+q_1q_2),1-2(q^2_2+q^2_3)) \end{bmatrix}$$

and in our case:

q0 = cos(theta/2)

q1 = sin(theta/2) * nx

q2 = sin(theta/2) * ny

q3 = sin(theta/2) * nz