psSphereTransformApply()

[gusciora@…]

I have included the code for the sphere transform below. This originated from a
personal email with Gene. I don't see this formula in the ADD. Can you verify
that this is correct?

x = coord->r;
y = coord->d;
dx = x - transform->phiP;
sinY = cos(y) * sin(dx) * transform->sinDeltaP +

sin(y) * transform->cosDeltaP;

cosY = sqrt(1.0 - sinY * sinY);
sinX = (cos(y) * sin(dx) * transform->cosDeltaP -

sin(y) * transform->sinDeltaP) / cos(y);

cosX = cos(y) * cos(dx) / cos(y);

out->r = atan2(sinX, cosX) + transform->alphaP;
out->d = atan2(sinY, cosY);

[Paul Price]

That's essentially equations 55-57 in the ADD. However, it uses different
signs. I know that the signs used in the equations in the ADD work with the
values for the Euler angles (i.e., alpha_p, delta_p, phi_p) for the conversions
(e.g., equatorial to Galactic) in the ADD. Your code (by which I mean the
choice of signs) is correct, but is relevant for a different specification of
the Euler angles than in the ADD.

[gusciora@…]

Gene's code differes from the ADD as far as signs are concerned, which you
mentioned. It also differs in the use of the sqrt() and atan2() function.
I have included my best guess at the current ADD algorithm. Is is correct? If
so, it would help to have the ADD include the calculation of the final (phi,
theta) pairs: this is one point where the ADD and Gene's code differ.

psF64 alpha = coord->r;
psF64 delta = coord->d;
psF64 alphaMinusAlphaP = alpha - transform->alphaP;

psF64 eq55 = (sin(delta) * transform->cosDeltaP) -

(cos(delta) * transform->sinDeltaP * sin(alphaMinusAlphaP));

psF64 eq56 = (cos(delta) * transform->cosDeltaP * sin(alphaMinusAlphaP)) +

(sin(delta) * transform->sinDeltaP);

psF64 eq57 = cos(delta) * cos(alphaMinusAlphaP);

psF64 theta = arcsin(eq55);
psF64 phi = atan2(eq56, eq57) + transform->alphaP;

out->r = phi;
out->d = theta;

[Paul Price]

Your code looks OK, though variable names "eq55" etc aren't very meaningful
mathematically.

The ADD specifies the algorithm quite clearly:

Since $\theta$ and $\delta$ have domains of $-\pi/2, \pi/2$, the value
of these angles are found by applying the arcsin to the sine of these
angles ($\theta = \arcsin \sin \theta$) which is always single-valued
and defined. The value of $\alpha-\alpha_p$ may be found from
\code{atan2(y,x)}, where $y = \cos \delta \sin (\alpha - \alpha_p)$
and $x = \cos \delta \cos (\alpha - \alpha_p)$; and similarly for
$\phi-\phi_p$.

[gusciora@…]

Okay.

My personal preference is that these formulas be specified explicitly, rather
than us having to derive it. I think we will be less likely to make mistakes if
that were to be, and the initial coding would be faster, and it would easier for
our code reviewers to verify that the function was coded correctly.