/* * photom: a package to perform photometric calibration * Copyright (C) 2001 Michael William Richmond * * Contact: Michael William Richmond * Physics Department * Rochester Institute of Technology * 85 Lomb Memorial Drive * Rochester, NY 14623-5603 * E-mail: mwrsps@rit.edu * * * This program is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License * as published by the Free Software Foundation; either version 2 * of the License, or (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * */ /* precess coordinates from some other epoch to J2000.0, or vice versa. the coordinates should be expressed in decimal degrees, i.e., RA 0 -> 360 Dec -90 -> +90 */ #include #include #include "precess.h" #ifndef M_PI #define M_PI 3.141592654 #endif /* multiply by the following factor to convert degrees to radians */ #define DEGTORAD (M_PI/180.0) /* #define MAIN */ /* define this if you want the 'main()' in this file */ /* the precession constants for conversions to/from J2000.0, from the Astronomical Almanac for 1989. */ double zeta1 = 0.6406161; double zeta2 = 0.0000839; double zeta3 = 0.0000050; double z1 = 0.6406161; double z2 = 0.0003041; double z3 = 0.0000051; double theta1 = 0.5567530; double theta2 = -0.0001185; double theta3 = -0.0000116; /* a silly little test 'main()' to see that I got the conversions correct */ #ifdef MAIN char progname[] = "precess"; main() { double oldra, olddec, epoch, newra, newdec; double hh, deg, seconds; int hours, minutes; char line[100]; printf("gimme an epoch other than J2000.0 "); scanf("%lf", &epoch); printf("gimme the RA and Dec is to be converted: RA (HH MM SS.sss) "); scanf("%d %d %lf", &hours, &minutes, &seconds); fflush(stdout); ratodeg(hours, minutes, seconds, &oldra); printf(" Dec (DD MM SS.sss) "); scanf("%d %d %lf", &hours, &minutes, &seconds); fflush(stdout); hmstohhh(hours, minutes, seconds, &olddec); to_j2000(epoch, oldra, olddec, &newra, &newdec); degtora(newra, &hours, &minutes, &seconds); format_hms(hours, minutes, seconds, line); printf("%6lf -> J2000.0, RA %s ", epoch, line); hhhtohms(newdec, &hours, &minutes, &seconds); format_hms(hours, minutes, seconds, line); printf(" Dec %s\n", line); from_j2000(epoch, oldra, olddec, &newra, &newdec); degtora(newra, &hours, &minutes, &seconds); format_hms(hours, minutes, seconds, line); printf("J2000.0 -> %6lf, RA %s ", epoch, line); hhhtohms(newdec, &hours, &minutes, &seconds); format_hms(hours, minutes, seconds, line); printf(" Dec %s\n", line); } #endif /* MAIN */ /* convert the coordinates, 'oldra' and 'olddec', for epoch J2000.0, to their values in epoch 'epoch'; the new values are placed in 'newra' and 'newdec'. */ void from_j2000(epoch, oldra, olddec, newra, newdec) double epoch, oldra, olddec, *newra, *newdec; { double T, zeta, z, theta, angle1, angle2; T = (epoch - 2000.0)/100.0; zeta = (zeta1*T + zeta2*T*T + zeta3*T*T*T)*DEGTORAD; z = (z1*T + z2*T*T + z3*T*T*T)*DEGTORAD; theta = (theta1*T + theta2*T*T + theta3*T*T*T)*DEGTORAD; /* convert RA and Dec from 360-degree to radian measure */ oldra *= DEGTORAD; olddec *= DEGTORAD; *newdec = asin ( cos(oldra + zeta)*sin(theta)*cos(olddec) + cos(theta)*sin(olddec) ); angle1 = asin ( sin(oldra + zeta)*cos(olddec)/cos(*newdec) ); angle2 = acos ( (cos(oldra + zeta)*cos(theta)*cos(olddec) - sin(theta)*sin(olddec))/cos(*newdec) ); /* put things back into degrees for quadrant calculations */ *newdec /= DEGTORAD; z /= DEGTORAD; angle1 /= DEGTORAD; angle2 /= DEGTORAD; /* use both angle1 and angle2 to decide which quadrant the new RA is in */ if ((angle1 > 0) && (angle1 < 90.0) && (angle2 > 0) && (angle2 < 90.0)) *newra = angle1 + z; /* first quadrant */ else if ((angle1 > 0) && (angle1 < 90) && (angle2 > 90) && (angle2 < 180)) *newra = angle2 + z; /* second quadrant */ else if ((angle1 < 0) && (angle1 > -90) && (angle2 > 90) && (angle2 < 180)) *newra = (180.0 - angle1) + z; /* third quadrant */ else *newra = (360.0 + angle1) + z; /* fourth quadrant */ if (*newra < 0) *newra += 360.0; } /* convert the coordinates, 'oldra' and 'olddec', for the given epoch, 'epoch', to their values in epoch J2000.0; the new values are placed in 'newra' and 'newdec'. */ void to_j2000(epoch, oldra, olddec, newra, newdec) double epoch, oldra, olddec, *newra, *newdec; { double T, zeta, z, theta, angle1, angle2; if (epoch == 2000.0) { *newra = oldra; *newdec = olddec; return; } T = (epoch - 2000.0)/100.0; zeta = (zeta1*T + zeta2*T*T + zeta3*T*T*T)*DEGTORAD; z = (z1*T + z2*T*T + z3*T*T*T)*DEGTORAD; theta = (theta1*T + theta2*T*T + theta3*T*T*T)*DEGTORAD; /* convert RA and Dec from 360-degree to radian measure */ oldra *= DEGTORAD; olddec *= DEGTORAD; *newdec = asin ( -cos(oldra - z)*sin(theta)*cos(olddec) + cos(theta)*sin(olddec) ); angle1 = asin ( sin(oldra - z)*cos(olddec)/cos(*newdec) ); angle2 = acos ( (cos(oldra - z)*cos(theta)*cos(olddec) + sin(theta)*sin(olddec))/cos(*newdec) ); /* put things back into degrees for quadrant calculations */ *newdec /= DEGTORAD; zeta /= DEGTORAD; angle1 /= DEGTORAD; angle2 /= DEGTORAD; /* use both angle1 and angle2 to decide which quadrant the new RA is in */ if ((angle1 > 0) && (angle1 < 90.0) && (angle2 > 0) && (angle2 < 90.0)) *newra = angle1 - zeta; /* first quadrant */ else if ((angle1 > 0) && (angle1 < 90) && (angle2 > 90) && (angle2 < 180)) *newra = angle2 - zeta; /* second quadrant */ else if ((angle1 < 0) && (angle1 > -90) && (angle2 > 90) && (angle2 < 180)) *newra = (180.0 - angle1) - zeta; /* third quadrant */ else *newra = (360.0 + angle1) - zeta; /* fourth quadrant */ if (*newra < 0) *newra += 360.0; if (*newra == 360.0) *newra = 0.0; }