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
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
|
define J2000 2000.0D0 # J2000
define JD2000 2451545.0D0 # J2000 Julian Date
define JYEAR 365.25D0 # Julian year
# AST_DATE_TO_EPOCH -- Convert Gregorian date and solar mean time to
# a Julian epoch. A Julian epoch has 365.25 days per year and 24
# hours per day.
procedure ast_date_to_epoch (year, month, day, ut, epoch)
int year # Year
int month # Month (1-12)
int day # Day of month
double ut # Universal time for date (mean solar day)
double epoch # Julian epoch
double jd, ast_date_to_julday()
begin
jd = ast_date_to_julday (year, month, day, ut)
epoch = J2000 + (jd - JD2000) / JYEAR
end
# AST_EPOCH_TO_DATE -- Convert a Julian epoch to year, month, day, and time.
procedure ast_epoch_to_date (epoch, year, month, day, ut)
double epoch # Julian epoch
int year # Year
int month # Month (1-12)
int day # Day of month
double ut # Universal time for date
double jd
begin
jd = JD2000 + (epoch - J2000) * JYEAR
call ast_julday_to_date (jd, year, month, day, ut)
end
# AST_DAY_OF_YEAR -- The day number for the given year is returned.
int procedure ast_day_of_year (year, month, day)
int year # Year
int month # Month (1-12)
int day # Day of month
int d
int bom[13] # Beginning of month
data bom/1,32,60,91,121,152,182,213,244,274,305,335,366/
begin
d = bom[month] + day - 1
if (month > 2 && mod (year, 4) == 0 &&
(mod (year, 100) != 0 || mod (year, 400) == 0))
d = d + 1
return (d)
end
# AST_DAY_OF_WEEK -- Return the day of the week for the given Julian day.
# The integer day of the week is 0=Sunday - 6=Saturday. The character string
# is the three character abbreviation for the day of the week. Note that
# the day of the week is for Greenwich if the standard UT is used.
procedure ast_day_of_week (jd, d, name, sz_name)
double jd # Julian date
int d # Day of the week (0=SUN)
char name[sz_name] # Name for day of the week
int sz_name # Size of name string
begin
d = mod (int (jd - 0.5) + 2, 7)
switch (d) {
case 0:
call strcpy ("SUN", name, sz_name)
case 1:
call strcpy ("MON", name, sz_name)
case 2:
call strcpy ("TUE", name, sz_name)
case 3:
call strcpy ("WED", name, sz_name)
case 4:
call strcpy ("THU", name, sz_name)
case 5:
call strcpy ("FRI", name, sz_name)
case 6:
call strcpy ("SAT", name, sz_name)
}
end
# AST_JULDAY -- Convert epoch to Julian day.
double procedure ast_julday (epoch)
double epoch # Epoch
double jd
begin
jd = JD2000 + (epoch - J2000) * JYEAR
return (jd)
end
# AST_DATE_TO_JULDAY -- Convert date to Julian day.
# This assumes dates after year 99.
double procedure ast_date_to_julday (year, month, day, t)
int year # Year
int month # Month (1-12)
int day # Day of month
double t # Time for date (mean solar day)
double jd
int y, m, d
begin
if (year < 100)
y = 1900 + year
else
y = year
if (month > 2)
m = month + 1
else {
m = month + 13
y = y - 1
}
jd = int (JYEAR * y) + int (30.6001 * m) + day + 1720995
if (day + 31 * (m + 12 * y) >= 588829) {
d = int (y / 100)
m = int (y / 400)
jd = jd + 2 - d + m
}
jd = jd - 0.5 + int (t * 360000. + 0.5) / 360000. / 24.
return (jd)
end
# AST_JULDAY_TO_DATE -- Convert Julian date to calendar date.
# This is taken from Numerical Receipes by Press, Flannery, Teukolsy, and
# Vetterling.
procedure ast_julday_to_date (j, year, month, day, t)
double j # Julian day
int year # Year
int month # Month (1-12)
int day # Day of month
double t # Time for date (mean solar day)
int ja, jb, jc, jd, je
begin
ja = nint (j)
t = 24. * (j - ja + 0.5)
if (ja >= 2299161) {
jb = int (((ja - 1867216) - 0.25) / 36524.25)
ja = ja + 1 + jb - int (jb / 4)
}
jb = ja + 1524
jc = int (6680. + ((jb - 2439870) - 122.1) / JYEAR)
jd = 365 * jc + int (jc / 4)
je = int ((jb - jd) / 30.6001)
day = jb - jd - int (30.6001 * je)
month = je - 1
if (month > 12)
month = month - 12
year = jc - 4715
if (month > 2)
year = year - 1
if (year < 0)
year = year - 1
end
# AST_MST -- Mean sidereal time of the epoch at the given longitude.
# This procedure may be used to optain Greenwich Mean Sidereal Time (GMST)
# by setting the longitude to 0.
double procedure ast_mst (epoch, longitude)
double epoch # Epoch
double longitude # Longitude in degrees
double jd, ut, t, st
double ast_julday()
begin
# Determine JD and UT, and T (JD in centuries from J2000.0).
jd = ast_julday (epoch)
ut = (jd - int (jd) - 0.5) * 24.
t = (jd - 2451545.d0) / 36525.d0
# The GMST at 0 UT in seconds is a power series in T.
st = 24110.54841d0 +
t * (8640184.812866d0 + t * (0.093104d0 - t * 6.2d-6))
# Correct for longitude and convert to standard hours.
st = mod (st / 3600. + ut - longitude / 15., 24.0D0)
if (st < 0)
st = st + 24
return (st)
end
|
