Euler angles don't have convention specified

[Paul Price]

Rotations, as well as the quaternion representation of a rotation, depends on
the convention adopted, and there is no clear information in the SDRS or ADD
that directly specifies what convention we adopt.

Figure 1 on p23 of the ADD (PSDC-430-006-011) seems to specify a left-handed
system with a zyz convention, but most of the information available on the web
are for right-handed systems. For example, see
​http://www.euclideanspace.com/maths/geometry/rotations/euler/index.htm , where I
think the "NASA Standard Aerospace" is a zyz convention in a right-handed
system, but the matrices don't match the "y convention" at Mathworld
​http://mathworld.wolfram.com/EulerAngles.html

Since this is an area abounding in confusion, we need to be be very clear about
the conventions we adopt and specify all relevant formulae (e.g., it is not
immediately clear how to convert the Euler angles fed into psSphereRotAlloc into
the quaternion representation). We also need to be careful that we are
supplying the correct formulae, since there are so many.

Perhaps Ed could provide the appropriate formulae?

[Paul Price]

OK, after much banging of my head, I think I've got this worked out. The ADD is
not clear on some things; I will mention these here, and then clean them up in
the document. Here's what I've got:

1. To Convert a position (alpha,delta) to a quaternion:

p0 = cos(alpha) cos(delta)
p1 = sin(alpha) cos(delta)
p2 = sin(delta)
p3 = 0

2. The appropriate Euler angles, in degrees, are:

(a) ICRS to Galactic:

alphaP = 180 - 192.85948
deltaP = 90 - 27.12825
phiP = 90 + 32.93192

(b) ICRS to Ecliptic:

alphaP = 270;
deltaP = 23o27'8.26" - 46.845" T - 0.0059" T^{2 + 0.00181" T}3
phiP = 90;

where T is the number of Julian centuries from 1900
(c) Precession

alphaP = 180 + 0.6406161 T + 0.0000839 T^{2 + 0.0000050 T}3
deltaP = 0.5567530 T - 0.0001185 T^{2 - 0.0000116 T}3
phiP = 180 + 0.6406161 T + 0.0003041 T^{2 + 0.0000051 T}3

where T is the number of Julian centuries between the two dates of interest.
These values replace those given in 3.4.1, 3.4.2 and 3.4.3.

3. Quaternion multiplication in the ADD is correct ("Combining two rotations" in

the ADD, section 3.3.2)

4. The ADD is incorrect about the use of quaternions to do rotation. Given a

rotation quaternion, r, and a position, p, then the resultant position is:

q = r p rbar

where rbar is the quaternion conjugate. Quaternion multiplication is
associative (whether you do the left pair or the right pair first doesn't
matter) but not commutative (you can't switch the order of the operands).

5. Conjugation: if a quaterion is defined,

q = (q0, q1, q2, q3),

where (q0, q1, q2) comprise the vector component, and q3 the scalar component,
then the quaternion conjugate is qbar = (-q0, -q1, -q2, q3).

6. The rotation is performed successively around the z, y, z axes

(a) Z axis: the rotation quaternion is:

r0 = 0
r1 = 0
r2 = sin(alphaP/2)
r3 = cos(alphaP/2)

(b) Y axis: the rotation quaternion is:

r0 = 0
r1 = sin(deltaP/2)
r2 = 0
r3 = cos(deltaP/2)

(c) Z axis again: the rotation quaternion is:

r0 = 0
r1 = 0
r2 = sin(phiP/2)
r3 = cos(phiP/2)

Each of these rotations is applied in turn, according to the rule in #4 above.

7. The resultant position is calculated from the quaternion, (q0, q1, q2, q3),

as follows:

phi = atan2(q1, q0)
theta = asin(q2)

where phi is the longitude and theta is the latitude.

8. I haven't checked 3.3.4, but 3.3.5 seems to work (though not extensively tested).

I will attach demonstration Perl code.

[Paul Price]

Demonstration perl code for rotations ICRS --> Galactic, Ecliptic and simple precession

[Paul Price]

Updated ADD entry

[Paul Price]

I think this is done. If there are any problems now, please let me know.

[David.Robbins@…]

I've been working on the SphereRot functions, trying to get just the Alloc
verified. However, based on the formulas given, the values do NOT match the
expected.
I'm not sure at this point whether the problem is the calculation of the
quaternions by the proposed Alloc method or from the verifying method. Note,
the alloc method uses the q = r * p * rbar method. The verifying method uses
n-hat * sin(angle/2) formula described by EulerAngles at mathworld and step one
of comment 1 above. These values should match but don't.
A couple of notes on the described methods. In step 5 above, it is noted that
the conjugate of a quaternion should be (-q0, -q1, -q2, q3). For whatever
reason, the rotation quaternion for rbar does not match this formula. If it
did, it would have r2= -sin(alphaP/2) as opposed to sin(phiP/2) and
r3=cos(alphaP/2) as opposed to cos(phiP/2). Also, the ADD appears unclear (to
me) about the order of combination of these three quaternions. It says first
that the formula is r*p*rbar so r*p should be done first. Afterwards it gives
the order as tsr (which would correspond to rbar*p*r). While the operation is
associative, it is Not commutative, so it seems to me that there is some
discrepancy here. I would assume the first way is correct.

[David.Robbins@…]

I think I figured one thing out. I'll run it by you tomorrow.

[David.Robbins@…]

While most of my comments from 5 above are still valid, I have found some useful
information from mathworld. For a given rotation, a quaternion can be expressed
in terms of the angles (alpha, delta, phi) as:
q0 = sin[1/2(alpha-phi)] * sin(1/2*delta)
q1 = cos[1/2(alpha-phi)] * sin(1/2*delta)
q2 = sin[1/2(alpha+phi)] * cos(1/2*delta)
q3 = cos[1/2(alpha+phi)] * cos(1/2*delta)

Implemented correctly, these values match those from the method given in ADD for
rotating a vector by combination of quaternions (p26). I'm still not exactly
sure why r-bar is not equal to the conjugate of r here, but this seems to be
correct.

[Paul Price]

Bear in mind that Mathworld define their quaternions differently (i.e., they
swap q0 and q3, or similar).

There are test cases in the ADD that should indicate if the rotations are being
done correctly.

[Paul Price]

How was this bug fixed?

[David.Robbins@…]

The bug was fixed by using the formulas from comment 7 to verify the quaternion
values resulting from SphereRotAlloc. Along with implementing the changes sent
to me concerning the inverse function and the conjugate function, all of the
values from every SphereRot test match the expected.