« back to all changes in this revision
Viewing changes to CalcEphem.c
- browse files at revision 1
- compare with another revision
- download diff
- download tarball
- view history from revision 1
- Committer: Bazaar Package Importer
- Author(s): Martin Zobel-Helas
- Date: 2004-11-19 22:55:21 UTC
- Revision ID: james.westby@ubuntu.com-20041119225521-pmuy6kl327g5wevc
Tags: upstream-0.6
Import upstream version 0.6
added
removed
CalcEphem.c
1
/* (C) Mike Henderson <mghenderson@lanl.gov>.
2
*
3
* I've converted to glib types, removed unused variables and
4
* piped the whole thing through indent.
5
*
6
* Josh Buhl <jbuhl@users.sourceforge.net>
7
*/
8
#ifdef HAVE_CONFIG_H
9
#include <config.h>
10
#endif
11
12
#include "CalcEphem.h"
13
#include "Moon.h"
14
15
16
static gdouble kepler(gdouble M, gdouble e)
17
{
18
gint n = 0;
19
gdouble E, Eold, eps = 1.0e-8;
20
21
E = M + e * sin(M);
22
do {
23
Eold = E;
24
E = Eold + (M - Eold + e * sin(Eold))
25
/ (1.0 - e * cos(Eold));
26
++n;
27
} while ((fabs(E - Eold) > eps) && (n < 100));
28
return (E);
29
}
30
31
static gint DayofYear(gint year, gint month, gint day)
32
{
33
gdouble jd();
34
return ((gint) (jd(year, month, day, 0.0) - jd(year, 1, 0, 0.0)));
35
}
36
37
38
static gint DayofWeek(gint year, gint month, gint day, gchar dowstr[])
39
{
40
gdouble JD, A, Afrac, jd();
41
gint n, iA;
42
43
JD = jd(year, month, day, 0.0);
44
A = (JD + 1.5) / 7.0;
45
iA = (gint) A;
46
Afrac = A - (gdouble) iA;
47
n = (gint) (Afrac * 7.0 + 0.5);
48
switch (n) {
49
case 0:
50
strcpy(dowstr, "Sunday");
51
break;
52
case 1:
53
strcpy(dowstr, "Monday");
54
break;
55
case 2:
56
strcpy(dowstr, "Tuesday");
57
break;
58
case 3:
59
strcpy(dowstr, "Wednesday");
60
break;
61
case 4:
62
strcpy(dowstr, "Thursday");
63
break;
64
case 5:
65
strcpy(dowstr, "Friday");
66
break;
67
case 6:
68
strcpy(dowstr, "Saturday");
69
break;
70
}
71
return (n);
72
}
73
74
/*
75
* Compute the Julian Day number for the given date.
76
* Julian Date is the number of days since noon of Jan 1 4713 B.C.
77
*/
78
gdouble jd(gint ny, gint nm, gint nd, gdouble UT)
79
{
80
gdouble A, B, C, D, JD, day;
81
82
day = nd + UT / 24.0;
83
84
85
if ((nm == 1) || (nm == 2)) {
86
ny = ny - 1;
87
nm = nm + 12;
88
}
89
90
if (((gdouble) ny + nm / 12.0 + day / 365.25) >=
91
(1582.0 + 10.0 / 12.0 + 15.0 / 365.25)) {
92
A = ((gint) (ny / 100.0));
93
B = 2.0 - A + (gint) (A / 4.0);
94
} else {
95
B = 0.0;
96
}
97
98
if (ny < 0.0) {
99
C = (gint) ((365.25 * (gdouble) ny) - 0.75);
100
} else {
101
C = (gint) (365.25 * (gdouble) ny);
102
}
103
104
D = (gint) (30.6001 * (gdouble) (nm + 1));
105
106
107
JD = B + C + D + day + 1720994.5;
108
return (JD);
109
}
110
111
gdouble hour24(gdouble hour)
112
{
113
gint n;
114
115
if (hour < 0.0) {
116
n = (gint) (hour / 24.0) - 1;
117
return (hour - n * 24.0);
118
} else if (hour > 24.0) {
119
n = (gint) (hour / 24.0);
120
return (hour - n * 24.0);
121
} else {
122
return (hour);
123
}
124
}
125
126
gdouble angle2pi(gdouble angle)
127
{
128
gint n;
129
gdouble a;
130
a = 2.0 * M_PI;
131
132
if (angle < 0.0) {
133
n = (gint) (angle / a) - 1;
134
return (angle - n * a);
135
} else if (angle > a) {
136
n = (gint) (angle / a);
137
return (angle - n * a);
138
} else {
139
return (angle);
140
}
141
}
142
143
static gdouble angle360(gdouble angle)
144
{
145
gint n;
146
147
if (angle < 0.0) {
148
n = (gint) (angle / 360.0) - 1;
149
return (angle - n * 360.0);
150
} else if (angle > 360.0) {
151
n = (gint) (angle / 360.0);
152
return (angle - n * 360.0);
153
} else {
154
return (angle);
155
}
156
}
157
158
#if 0
159
static void Radec_to_Cart(gdouble ra, gdouble dec, Vector * r)
160
{
161
162
/*
163
* Convert ra/dec from degrees to radians
164
*/
165
ra *= RadPerDeg;
166
dec *= RadPerDeg;
167
168
169
/*
170
* Compute cartesian coordinates (in GEI)
171
*/
172
r->x = cos(dec) * cos(ra);
173
r->y = cos(dec) * sin(ra);
174
r->z = sin(dec);
175
176
return;
177
}
178
#endif
179
180
#if 0
181
static gint LeapYear(gint year)
182
{
183
if ((year % 100 == 0) && (year % 400 != 0))
184
return (0);
185
else if (year % 4 == 0)
186
return (1);
187
else
188
return (0);
189
}
190
#endif
191
192
void CalcEphem(glong date, gdouble UT, CTrans * c)
193
{
194
gint year, month, day;
195
gdouble TU, TU2, TU3, T0, gmst;
196
gdouble varep, varpi;
197
gdouble eccen, epsilon;
198
gdouble days, M, E, nu, lambnew;
199
gdouble r0, earth_sun_distance;
200
gdouble RA, DEC, RA_Moon, DEC_Moon;
201
gdouble TDT;
202
gdouble AGE;
203
gdouble LambdaMoon, BetaMoon, R;
204
gdouble Ta, Tb, Tc, frac();
205
gdouble SinGlat, CosGlat, SinGlon, CosGlon, Tau, lmst, x, y, z;
206
gdouble SinTau, CosTau, SinDec, CosDec;
207
208
c->UT = UT;
209
210
year = (gint) (date / 10000);
211
month = (gint) ((date - year * 10000) / 100);
212
day = (gint) (date - year * 10000 - month * 100);
213
c->year = year;
214
c->month = month;
215
c->day = day;
216
217
c->doy = DayofYear(year, month, day);
218
c->dow = DayofWeek(year, month, day, c->dowstr);
219
220
221
/*
222
* Compute Greenwich Mean Sidereal Time (gmst)
223
* The TU here is number of Julian centuries
224
* since 2000 January 1.5
225
* From the 1996 astronomical almanac
226
*/
227
TU = (jd(year, month, day, 0.0) - 2451545.0) / 36525.0;
228
TU2 = TU * TU;
229
TU3 = TU2 * TU;
230
T0 =
231
(6.0 + 41.0 / 60.0 + 50.54841 / 3600.0) +
232
8640184.812866 / 3600.0 * TU + 0.093104 / 3600.0 * TU2 -
233
6.2e-6 / 3600.0 * TU3;
234
T0 = hour24(T0);
235
236
c->gmst = hour24(T0 + UT * 1.002737909);
237
238
/* convert to radians for ease later on */
239
gmst = c->gmst * 15.0 * M_PI / 180.0;
240
241
lmst = 24.0 * frac((c->gmst - c->Glon / 15.0) / 24.0);
242
243
/*
244
245
* Construct Transformation Matrix from GEI to GSE systems
246
*
247
*
248
* First compute:
249
* mean ecliptic longitude of sun at epoch TU (varep)
250
* elciptic longitude of perigee at epoch TU (varpi)
251
* eccentricity of orbit at epoch TU (eccen)
252
*
253
* The TU here is the number of Julian centuries since
254
* 1900 January 0.0 (= 2415020.0)
255
*/
256
TDT = UT + 59.0 / 3600.0;
257
TU = (jd(year, month, day, TDT) - 2415020.0) / 36525.0;
258
varep =
259
(279.6966778 + 36000.76892 * TU +
260
0.0003025 * TU * TU) * RadPerDeg;
261
varpi =
262
(281.2208444 + 1.719175 * TU +
263
0.000452778 * TU * TU) * RadPerDeg;
264
eccen = 0.01675104 - 0.0000418 * TU - 0.000000126 * TU * TU;
265
c->eccentricity = eccen;
266
267
/*
268
* Compute the Obliquity of the Ecliptic at epoch TU
269
* The TU in this formula is the number of Julian
270
* centuries since epoch 2000 January 1.5
271
*/
272
TU = (jd(year, month, day, TDT) - jd(2000, 1, 1, 12.0)) / 36525.0;
273
epsilon = (23.43929167 - 0.013004166 * TU - 1.6666667e-7 * TU * TU
274
- 5.0277777778e-7 * TU * TU * TU) * RadPerDeg;
275
c->epsilon = epsilon;
276
277
/*
278
* Compute:
279
* Number of Days since epoch 1990.0 (days)
280
* The Mean Anomaly (M)
281
* The True Anomaly (nu)
282
* The Eccentric Anomaly via Keplers equation (E)
283
*
284
*
285
*/
286
days = jd(year, month, day, TDT) - jd(year, month, day, TDT);
287
M = angle2pi(2.0 * M_PI / 365.242191 * days + varep - varpi);
288
E = kepler(M, eccen);
289
nu =
290
2.0 * atan(sqrt((1.0 + eccen) / (1.0 - eccen)) * tan(E / 2.0));
291
lambnew = angle2pi(nu + varpi);
292
c->lambda_sun = lambnew;
293
294
/*
295
* Compute distance from earth to the sun
296
*/
297
r0 = 1.495985e8; /* in km */
298
earth_sun_distance =
299
r0 * (1 - eccen * eccen) / (1.0 + eccen * cos(nu)) / 6371.2;
300
c->earth_sun_dist = earth_sun_distance;
301
302
/*
303
* Compute Right Ascension and Declination of the Sun
304
*/
305
RA =
306
angle360(atan2(sin(lambnew) * cos(epsilon), cos(lambnew)) *
307
180.0 / M_PI);
308
DEC = asin(sin(epsilon) * sin(lambnew)) * 180.0 / M_PI;
309
c->RA_sun = RA;
310
c->DEC_sun = DEC;
311
312
/*
313
* Compute Moon Phase and AGE Stuff. The AGE that comes out of Moon()
314
* is actually the Phase converted to days. Since AGE is actually defined
315
* to be time since last NewMoon, we need to figure out what the JD of the
316
* last new moon was. Thats done below....
317
*/
318
TU = (jd(year, month, day, TDT) - 2451545.0) / 36525.0;
319
c->MoonPhase = Moon(TU, &LambdaMoon, &BetaMoon, &R, &AGE);
320
LambdaMoon *= RadPerDeg;
321
BetaMoon *= RadPerDeg;
322
323
324
RA_Moon =
325
angle360(atan2
326
(sin(LambdaMoon) * cos(epsilon) -
327
tan(BetaMoon) * sin(epsilon),
328
cos(LambdaMoon)) * DegPerRad);
329
DEC_Moon =
330
asin(sin(BetaMoon) * cos(epsilon) +
331
cos(BetaMoon) * sin(epsilon) * sin(LambdaMoon)) *
332
DegPerRad;
333
c->RA_moon = RA_Moon;
334
c->DEC_moon = DEC_Moon;
335
336
/*
337
* Compute Alt/Az coords
338
*/
339
Tau = (15.0 * lmst - RA_Moon) * RadPerDeg;
340
CosGlat = cos(c->Glat * RadPerDeg);
341
SinGlat = sin(c->Glat * RadPerDeg);
342
CosGlon = cos(c->Glon * RadPerDeg);
343
SinGlon = sin(c->Glon * RadPerDeg);
344
CosTau = cos(Tau);
345
SinTau = sin(Tau);
346
SinDec = sin(DEC_Moon * RadPerDeg);
347
CosDec = cos(DEC_Moon * RadPerDeg);
348
x = CosDec * CosTau * SinGlat - SinDec * CosGlat;
349
y = CosDec * SinTau;
350
z = CosDec * CosTau * CosGlat + SinDec * SinGlat;
351
c->A_moon = DegPerRad * atan2(y, x) + 180.0;
352
c->h_moon = DegPerRad * asin(z);
353
c->Visible = (c->h_moon < 0.0) ? 0 : 1;
354
355
/*
356
* Compute accurate AGE of the Moon
357
*/
358
Tb = TU - AGE / 36525.0; /* should be very close to minimum */
359
Ta = Tb - 0.4 / 36525.0;
360
Tc = Tb + 0.4 / 36525.0;
361
c->MoonAge = (TU - NewMoon(Ta, Tb, Tc)) * 36525.0;
362
363
/*
364
* Set Earth-Moon distance (calc above w/ func Moon)
365
*/
366
c->EarthMoonDistance = R;
367
368
369
c->SinGlat = SinGlat;
370
c->CosGlat = CosGlat;
371
}
Loggerhead is a web-based interface for Breezy
Version: 2.0.1