1 #include <stdio.h>
2 #include <math.h>
3
4 #include "astro.h"
5
6 static void precess_hiprec (double mjd1, double mjd2, double *ra, double *dec);
7
8
9 #define DCOS(x) cos(degrad(x))
10 #define DSIN(x) sin(degrad(x))
11 #define DASIN(x) raddeg(asin(x))
12 #define DATAN2(y,x) raddeg(atan2((y),(x)))
13
14 /* corrects ra and dec, both in radians, for precession from epoch 1 to epoch 2.
15 * the epochs are given by their modified JDs, mjd1 and mjd2, respectively.
16 * N.B. ra and dec are modifed IN PLACE.
17 */
18 void
precess(double mjd1,double mjd2,double * ra,double * dec)19 precess (
20 double mjd1, double mjd2, /* initial and final epoch modified JDs */
21 double *ra, double *dec) /* ra/dec for mjd1 in, for mjd2 out */
22 {
23 precess_hiprec (mjd1, mjd2, ra, dec);
24 }
25
26 /*
27 * Copyright (c) 1990 by Craig Counterman. All rights reserved.
28 *
29 * This copy of the precess_hiprec() routine is, by permission of the
30 * author, licensed under the same license as the PyEphem package in
31 * which it is here included.
32 *
33 * Rigorous precession. From Astronomical Ephemeris 1989, p. B18
34 *
35 * 96-06-20 Hayo Hase <hase@wettzell.ifag.de>: theta_a corrected
36 */
37 static void
precess_hiprec(double mjd1,double mjd2,double * ra,double * dec)38 precess_hiprec (
39 double mjd1, double mjd2, /* initial and final epoch modified JDs */
40 double *ra, double *dec) /* ra/dec for mjd1 in, for mjd2 out */
41 {
42 static double last_mjd1 = -213.432, last_from;
43 static double last_mjd2 = -213.432, last_to;
44 double zeta_A, z_A, theta_A;
45 double T;
46 double A, B, C;
47 double alpha, delta;
48 double alpha_in, delta_in;
49 double from_equinox, to_equinox;
50 double alpha2000, delta2000;
51
52 /* convert mjds to years;
53 * avoid the remarkably expensive calls to mjd_year()
54 */
55 if (last_mjd1 == mjd1)
56 from_equinox = last_from;
57 else {
58 mjd_year (mjd1, &from_equinox);
59 last_mjd1 = mjd1;
60 last_from = from_equinox;
61 }
62 if (last_mjd2 == mjd2)
63 to_equinox = last_to;
64 else {
65 mjd_year (mjd2, &to_equinox);
66 last_mjd2 = mjd2;
67 last_to = to_equinox;
68 }
69
70 /* convert coords in rads to degs */
71 alpha_in = raddeg(*ra);
72 delta_in = raddeg(*dec);
73
74 /* precession progresses about 1 arc second in .047 years */
75 /* From from_equinox to 2000.0 */
76 if (fabs (from_equinox-2000.0) > .02) {
77 T = (from_equinox - 2000.0)/100.0;
78 zeta_A = 0.6406161* T + 0.0000839* T*T + 0.0000050* T*T*T;
79 z_A = 0.6406161* T + 0.0003041* T*T + 0.0000051* T*T*T;
80 theta_A = 0.5567530* T - 0.0001185* T*T - 0.0000116* T*T*T;
81
82 A = DSIN(alpha_in - z_A) * DCOS(delta_in);
83 B = DCOS(alpha_in - z_A) * DCOS(theta_A) * DCOS(delta_in)
84 + DSIN(theta_A) * DSIN(delta_in);
85 C = -DCOS(alpha_in - z_A) * DSIN(theta_A) * DCOS(delta_in)
86 + DCOS(theta_A) * DSIN(delta_in);
87
88 alpha2000 = DATAN2(A,B) - zeta_A;
89 range (&alpha2000, 360.0);
90 delta2000 = DASIN(C);
91 } else {
92 /* should get the same answer, but this could improve accruacy */
93 alpha2000 = alpha_in;
94 delta2000 = delta_in;
95 };
96
97
98 /* From 2000.0 to to_equinox */
99 if (fabs (to_equinox - 2000.0) > .02) {
100 T = (to_equinox - 2000.0)/100.0;
101 zeta_A = 0.6406161* T + 0.0000839* T*T + 0.0000050* T*T*T;
102 z_A = 0.6406161* T + 0.0003041* T*T + 0.0000051* T*T*T;
103 theta_A = 0.5567530* T - 0.0001185* T*T - 0.0000116* T*T*T;
104
105 A = DSIN(alpha2000 + zeta_A) * DCOS(delta2000);
106 B = DCOS(alpha2000 + zeta_A) * DCOS(theta_A) * DCOS(delta2000)
107 - DSIN(theta_A) * DSIN(delta2000);
108 C = DCOS(alpha2000 + zeta_A) * DSIN(theta_A) * DCOS(delta2000)
109 + DCOS(theta_A) * DSIN(delta2000);
110
111 alpha = DATAN2(A,B) + z_A;
112 range(&alpha, 360.0);
113 delta = DASIN(C);
114 } else {
115 /* should get the same answer, but this could improve accruacy */
116 alpha = alpha2000;
117 delta = delta2000;
118 };
119
120 *ra = degrad(alpha);
121 *dec = degrad(delta);
122 }
123
124 #if 0
125 static void
126 precess_fast (
127 double mjd1, double mjd2, /* initial and final epoch modified JDs */
128 double *ra, double *dec) /* ra/dec for mjd1 in, for mjd2 out */
129 {
130 #define N degrad (20.0468/3600.0)
131 #define M hrrad (3.07234/3600.0)
132 double nyrs;
133
134 nyrs = (mjd2 - mjd1)/365.2425;
135 *dec += N * cos(*ra) * nyrs;
136 *ra += (M + (N * sin(*ra) * tan(*dec))) * nyrs;
137 range (ra, 2.0*PI);
138 }
139 #endif
140
141