A) The easiest way, although with non-constant speed: Interpolate the positions directly in a linear fashion but re-scale the vector to the appropriate length. This works best if the start and end point are relatively close. This method fails miserably if start and end position are at most on opposites sides of the sphere.
With start position s and end position e, both in world co-ordinates, compute their counterparts relative to sphere center c as
s' := s - c
e' := e - c
Compute the difference vector from start to end
d' := e' - s'
Interpolate the difference vector in N steps to get the linear in-between positions
b'n := s' + d' * n / N w/ 0<= n <= N
Rescale the vector to a length equal to the radius r of the sphere
p'n := b'n / | b'n | * r
Eventually transform the current position back into world space
pn := p'n + c
B) Utilizing a rotation, getting constant speed: With the start and end position given as vectors relative to the sphere's origin, compute the axis of rotation as cross-product of the 2 vectors and the difference angle from the dot-product formula. Then interpolate the angle in as many steps as you want, and use it together with the axis to build a rotation matrix. Apply the rotation onto the start position vector. This method fails (as all do more or less) if start and end position are at most on opposites sides of the sphere.
With start position s and end position e, both in world co-ordinates, compute their counterparts relative to sphere center c as
s' := s - c
e' := e - c
Compute the axis of rotation r as normalized cross-product
r := ( s' x e' ) / | s' x e' |
Compute the angle of rotation a from the dot-product
a := arccos( ( s' . e' ) / ( | s' | * | e' | ) )
Interpolate the angle in N steps to get the in-between angle b
bn := a * n / N w/ 0<= n <= N
Compute the rotation matrix R from axis and in-between angle
Rn( r, bn )
Apply the rotation onto the start position to yield in the current position
p'n := Rn * s'
Eventually transform the current position back into world space
pn := p'n + c
Another way (as shown in the cited post in the OP) is to convert the positions into quaternions and using nlerp or slerp, but that is IMHO less clear from a mathematical point of view.