sot / Quaternion

Quaternion manipulation
https://sot.github.io/Quaternion
BSD 3-Clause "New" or "Revised" License
3 stars 7 forks source link

Warning for equatorial calculations above 89.99999 #10

Open taldcroft opened 5 years ago

taldcroft commented 5 years ago

We should include warnings for any attempt to initialize a quaternion, or compute equatorial, for abs(dec) > 89.99999. Below that it seems to be numerically stable.

The transforms to and from equatorial are numerically unstable for abs(dec) very near 90:

>>> q_bad = normalize([0.24184476, -0.66446302,  0.24184476,  0.66446302])
>>> q4 = Quat(q=q_bad)
>>> q4.equatorial
array([  0.        ,  89.99999915,   0.        ])

>>> q5 = Quat(equatorial=q4.equatorial)
>>> q5.q
array([ 0.        , -0.70710678,  0.        ,  0.70710679])

Teasing out the exact numerical problem is a little tricky, but for the q4.equatorial statement I include debug print statements in _quat2equatorial and got:

        xa = q2[0] - q2[1] - q2[2] + q2[3]
        xb = 2 * (q[0] * q[1] + q[2] * q[3])
        xn = 2 * (q[0] * q[2] - q[1] * q[3])
        yn = 2 * (q[1] * q[2] + q[0] * q[3])
        zn = q2[3] + q2[2] - q2[0] - q2[1]
...
        print(xn, yn, zn, one_minus_xn2)                                                                            
        print(xa, xb)                                                                                               
        ra = degrees(atan2(xb, xa))
        dec = degrees(atan2(xn, sqrt(one_minus_xn2)))
        roll = degrees(atan2(yn, zn))

[ 1.] [ 0.] [ 0.] [  2.22044605e-16]
[ 0.] [ 0.]

So basically ra = atan2(0, 0) and roll = atan2(0, 0).

jeanconn commented 4 years ago

Do we still want/need to do this? I'm concerned it has been so long we will misremember and think there's a warning in 3.5.0 .