So one day I decided that I wanted to write my own astronomy program to make sky charts and to plot deep sky objects. The process seemed simple enough: take the celestial coordinates of whatever object interests you (say, the star Deneb in the constellation Cygnus), convert it to x/y pixels for your screen and then plot it! And, actually, it pretty much is that straight forward.
Just about any development language for a computer will provide you some means of controlling the pixels on your screen. Most window-based environments (MS Windows, Mac, Linux running a GUI) provide libraries or APIs to facilitate this. For this example, I am working in Windows, and I will use the Graphics Display Interface (GDI) to do my drawing. To get access to this API, you just need to include the Windows header file in your source code - you will have done this already since you are likely building a Windows application! The bottom line is if you can set a pixel in your language of choice, you can build an astronomy program!
The star Deneb has the following coordinates in the sky (J2000):
Let's assume we have a window with a client area 800 pixels wide by 600 pixels high. How are we going to figure out which pixel should represent Deneb? We need to know a few more things, besides the star's coordinates. We need to know the RA and Dec that our field of view is centered on, and we need to know how wide (in degrees) our field of view is. If we have this information, we can plot our star!
We'll use a custom coordinates structure to return our x/y values in. Perhaps something like:
struct COORDS {
public int x;
public int y;
};
Now, we can build a function that will return a populated data structure:
static COORDS EQConvert(double RA, double Dec) {
/*
* NOTE: This function accepts the RA and Dec of an object (in degrees), and
* then returns the appropriate x/y pixels based on the current display.
*/
/* some handy numbers that could be constants elsewhere in the program */
double PI = 3.1415926535897932384626433832795; /* reasonable value for pi */
double TwoPI = 6.283185307179586476925286766559; /* double pi! */
double PI2 = 1.5707963267948966192313216916398; /* half pi */
double Rads = 0.01745329251994329576923690768489; /* radians to degrees and back */
/* these variables will contain our display's parameters - they would normally be
retrieved from somewhere else like a global struct variable or class */
double displayRA = 300.0; /* Right Ascension for display in degrees */
double displayDec = 40.0; /* Declination for display in degrees */
double displayFOV = 60.0; /* Field of View for the display in degrees */
double displayROT = 0.0; /* Field Rotation in degrees */
int screenWidth = 800; /* Client Area width in pixels */
int screenHeight = 600; /* Client Area height in pixels */
COORDS myCoords; /* Coordinates structure so we can return x/y */
/* these are local working variables for this function */
double Az, Alt, tmpval, nz, az2, tx, ty, xyscale;
double fieldRA, fieldDec, fieldFOV, fieldROT;
/* first, convert all our values to radians (from degrees) */
RA = RA * Rads;
Dec = Dec * Rads;
fieldRA = displayRA * Rads;
fieldDec = displayDec * Rads;
fieldFOV = displayFOV * Rads;
fieldROT = displayROT * Rads;
/* second, let's convert our RA/Dec coordinates to Alt/Az coordinates */
Az = Math.Atan2(Math.Sin(fieldRA - RA), Math.Cos(fieldRA - RA) * Math.Sin(fieldDec)
- Math.Tan(Dec) * Math.Cos(fieldDec));
tmpval = Math.Sin(fieldDec) * Math.Sin(Dec) + Math.Cos(fieldDec) * Math.Cos(Dec)
* Math.Cos(fieldRA - RA);
if(tmpval >= 1.0) {
Alt = PI2;
} else {
Alt = Math.Asin(tmpval);
}
/* third, we can convert from Alt/Az to our x/y pixels */
nz = 1.0 - 2.0 * Alt / PI;
az2 = Az - PI2 + fieldROT;
tx = (nz * Math.Cos(az2)) * PI / fieldFOV;
ty = -(nz * Math.Sin(az2)) * PI / fieldFOV;
xyscale = ((double)screenWidth / fieldFOV) / (120.0 / displayFOV);
myCoords.x = (int)(((double)screenWidth / 2.0) + (tx * xyscale));
myCoords.y = (int)(((double)screenHeight / 2.0) + (ty * xyscale));
return(myCoords);
}
The math for the EQConvert() function was drawn from Astronomical Algorithms by Jean Meeus, as well as help from Philippe Deverchere (author of C2A - Computer Aided Astronomy) - many thanks Philippe! Now that we can convert our coordinates, we can call our function from elsewhere in our program (perhaps the WM_PAINT area?)...
void somefunction() {
COORDS skyCoords;
/* convert coordinates for star Deneb in constellation Cygnus */
skyCoords = EQConvert(310.3579167, 45.28027778);
/* now turn on the pixel! */
SetPixel(hDc, skyCoords.x, skyCoords.y, RGB(255, 255, 255));
/* etc... */
}
With a little reverse engineering, you can build the opposite function, which will take x/y coordinates (say, from a mouse click?) and convert them back to RA/Dec. That is what I did so that the center of the field of view on the display changes when you click the mouse.