Dealing with singular points in equatorial coordinates

[javierggt]

The following test (which explores an occurrence of gimbal lock) currently fails:

q1 = Quat([90, 90, 0])
q2 = Quat([0, 90, 0])
assert np.allclose(q1.q, q2.q)

This is not a problematic issue, because Chandra never operates close to declination 90 or -90. Still, this is wrong and kind of bugs me. Let's see if we can fix it.

Tom added that the issue is worse, because the numerical precision worsens as declination approaches 90 or -90. What would be the point of fixing this issue if we still have floating point precision issues? It remains to be seen whether we can fix this.

A possible solution to the first issue

This is my current understanding of the issue. The rotation matrix in equatorial coordinates is given by:

[cos(alpha)*cos(delta), -sin(alpha)*cos(rho) - sin(delta)*sin(rho)*cos(alpha),  sin(alpha)*sin(rho) - sin(delta)*cos(alpha)*cos(rho)],
[sin(alpha)*cos(delta), -sin(alpha)*sin(delta)*sin(rho) + cos(alpha)*cos(rho), -sin(alpha)*sin(delta)*cos(rho) - sin(rho)*cos(alpha)],
[           sin(delta),                                   sin(rho)*cos(delta),                                   cos(delta)*cos(rho)]

where (alpha, delta, rho) is (RA, dec, roll).

alpha and rho are currently determined by doing:

alpha = arctan2(m[0,1], m[0,0])
rho = arctan2(m[2,2], m[2,1])

and this fails when m[0,1] = m[0,0] = m[0,1] = m[0,0] = 0). One could do instead (by convention, choose rho=0 to break the degeneracy):

• if m[0,2] == 1:

  rho = 0
  alpha = arctan2(m[0,1]+m[1,2], m[0,2]-m[1,1])

• if m[0,2] == -1:

  rho = 0
  alpha = arctan2(m[0,1]-m[1,2], m[0,2]+m[1,1])