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
146
147
148
149
150
151
152
153
154
|
include <math.h>
# PRECESSMGB -- Calculate the apparent change in right ascension and
# declination between two epochs. Corrections are made for precession,
# aberration, and some terms of nutation.
#
# Based on the work of R. N. Manchester, M. A. Gordon, J. A. Ball (NRAO)
# January 1970 (DOPSET, IBM360)
# Modified by E. Anderson (NOAO) April 1985 (DOPSET, VMS/VAX)
# Converted to IRAF by F. Valdes (NOAO) March 1986.
procedure precessmgb (ra1, dec1, epoch1, ra2, dec2, epoch2)
double ra1, dec1, epoch1 # First epoch coordinates
double ra2, dec2, epoch2 # Second epoch coordinates
double ra, dec, epoch, epoch3, dc
bool fp_equald()
begin
# Return if the epochs are the same.
if (fp_equald (epoch1, epoch2)) {
ra2 = ra1
dec2 = dec1
return
}
ra = ra1
dec = dec1
epoch = epoch1
# Check if epoch1 is the beginning of a year. If not make correction
# to the beginning of the year.
if (epoch - int (epoch) > 1E-6) {
epoch3 = epoch
epoch = int (epoch)
call mgb_precess (ra, dec, epoch, ra2, dec2, epoch3, dc)
ra = ra - (ra2 - ra)
dec = dec - (dec2 - dec)
}
# Precess to epoch2.
call mgb_precess (ra, dec, epoch, ra2, dec2, epoch2, dc)
end
# MGB_PRECESS -- Calculate the apparent change in right ascension and
# declination between two epochs. The first epoch is assumed to be the
# beginning of a year (eg 1950.0). Also calculate the equation of the
# equinoxes (in minutes of time) which may be added to the mean siderial
# time to give the apparent siderial time (AENA-469). Corrections are
# made for precession, aberration, and some terms of nutation.
# AENA - The American Ephemeris and Nautical Almanac (the blue book)
# ESE - The explanatory suplement to AENA (the green book)
procedure mgb_precess (ra1, dec1, epoch1, ra2, dec2, epoch2, dc)
double ra1, dec1, epoch1 # First epoch coordinates
double ra2, dec2, epoch2 # Second epoch coordinates
double dc # Siderial time correction
double ra, dec, delr, deld
double t1, t2, nday2
double theta, zeta, z, am, an, al
double csr, snr, csd, snd, tnd, csl, snl
double omega, dlong, doblq
bool fp_equald()
begin
# Check if epochs are the same.
if (fp_equald (epoch1, epoch2)) {
ra2 = ra1
dec2 = dec1
return
}
# Convert input coordinates to radians.
ra = DEGTORAD (ra1 * 15.)
dec = DEGTORAD (dec1)
# Compute sines, cosines, and tangents.
csr = cos (ra)
snr = sin (ra)
csd = cos (dec)
snd = sin (dec)
tnd = snd / csd
# T1 is the time from 1900 to epoch1 (centuries),
# t2 is the time from epoch1 to epoch2 (centuries), and
# nday2 is the number of days since the beginning of the year
# for epoch2. The number of ephemeris days in a tropical
# year is 365.2421988.
t1 = (epoch1 - 1900.) / 100.
t2 = (epoch2 - epoch1) / 100.
nday2 = (epoch2 - int (epoch2)) * 365.2421988
# Theta, zeta, and z are precessional angles from ESE-29 (arcseconds).
theta = t2 * ((2004.682 - 0.853 * t1) - t2 * (0.426 + t2 * 0.042))
zeta = t2 * ((2304.250 + 1.396 * t1) + t2 * (0.302 + t2 * 0.018))
z = zeta + 0.791 * t2 ** 2
# am and an are the M and N precessional numbers (see AENA-50, 474)
# (radians) and alam is an approximate mean longitude for the sun
# (AENA-50) (radians)
am = DEGTORAD ((zeta + z) / 3600.)
an = DEGTORAD (theta / 3600.)
al = DEGTORAD (0.985647 * nday2 + 278.5)
snl = sin (al)
csl = cos (al)
# Delr and deld are the annual aberation term in ra and dec (radians)
# (ESE-47,48) (0.91745051 cos (obliquity of ecliptic))
# (0.39784993 sin (obliquity of ecliptic))
# (-9.92413605E-5 K 20.47 ARCSECONDS constant of aberration)
# (ESE) plus precession terms (see AENA-50 and ESE-39).
delr = -9.92413605e-5 * (snl * snr + 0.91745051 * csl * csr) / csd +
am + an * snr * tnd
deld = -9.92413605e-5 * (snl * csr * snd - 0.91745051 * csl * snr *
snd + 0.39784993 * csl * csd) + an * csr
# The following calculates the nutation (approximately) (ESE-41,45)
# Omega is the angle of the first term of nutation (ESE-44)
# (approximate formula) (radians).
# Dlong is the nutation in longitude (delta-psi) (radians)
# Doblq is the nutation in obliquity (delta-epsilon) (radians)
omega = DEGTORAD (259.183275 - 1934.142 * (t1 + t2))
dlong = -8.3597e-5 * sin (omega)
doblq = 4.4678e-5 * cos (omega)
# Add nutation into delr and deld (ESE-43).
delr = delr + dlong * (0.91745051 + 0.39784993 * snr * tnd) -
csr * tnd * doblq
deld = deld + 0.39784993 * csr * dlong + snr * doblq
# Compute new position and the equation of the equinoxes
# (dc in minutes of tim, ESE-43)
ra2 = ra1 + RADTODEG (delr / 15.)
dec2 = dec1 + RADTODEG (deld)
dc = dlong * 210.264169
end
|
