I think I wasn’t clear last time. What I was saying was that you should construct your own matrix. You are carrying out the multiplications for the single point. What happens if you want to rotate many points? Constructing a matrix like the one that glRotate does and then using glMultMatrixd or glLoadMatrixd is a more general way to handle the rotation.
However, if you still want to do it yourself, here is a way to think about it. The matrices that define rotations about the x, y and z axes are well known (look in the redbook or any linear algebra book). They are called the elementary rotations and accomplish a rotation about the x, y or z axis by an arbitrary angle, say theta.
What you want to do is rotation about an arbitrary vector by an arbitrary angle. This problem is only slightly harder than the elementary rotation problem. You just need to decompose the problem into one you already know (specifically, the elementary rotation problem).
Suppose the arbitrary axis that you want to rotate on is denote as: V
One way to achieve the rotation of theta about V is to do find a rotation (we’ll call it R) that rotates V into the elementary Z axis. Apply the elementary rotation that rotates an angle theta about the Z axis (this is the easy problem). Then apply the inverse of R to rotate the vector V back to it’s original direction. Ez is elementary rotation matrix of theta about the z axis, and R’ is the inverse of R, then the operations that you will apply to define the final matrix are:
(R’) (Ez) ® = C (where C is the matrix that you are looking for).
This leaves the problem of how to find that pesky R matrix (the inverse of it is easy). For this, you just need to define a coordinate system where V is the Z-axis. You’ll need to define a local X and Y axis for this local coordinate system, so pick any old vectors that make the cross products. Also, make sure that they are all unit vectors. What you can do is slap the local X, Y and your unitized V vector into a 4x4 matrix as columns like this:
X0 Y0 V0 0
X1 Y1 V1 0
X2 Y2 V2 0
0 0 0 1
This matrix (or it’s transpose, depending on how you think of the problem) is your R matrix. However you are think of the problem, the transpose of your R matrix is also the inverse by virtue of correctly picking your local X and Y vectors and making sure that all three (local X, local Y and V) are unit vectors.
Well, this is just one way to think of the problem. The are, for sure, others. Just remember that anytime you take three orthonormal vectors and slap them into a matrix as I did above, you are defining a transformation from a local coordinate system (whose origin happens to be 0,0,0) to the global coordinate system. This helps enable you to decompose a hard problem into simpler problems that involve simple transformations from one coordinate system to another.
Hope this helps.