1 /* code to convert between equitorial and galactic coordinates */
2
3 #include <stdio.h>
4 #include <math.h>
5
6 #include "astro.h"
7
8 static void galeq_aux (int sw, double x, double y, double *p, double *q);
9 static void galeq_init (void);
10
11 #define EQtoGAL 1
12 #define GALtoEQ (-1)
13 #define SMALL (1e-20)
14
15 static double an = degrad(32.93192); /* G lng of asc node on equator */
16 static double gpr = degrad(192.85948); /* RA of North Gal Pole, 2000 */
17 static double gpd = degrad(27.12825); /* Dec of " */
18 static double cgpd, sgpd; /* cos() and sin() of gpd */
19 static double mj2000; /* mj of 2000 */
20 static int before; /* whether these have been set yet */
21
22 /* given ra and dec, each in radians, for the given epoch, find the
23 * corresponding galactic latitude, *lt, and longititude, *lg, also each in
24 * radians.
25 */
26 void
eq_gal(double mj,double ra,double dec,double * lt,double * lg)27 eq_gal (double mj, double ra, double dec, double *lt, double *lg)
28 {
29 galeq_init();
30 precess (mj, mj2000, &ra, &dec);
31 galeq_aux (EQtoGAL, ra, dec, lg, lt);
32 }
33
34 /* given galactic latitude, lt, and longititude, lg, each in radians, find
35 * the corresponding equitorial ra and dec, also each in radians, at the
36 * given epoch.
37 */
38 void
gal_eq(double mj,double lt,double lg,double * ra,double * dec)39 gal_eq (double mj, double lt, double lg, double *ra, double *dec)
40 {
41 galeq_init();
42 galeq_aux (GALtoEQ, lg, lt, ra, dec);
43 precess (mj2000, mj, ra, dec);
44 }
45
46 static void
galeq_aux(int sw,double x,double y,double * p,double * q)47 galeq_aux (
48 int sw, /* +1 for eq to gal, -1 for vv. */
49 double x, double y, /* sw==1: x==ra, y==dec. sw==-1: x==lg, y==lt. */
50 double *p, double *q) /* sw==1: p==lg, q==lt. sw==-1: p==ra, q==dec. */
51 {
52 double sy, cy, a, ca, sa, b, sq, c, d;
53
54 cy = cos(y);
55 sy = sin(y);
56 a = x - an;
57 if (sw == EQtoGAL)
58 a = x - gpr;
59 ca = cos(a);
60 sa = sin(a);
61 b = sa;
62 if (sw == EQtoGAL)
63 b = ca;
64 sq = (cy*cgpd*b) + (sy*sgpd);
65 *q = asin (sq);
66
67 if (sw == GALtoEQ) {
68 c = cy*ca;
69 d = (sy*cgpd) - (cy*sgpd*sa);
70 if (fabs(d) < SMALL)
71 d = SMALL;
72 *p = atan (c/d) + gpr;
73 } else {
74 c = sy - (sq*sgpd);
75 d = cy*sa*cgpd;
76 if (fabs(d) < SMALL)
77 d = SMALL;
78 *p = atan (c/d) + an;
79 }
80
81 if (d < 0) *p += PI;
82 if (*p < 0) *p += 2*PI;
83 if (*p > 2*PI) *p -= 2*PI;
84 }
85
86 /* set up the definitions */
87 static void
galeq_init()88 galeq_init()
89 {
90 if (!before) {
91 cgpd = cos (gpd);
92 sgpd = sin (gpd);
93 mj2000 = J2000;
94 before = 1;
95 }
96 }
97
98