# SPHERICAL/ASTRONOMICAL COORDINATE TRANSFORMATION

Carl Heiles

In our never-ending attempt to make your life easier, we present you with the quickest of quick summaries of spherical coordinate transformation with matrix techniques. For example, you might need to point the telescope at some particular position in the Galaxy. In other words, you need to convert from Galactic coordinates to altitude and azimuth. This involves three separate coordinate transformations: Galactic longitude and latitude to equatorial right ascension and declination [$(\ell , b) \rightarrow (\alpha , \delta )$]; to hour angle and declination [$(\alpha , \delta ) \rightarrow (ha, \delta )$]; to azimuth and altitude [$(ha, \delta ) \rightarrow (az, alt)$]. Or maybe you need to go all or partway in the the other direction, i.e. [$(az, alt) \rightarrow [(\alpha , \delta )$] or maybe just [$(az, alt) \rightarrow [(ha, \delta )$].

These are conversions among four spherical coordinate systems. Such conversions involve all those complicated combinations of trig functions (sigh). You can write down all these trig functions for the various possible transformation—twelve possibilities in all if you include going both directions.

But there is a much easier, more elegant, and more politically correct way: using rotation matrices. In this method, you generate a vector in the original coordinate system; convert the vector to another coordinate system by rotating the coordinates using matrix multiplication; and convert the vector to the angles of the new coordinate system.

There are two big advantages with this method. First, you can apply several transformations in succession by multiplying the rotation matrices in succession, so you break the process down into single transformations, each with its own rotation matrix. Second, it’s easy to go “backwards”—you just use the inverse of the matrix.

The method is general and can be applied to any coordinate transformation. Spherical coordinates are characterized by two angles. One “goes around the z-axis”—it is like longitude on the earth. The other “goes up and down” and is like latitude on the earth. These angles are “longitude-like” and “latitude-like”, and we’ll denote them by long and lat. Thus, for Galactic coordinates, $\ell$ is the “longitude-like” long and b is the “latitude-like” lat; for equatorial coordinates, it’s α (or ha) and δ; for terrestrial coordinates, it’s az and alt.

One more thing before we get into details. Our discussion is oriented towards astronomy, but the method works for any type of spherical coordinate transformation. There is an excellent, short discussion of the general situation in Goldstein’s Classical Mechanics, §4.4.

## 1 ROTATION MATRICES: THE METHOD

To restate the problem: we start with (long,lat) in one coordinate system and want to convert to $(long^\prime , lat^\prime )$ in some other coordinate system. Here’s the prescription:

Step 1. First, convert the angles to rectangular coordinates. One would usually call these (x,y,z); here, to emphasize the vector/matrix flavor, we call them (x0,x1,x2) and denote the 3-element vector ${\vec x}$. To accomplish this conversion in Python:

x = numpy.array([0.,0,0])
x[0] = numpy.cos(lat) * numpy.cos(long)
x[1] = numpy.cos(lat) * numpy.sin(long)
x[2] = nump.sin(lat)

Step 2. Apply the rotation matrix R (we’ll discuss its definitions below):

${\vec x^\prime } = {R \cdot \vec x} \; . \,\!$

In Python, we’ll use numpy arrays to represent $\vec x$ and ${\vec x^\prime }$, and you do matrix multiplication with numpy.dot:

xp = numpy.dot(R, x)

Step 3. Convert the primed rectangular coordinates to the new set of spherical coordinates $(long^\prime , lat^\prime )$:

\begin{align} long^\prime & = \tan ^{-1} \left({{\vec x^\prime }_1 \over {\vec x^\prime }_0}\right) \; , \\ lat^\prime & = \sin ^{-1} ({\vec x^\prime }_2) \; . \end{align}\,\!

To do these in Python, where we’ll write longp, latp:

longp = numpy.arctan2(xp[1], xp[0])
latp = numpy.arcsin(xp[2])

That’s it! If you want to “go backwards”, you just apply the matrix multiplications using the inverse matrices. For rotation matrices, the inverse is always equal to the transpose(!)If you don’t believe this, check on the R’s below in §2!—symbolically for the matrix R, R − 1 = RT. So to go from the primed system to the unprimed, you switch the primed and unprimed in equation (2) and use the inverse rotation matrix

${\vec x} = {R^{-1} \cdot x^\prime } = {R^ T \cdot x^\prime }\; . \,\!$

and in Python

iR = numpy.linalg.inv(R)
x = numpy.dot(iR, xp)

## 2 ROTATION MATRICES: SPECIFICS FOR OUR PROBLEM

OK, what are these rotation matrices? We’ll do you a big favor and tell you.

### 2.1 (RA, DEC) to (HA, DEC)—any epoch.

Converting from $(\alpha , \delta ) \rightarrow (ha, \delta )$ keeps the declination the same and uses the relationship ha = LST − α. It is easiest to think of this in two steps, so we express

${R}_{(\alpha , \delta ) \rightarrow (ha, \delta )} = {R}_{(\alpha , \delta ) \rightarrow (ha, \delta ), 2} \cdot {R}_{(\alpha , \delta ) \rightarrow (ha, \delta ), 1} \; . \,\!$

First, we rotate around the equatorial pole by an angle equal to the Local Sidereal Time (LST), which does $\alpha \rightarrow (\alpha - LST)$:

\begin{align} {R}_{(\alpha , \delta ) \rightarrow (ha, \delta ),1} = \left[ \begin{array}{ccc} \cos (LST) & \sin (LST) & 0 \\ -\sin (LST) & \cos (LST) & 0 \\ 0 & 0 & 1 \\ \end{array} \; \right] \; . \end{align}\,\!

Next, the ha and α go in opposite directions, which is equivalent to converting from the original left-handed to a right-handed coordinate system, so the second step is just to perform this reversal:

\begin{align} {R}_{(\alpha , \delta ) \rightarrow (ha, \delta ), 2} = \left[ \begin{array}{ccc} 1 & 0 & 0 \\ 0 & -1& 0 \\ 0 & 0& 1 \\ \end{array} \; \right] \; . \end{align}\,\!

The full rotation matrix is the matrix product ${ R}_{(\alpha , \delta ) \rightarrow (ha, \delta ),2} { \cdot } { R}_{(\alpha , \delta ) \rightarrow (ha, \delta ), 1}$. Note the order! Applying ${R}_{(\alpha , \delta ) \rightarrow (ha, \delta ),2}$ at the beginning in the matrix product means that it operates last on the vector ${\vec x}$, which is what we want. So we have as the product

\begin{align} {R}_{(\alpha , \delta ) \rightarrow (ha, \delta )} = \left[ \begin{array}{ccc} \cos (LST) & \sin (LST) & 0 \\ \sin (LST) & -\cos (LST) & 0 \\ 0 & 0 & 1 \\ \end{array} \; \right] \; . \end{align}\,\!

### 2.2 (HA, DEC) to (AZIMUTH, ALTITUDE).

Because (ha,δ) are Earth-based coordinates, this conversion depends only on your terrestrial latitude φ:

\begin{align} {R}_{(ha, \delta ) \rightarrow (az, alt)} = \left[ \begin{array}{ccc} -\sin \phi & 0 & \cos \phi \\ 0 & -1 & 0 \\ \cos \phi & 0 & \sin \phi \\ \end{array} \; \right] \; . \end{align}\,\!

### 2.3 EQUATORIAL to GALACTIC.

\begin{align} {R}_{(\alpha , \delta )_{1950} \rightarrow (\ell , b)} = \left[ \begin{array}{rrr} -0.066989 & -0.872756 & -0.483539 \\ 0.492728 & -0.450347 & 0.744585 \\ -0.867601 & -0.188375 & 0.460200 \\ \end{array} \; \right] \; . \end{align}\,\!

This is from Green’s Spherical Astronomy, chapter 14.6, problem 14.6 (answers in back of book). We should’ve made you derive this, but we’re softies. You really should at least glance at Green’s chapter 2.7, which defines Galactic coordinates. In truth, the precession of the equatorial coordinate system makes this matrix a function of time: the equatorial coordinates move around the sky, but the Galactic ones do not. Precession amounts to nearly an arcminute per year! In principle, you should derive the matrix for the current epoch. In practice, you may not need such high accuracy. Better than the 1950 version are the numbers for epoch 2000, which are in Green’s equation (14.55); epoch 2000 is lots closer to the present than is epoch 1950:

\begin{align} {R}_{(\alpha , \delta )_{2000} \rightarrow (\ell , b)} = \left[ \begin{array}{rrr} -0.054876 & -0.873437 & -0.483835 \\ 0.494109 & -0.444830 & 0.746982 \\ -0.867666 & -0.198076 & 0.455984 \\ \end{array} \; \right] \; . \end{align}\,\!

### 2.4 Precession—converting equatorial between epochs.

We won’t need these for the lab course, but we give you the info for the sake of completeness. Generating the rotation matrix for precession is a bit tedious and we won’t give the explicit formulae here. They are in Green’s book. The elements of the matrix are in equation (9.31). These elements contain angles, which depend on time as in equation (9.23) if you are converting from epoch 2000 to some other epoch. Precession isn’t all there is; for precision exceeding ˜10 you also need to account for nutation of the Earth, which has a random component and is not completely predictable. For the complete story, see Green’s chapter 9 and The Astronomical Almanac 1998, pages B39-B43—for interested parties only!

## 3 DOING ALL THIS IN Python

Obviously, all this stuff is simple in Python/Numpy, which deals easily with matrices. Before beginning, though, a cautionary note about 2-D arrays on computers:

### 3.1 A CRUCIAL PRELIMINARY: 2-D arrays

In a computer, a multidimensional data set can be indexed in two ways, the column-major and row-major formats. IDL and Fortran use the row-major format; C and Python use column-major. Suppose you have a $2 \times 2$ matrix called A. In IDL’s row-major format, when you type [print, A] IDL prints

\begin{align} \left[ \begin{array}{rrr} A_{0,0} & A_{1,0} \\ A_{0,1} & A_{1,1} \\ \end{array} \; \right] \; , \end{align}\,\!

which is different from what you are used to seeing in standard matrix notation (and Python and C) which is the column-major format

\begin{align} \left[ \begin{array}{rrr} A_{0,0} & A_{0,1} \\ A_{1,0} & A_{1,1} \\ \end{array} \; \right] \; . \end{align}\,\!

Thank goodness we are using Python! Why does IDL do this nonstandard thing? It’s because it’s more straightforward for image processing, in which traditionally the images are scanned row-by-row (as in a TV set) instead of column-by-column. And IDL’s origins are image processing, not matrix math.

### 3.2 Try the following examples in Python.

Test ${R}_{(ha, \delta ) \rightarrow (az, alt)}$, going both forwards and backwards. For an observatory at latitude $41.36^\circ$, $(az,alt)=(137.60^\circ ,32.43^\circ )$ transforms to $(ha,\delta )=(325.05s^\circ , -6.52^\circ )$. Hour angle is usually given in hours using sexagesimal notation: ha = 21h40m12s.) See the PyEphem module for functions that go back and forth between decimal and sexagesimal notations.

Test ${R}_{(\alpha , \delta ) \rightarrow (ha, \delta )}$ by making up your own example, using the fact that ha = LST − α.

Test ${R}_{(\alpha , \delta )_{1950} \rightarrow (\ell , b)}$ for the Crab Nebula. The Crab has 1950 equatorial coordinates $(\alpha , \delta )=(05^ h31^ m.5, 21^\circ 59’)$ and Galactic coordinates $(\ell ,b)=(184^\circ 33’, -5^\circ 47’)$.

Finally, put them all together and make sure that works, too. For all of these, make sure you know how to go backwards! For example, suppose you want to convert $(az, alt) \rightarrow (\alpha , \delta )$. You need to first apply ${R}_{(ha, \delta ) \rightarrow (az, alt)}^{-1}$ and then ${R}_{(\alpha , \delta ) \rightarrow (ha, \delta )}^{-1}$. So the full rotation matrix in this case is

${R}_{(az, alt) \rightarrow (\alpha , \delta )} = {R}_{(\alpha , \delta ) \rightarrow (ha, \delta )}^{-1} \cdot {R}_{(ha, \delta ) \rightarrow (az, alt)}^{-1} \,\!$

Again, note the order! Applying ${R}_{(\alpha , \delta ) \rightarrow (ha, \delta )}^{-1}$ at the beginning in the matrix product means that it operates last on the vector ${\vec x}$.