For two unit vectors a and b:
Code: Select all
ang = 2 * atan2 (mag(a - b), mag(a + b);// use 360/pi for degrees
This is perfectly safe at any angle, and can never produce a nan unless either a or b has a nan in a component.
If either a or b may be a non-unit vector (even zero!!!), then with a slight modification the angle can still be safely determined:
Code: Select all
ang = 2 * atan2(mag(mag(b)*a - mag(a)*b), mag(mag(b)*a + mag(a)*b));
There are only two disadvantages to using this: it's a little slower (atan and thus atan2 (which is just a very smart wrapper for atan) is inherently slower than acos or asin), and you lose the direction of the angle. However, this can be retrieved by finding the cross-product of a and b and dotting that with a reference direction to get the sign of the angle.
Some nice bonus information...
For unit vectors, mag(a - b) == 2*sin(ang(a,b)/2), ie, twice the sin of the half-angle between a and b, and mag(a + b) == 2*cos(ang(a,b)/2) (twice the cos of the half-angle)
Now, this is very handy when it comes to finding the quaternion describing the rotation between the two vectors. When I implemented the TINU (There Is No Up) mod for KSP, I discovered that Unity's Quaternion.FromToRotation was not usable for very small angles (think per-frame rotation while in low orbit, worse in high orbit). The reason for this is, I assume, due to using epsilons and guarding against going outside the -1..1 range for acos.
I'll leave the derivation as an exercise (hint, it takes advantage of trig identities), but this is taken directly from the code in TINU:
Code: Select all
Quaternion fromtorot(Vector3 a, Vector3 b)
{
float ma = a.magnitude;
float mb = b.magnitude;
Vector3 mb_a = mb * a;
Vector3 ma_b = ma * b;
float den = 2 * ma * mb;
float mba_mab = (mb_a + ma_b).magnitude;
float c = mba_mab / den;
Vector3 v = Vector3.Cross (a, b) / mba_mab;
return new Quaternion(v.x, v.y, v.z, c);
}