Plotting Equatorial coordinates to X,Y plane (simulating telescope/camera view)

Question

I would like to simulate a telescope camera view, projecting (RA, Dec) coordinate sets to a 2D plane. So effectively the need is to convert Equatorial coordinates to (x, y) coordinates, when we can know how many degrees a single pixel spans. For convenience I'm using degrees for RA and Dec here, and the code is just pseudocode.

The following base information is provided as an example: Imaging equipment:

• Telescope focal length in mm (flen): 420
• Image resolution (widthPx, heightPx): 4000 * 3000
• Pixel size in microns (pxSize): 8
• Binning (bin): 1

So from my understanding we can calculate a few things from this for starters. The field of view in degrees should be:

fovW = radToDeg( 2 * atan((pxSize * bin * widthPx / 1000) / (2 * flen)) )
fovH = radToDeg( 2 * atan((pxSize * bin * heightPx / 1000) / (2 * flen)) )

And from this we could derive:

pixelsPerDegWidth = widthPx / fovW;
pixelsPerDegHeight = heightPx / fovH;

For locations closer to the equator I can get some valid looking results using this formula:

deltaRa = starRa - imageCenterRa
deltaRa = deltaRa > 180 ? 360 + deltaRa : deltaRa
deltaDec = -1 * (starDec - imageCenterDec)

posX = deltaRa * cosdeg(starDec) * pixelsPerDegWidth + 0.5 * widthPx
posY = deltaDec * pixelsPerDegHeight + 0.5 * heightPx

With some catalog star positions projected using this formula, I can even get Astrometry.net solver to correctly identify the center position. However when we move to the poles this comes crashing down.

So let's go ahead and use a target center point near the pole:

• RA 30 degrees (imageCenterRa)
• Dec 89 degrees (imageCenterDec)

Let's also say that we have a hypothetical star at RA, Dec (220, 89.2) - this is close to the pole, and both the star and the pole gets included in the image. What kind of formula should we use to get the correct location and to draw the area correctly even in the polar regions?

I've been looking for an example for a while but haven't come up with any. I've seen some formulas that may be just the thing I'm looking for but my thick skull loves examples - hence why I included some exact numbers here, as I would love to see how this would be calculated in practice.

So to summarize: how would we calculate the (X, Y) coordinate in the image for the star at RA, Dec (220, 89.2)?

Answer 1

This isn’t code for a full implementation, but you could do something like this:

1. Transform your RA, Dec coordinates (which are spherical polar coordinates, just all with the same $r$) to Cartesian $x, y, z$ coordinates. (The $r$ value you would use here is arbitrary, for example $r = 1$, and should be the same for all stars, not their real distances; you're just projecting them onto a sphere so you can do the transformation.)

2. Apply a coordinate frame rotation using your field-center RA, Dec coordinates (or perhaps their complements) as the rotation angles, so that the field center direction is the new $z$ axis.

3. Plot the resulting new $x, y$ coordinates, ignoring $z$.

You might have to do a further rotation of those $x, y$ coordinates to get the plot oriented the way you want, or maybe there’s a clever way to do that with the proper choice of angles in step 2.

Answer 2

Ok, it needed some searching with the right words (astrometry, standard coordinates) to actually find the equations I needed. I'll try to go through the process here to make this a complete answer to my own question.

The key equation here is the conversion of Celestial Coordinates (RA, Dec) to Standard Coordinates (X, Y) on a tangent plane of the sky. This figure from the article Astrometry: The Foundation for Observational Astronomy (Amith Govind et al) explains what is needed.

From the same article we can find the equations for calculating the Standard Coordinates:

$$X = \frac{cos(\delta)sin(\alpha - \alpha_0)}{cos(\delta_0)cos(\delta)cos(\alpha - \alpha_0) + sin(\delta)sin(\delta_0)}$$ $$Y = \frac{sin(\delta_0)cos(\delta)cos(\alpha - \alpha_0) - cos(\delta_0)sin(\delta)}{cos(\delta_0)cos(\delta)cos(\alpha - \alpha_0) + sin(\delta)sin(\delta_0)}$$

where $\alpha$ is the angle RA of the star and $\delta$ is the angle Dec of the star, and $\alpha_0$ is the field center RA and $\delta_0$ the field center Dec, all in radians.

With these equations we get the Standard Coordinates in radians in respect to the center of the field. To get the pixel coordinates, those coordinates will still need to be multiplied by the amount of pixels per radians and those rations can be calculated from the focal length, binning, image resolution and the camera sensor pixel size. If we want to account for a possible field rotation, we will need to factor in the rotation of the coordinates by forming the coordinates from X and Y components with sine and cosine factors. In addition, as the above illustration shows the X coordinate will be mirrored, so we'll need to flip the final resulting X coordinate.

So proceeding on with pseudocode, we have:

imageWidthRad = 2 * atan((pxSize * bin * widthPx / 1000.0) / (2 * fLen))
imageHeightRad = 2 * atan((pxSize * bin * heightPx / 1000.0) / (2 * fLen))

pixelsPerRadW = widthPx / imageWidthRad
pixelsPerRadH = heightPx / imageHeightRad

imageCenterX = widthPx / 2
imageCenterY = heightPx / 2

(note that pixel size (in microns) and binning is assumed to be the same for both width and height here)

With star RA, Dec (s_ra, s_dec) and field center RA, Dec (c_ra, c_dec) in radians we get the standard coordinates:

stdX = cos(s_dec) * sin(s_ra - c_ra) / (cos(c_dec) * cos(s_dec) * cos(s_ra - c_ra) + sin(c_dec) * sin(s_dec))
stdY = (sin(c_dec) * cos(s_dec) * cos(s_ra - c_ra) - cos(c_dec) * sin(s_dec)) / (cos(c_dec) * cos(s_dec) * cos(s_ra - c_ra) + sin(c_dec) * sin(s_dec))

If rotation is assumed as zero, we now have everything we need to calculate the star positions:

starPixelX = pixelsPerRadW * stdX + imageCenterX
starPixelY = pixelsPerRadH * stdY + imageCenterY
starPixelX = width - starPixelX

To account for field rotation, we need to change it a little. Let's add the field rotation angle $r$.

starPixelX =  cos(r) * pixelsPerRadW * stdX + sin(r) * pixelsPerRadH * stdY + imageCenterX
starPixelY = -sin(r) * pixelsPerRadW * stdX + cos(r) * pixelsPerRadH * stdY + imageCenterY
starPixelX = width - starPixelX

With this, the star positions get rotated accordingly. At 90 degrees (cos(90deg) == 0, sin(90deg) == 1) we see that X effectively becomes Y and vice versa.

Conclusion

With a program that reads the RA, Dec coordinates of Tycho-2 catalog stars I drew the stars to an image using these equations, using the following parameters:

• Pixel size (pxSize) = 3.8 microns
• Image Width (widthPx) = 4656
• Image Height (heightPx) = 3520
• Center RA (c_ra) = 194.464 deg
• Center Dec (c_dec) = 71.217 deg
• Binning (bin) = 1
• Focal length (fLen) = 25mm
• Rotation (r) = -268 deg

Resulting in this image:

I then ran this through http://nova.astrometry.net astrometric solver, and successfully solved the image with the calculated (RA, Dec) center being: (194.461, 71.223).

As astrometry.net annotated with grid:

So there it is, projected stars from sphere to a flat plane.