/* goto.c -- crude program for doing two-star alignment */ /* */ /* Notes: This program must be compiled with a math */ /* library. It currently only accepts decimal */ /* coordinates, and I have no immediate plans to allow */ /* sexagesimal input. It does, however, output in */ /* both decimal and sexagesimal. */ /* */ /* Update (1 Sep 1999): Updated code with new code */ /* based on dot (inner) products for angular */ /* separation, rather than Euclidean distance. The */ /* resulting algorithm isn't substantially different, */ /* but it is cleaner and easier to follow. */ /* */ /* The code also uses a new time correction code. A */ /* star that is initially pointed to may drift out of */ /* the field of view because the scope is not aligned */ /* with the celestial pole. There is a clever way to */ /* handle this by creating an imaginary star that is */ /* x hours to the west in RA, and determining where */ /* that imaginary star would have been x hours ago. */ /* This correctly locates the target star x hours */ /* after the initial two-star alignment. This new */ /* correction uses gettimeofday and might not port. */ /* */ /* Update (23 Feb 2000): Re-ordered alignment stars */ /* by right ascension. This makes it easier to pick */ /* appropriate stars. Ideally, the stars should be */ /* 90 degrees apart in the sky. This will generally */ /* be *more* than 90 degrees of RA, but with the */ /* exception of the circumpolar stars, choosing stars */ /* that are roughly 90 degrees (6 hours) of RA apart */ /* should be just fine. */ /* */ /* I also added the ability to use two user-defined */ /* stars. This should be useful if you have the */ /* precise positions of alignment stars that aren't */ /* on the prepared list. */ /* */ /* E-mail me should you want an explanation of this */ /* code. Alternatively, I plan to write an article */ /* on mounts and two-star alignments for my July 2000 */ /* Astronomical Games column, and you can wait for */ /* that. */ /* */ /* Copyright (c) 1998, 1999, 2000 Brian Tung */ #include #include #include #include #ifndef PI #define PI (3.1415926535) #endif #define DEGREES (PI/180.0) typedef struct { char name[80]; double ra; double dec; } stardata; stardata table[] = { {"alpha Eridani (Achernar)" , 24.428, -57.237}, {"alpha Ursae Minoris (Polaris)" , 37.946, +89.264}, {"beta Persei (Algol)" , 47.042, +40.956}, {"eta Tauri (Alcyone)" , 56.871, +24.105}, {"alpha Tauri (Aldebaran)" , 68.980, +16.510}, {"beta Orionis (Rigel)" , 78.634, -08.202}, {"alpha Aurigae (Capella)" , 79.172, +45.999}, {"alpha Orionis (Betelgeuse)" , 88.793, +07.407}, {"alpha Carinae (Canopus)" , 95.988, -52.696}, {"alpha Canis Majoris (Sirius)" , 101.289, -16.713}, {"epsilon Canis Majoris (Adhara)" , 104.656, -28.972}, {"alpha Geminorum (Castor)" , 113.650, +31.889}, {"alpha Canis Minoris (Procyon)" , 114.827, +05.228}, {"beta Geminorum (Pollux)" , 116.331, +28.026}, {"alpha Leonis (Regulus)" , 152.094, +11.967}, {"alpha Ursae Majoris (Dubhe)" , 165.933, +61.751}, {"beta Leonis (Denebola)" , 177.266, +14.572}, {"alpha Crucis (Acrux)" , 186.650, -63.099}, {"beta Crucis (Mimosa)" , 191.930, -59.689}, {"zeta Ursae Majoris (Mizar)" , 200.981, +54.925}, {"alpha Virginis (Spica)" , 201.298, -11.161}, {"beta Centauri (Hadar)" , 210.956, -60.373}, {"alpha Bootis (Arcturus)" , 213.918, +19.187}, {"alpha Centauri (Rigil Kentaurus)" , 219.920, -60.835}, {"alpha Scorpii (Antares)" , 247.352, -26.432}, {"alpha Lyrae (Vega)" , 279.234, +38.783}, {"beta Cygni (Albireo)" , 292.680, +27.960}, {"alpha Aquilae (Altair)" , 297.695, +08.867}, {"alpha Cygni (Deneb)" , 310.358, +45.280}, {"alpha Piscis Austrini (Fomalhaut)" , 344.412, -29.622}, {"User Defined Star 1" , 0.000, +00.000}, {"User Defined Star 2" , 0.000, +00.000} }; main () { struct timeval tv; double ra1, dec1; double ra2, dec2; double ra3, dec3; double ra31, dec31; double ra32, dec32; double ra33; double m31, m32, m33; double x1, y1, z1; double x2, y2, z2; double x31, y31, z31; double x32, y32, z32; double X1, Y1, Z1; double X2, Y2, Z2; double X3, Y3, Z3; double d1, d2; double Dxy, Dxz, Dyz, Dzy; double Ddy, Ddz; double my, by; double mz, bz; double a, b, c; double DXY, DYZ, DZX; double P3; double p31, p32; double a12; double A12, A13, A23; double d12; double tx1, ty1, tz1; double tx2, ty2, tz2; double ux1, uy1, uz1; double ux2, uy2, uz2; double skew; int starcount; int i; int h31, h32, h33; int s1, s2; int start_time; int curr_time; /* Read in the stellar database and choose alignment stars. */ starcount = sizeof (table)/sizeof (stardata); for (i = 0; i < starcount; i++) printf ("%2d. %-36s %7.3f %7.3f\n", i+1, table[i].name, table[i].ra, table[i].dec); s1 = s2 = 0; while (s1 == s2) { printf ("\nEnter calibration stars by index number: "); scanf ("%d %d", &s1, &s2); if (s1 == s2) fprintf (stderr, "You must enter two different stars.\n"); } if (s1 > 30) { printf ("Enter the true RA and Dec of %s: ", table[s1-1].name); scanf ("%lf %lf", &table[s1-1].ra, &table[s1-1].dec); } if (s2 > 30) { printf ("Enter the true RA and Dec of %s: ", table[s2-1].name); scanf ("%lf %lf", &table[s2-1].ra, &table[s2-1].dec); } /* Get unit sphere rectangular coordinates for observed and real */ /* positions. */ printf ("\nPoint to %s\n", table[s1-1].name); printf (" Enter RA, Dec: "); scanf ("%lf %lf", &ra1, &dec1); x1 = cos(dec1*DEGREES)*cos(ra1*DEGREES); y1 = cos(dec1*DEGREES)*sin(ra1*DEGREES); z1 = sin(dec1*DEGREES); X1 = cos(table[s1-1].dec*DEGREES)*cos(table[s1-1].ra*DEGREES); Y1 = cos(table[s1-1].dec*DEGREES)*sin(table[s1-1].ra*DEGREES); Z1 = sin(table[s1-1].dec*DEGREES); printf ("\nPoint to %s\n", table[s2-1].name); printf (" Enter RA, Dec: "); scanf ("%lf %lf", &ra2, &dec2); x2 = cos(dec2*DEGREES)*cos(ra2*DEGREES); y2 = cos(dec2*DEGREES)*sin(ra2*DEGREES); z2 = sin(dec2*DEGREES); X2 = cos(table[s2-1].dec*DEGREES)*cos(table[s2-1].ra*DEGREES); Y2 = cos(table[s2-1].dec*DEGREES)*sin(table[s2-1].ra*DEGREES); Z2 = sin(table[s2-1].dec*DEGREES); /* Compute the observed and real separation angles. */ a12 = acos (x1*x2+y1*y2+z1*z2); A12 = acos (X1*X2+Y1*Y2+Z1*Z2); if (a12 <= A12) skew = 0.5*acos (a12/A12); else skew = -0.5*acos (A12/a12); gettimeofday (&tv, 0); start_time = tv.tv_sec; printf ("\nCalibration completed with skew factor = %.6f\n", skew/DEGREES); if (fabs (y1) < 0.000001 || fabs (y2) < 0.000001 || fabs (z1) < 0.000001 || fabs (z2) < 0.000001) { fprintf (stderr, "However, these values are likely to result in "); fprintf (stderr, "substantial errors.\n"); } printf ("\n"); while (1) { /* Get real target position and determine separation from each */ /* alignment stars. */ printf ("Enter desired RA, Dec: "); scanf ("%lf %lf", &ra3, &dec3); X3 = cos (dec3*DEGREES)*cos(ra3*DEGREES); Y3 = cos (dec3*DEGREES)*sin(ra3*DEGREES); Z3 = sin (dec3*DEGREES); A13 = acos (X1*X3+Y1*Y3+Z1*Z3); A23 = acos (X2*X3+Y2*Y3+Z2*Z3); /* We now advance the time and recompute coordinates. */ /* To account for the scope axis being off from the true */ /* celestial axis, we note that the target is where another */ /* star would have been x hours ago, but shifted down by x */ /* hours of right ascension. We then figure out where that */ /* star would have been with respect to the alignment stars at */ /* that time, then shift back up by x hours at the end. */ gettimeofday (&tv, 0); curr_time = tv.tv_sec; ra3 -= (double) (curr_time-start_time)/240.0; X3 = cos (dec3*DEGREES)*cos(ra3*DEGREES); Y3 = cos (dec3*DEGREES)*sin(ra3*DEGREES); Z3 = sin (dec3*DEGREES); d1 = X1*X3+Y1*Y3+Z1*Z3; d2 = X2*X3+Y2*Y3+Z2*Z3; d12 = sqrt ((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)+(z2-z1)*(z2-z1)); /* We also need to correct for errors in alignment. */ /* If the observed and true separation angles (a12 and A12) are */ /* off by even a small amount, this can magnify errors across */ /* the celestial globe, particularly along the great circle */ /* connecting the two alignment stars. This portion */ /* recalibrates the two alignment stars each time a new */ /* position is requested, dependent on how far the target is */ /* from each alignment star. */ tx1 = x1+2.0*A13/(A13+A23)*cos(a12/2.0)*(tan(a12/2.0)-tan(A12/2.0))* (x2-x1)/d12; tx2 = x2+2.0*A23/(A13+A23)*cos(a12/2.0)*(tan(a12/2.0)-tan(A12/2.0))* (x1-x2)/d12; ty1 = y1+2.0*A13/(A13+A23)*cos(a12/2.0)*(tan(a12/2.0)-tan(A12/2.0))* (y2-y1)/d12; ty2 = y2+2.0*A23/(A13+A23)*cos(a12/2.0)*(tan(a12/2.0)-tan(A12/2.0))* (y1-y2)/d12; tz1 = z1+2.0*A13/(A13+A23)*cos(a12/2.0)*(tan(a12/2.0)-tan(A12/2.0))* (z2-z1)/d12; tz2 = z2+2.0*A23/(A13+A23)*cos(a12/2.0)*(tan(a12/2.0)-tan(A12/2.0))* (z1-z2)/d12; ux1 = tx1/sqrt (tx1*tx1+ty1*ty1+tz1*tz1); uy1 = ty1/sqrt (tx1*tx1+ty1*ty1+tz1*tz1); uz1 = tz1/sqrt (tx1*tx1+ty1*ty1+tz1*tz1); ux2 = tx2/sqrt (tx2*tx2+ty2*ty2+tz2*tz2); uy2 = ty2/sqrt (tx2*tx2+ty2*ty2+tz2*tz2); uz2 = tz2/sqrt (tx2*tx2+ty2*ty2+tz2*tz2); /* The rest is matrix algebra to solve simultaneous equations. */ Dxy = ux1*uy2-ux2*uy1; Dxz = ux1*uz2-ux2*uz1; Dyz = uy1*uz2-uy2*uz1; Dzy = uz1*uy2-uz2*uy1; Ddy = d1*uy2-d2*uy1; Ddz = d1*uz2-d2*uz1; DXY = X1*Y2-X2*Y1; DYZ = Y1*Z2-Y2*Z1; DZX = Z1*X2-Z2*X1; P3 = DYZ*X3+DZX*Y3+DXY*Z3; my = Dxz/Dyz; by = Ddz/Dyz; mz = Dxy/Dzy; bz = Ddy/Dzy; /* We end up with a quadratic equation... */ a = my*my+mz*mz+1.0; b = -2.0*(my*by+mz*bz); c = by*by+bz*bz-1.0; /* ...which of course has two different solutions. */ /* Both solutions are equally distant from each alignment star. */ /* How do we distinguish between them? Let v1 and v2 be the */ /* position vectors for the two alignment stars, and v3 be the */ /* position vector for the target. Then we compute */ /* (v1 x v2) . v3, where x and . are cross and dot product, */ /* respectively. The result will be positive for one of the */ /* solutions and negative for the other--we just choose the */ /* one that most closely matches the result for the real */ /* positions. */ if (b*b > 4.0*a*c) x31 = (-b+sqrt (b*b-4.0*a*c))/(2.0*a); else x31 = (-b)/(2.0*a); y31 = by-my*x31; z31 = bz-mz*x31; p31 = Dyz*x31-Dxz*y31+Dxy*z31; dec31 = asin(z31)/DEGREES; ra31 = atan(y31/x31)/DEGREES; if (x31 < 0.0) ra31 += 180.0; if (ra31 < 0.0) ra31 += 360.0; ra31 += (double) (curr_time-start_time)/240.0; if (ra31 > 360.0) ra31 -= 360.0; h31 = (int) (ra31/15.0); m31 = (ra31-h31*15.0)*4.0; if (b*b > 4.0*a*c) x32 = (-b-sqrt (b*b-4.0*a*c))/(2.0*a); else x32 = (-b)/(2.0*a); y32 = by-my*x32; z32 = bz-mz*x32; p32 = Dyz*x32-Dxz*y32+Dxy*z32; dec32 = asin(z32)/DEGREES; ra32 = atan(y32/x32)/DEGREES; if (x32 < 0.0) ra32 += 180.0; if (ra32 < 0.0) ra32 += 360.0; ra32 += (double) (curr_time-start_time)/240.0; if (ra32 > 360.0) ra32 -= 360.0; h32 = (int) (ra32/15.0); m32 = (ra32-h32*15.0)*4.0; if (fabs (p32-P3) > fabs (p31-P3)) printf ("Point to RA %7.3f (%2dh %4.1fm), Dec %7.3f\n", ra31, h31, m31, dec31); else printf ("Point to RA %7.3f (%2dh %4.1fm), Dec %7.3f\n", ra32, h32, m32, dec32); } }