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