1 /* Adapted from "moontool.c" by John Walker, Release 2.5
2 **
3 ** Quoting from the original:
4 **
5 ** The algorithms used in this program to calculate the positions Sun and
6 ** Moon as seen from the Earth are given in the book "Practical Astronomy
7 ** With Your Calculator" by Peter Duffett-Smith, Second Edition,
8 ** Cambridge University Press, 1981. Ignore the word "Calculator" in the
9 ** title; this is an essential reference if you're interested in
10 ** developing software which calculates planetary positions, orbits,
11 ** eclipses, and the like. If you're interested in pursuing such
12 ** programming, you should also obtain:
13 **
14 ** "Astronomical Formulae for Calculators" by Jean Meeus, Third Edition,
15 ** Willmann-Bell, 1985. A must-have.
16 **
17 ** "Planetary Programs and Tables from -4000 to +2800" by Pierre
18 ** Bretagnon and Jean-Louis Simon, Willmann-Bell, 1986. If you want the
19 ** utmost (outside of JPL) accuracy for the planets, it's here.
20 **
21 ** "Celestial BASIC" by Eric Burgess, Revised Edition, Sybex, 1985. Very
22 ** cookbook oriented, and many of the algorithms are hard to dig out of
23 ** the turgid BASIC code, but you'll probably want it anyway.
24 **
25 ** See http://www.fourmilab.ch/moontool/
26 */
27
28 #include <stdlib.h>
29 #include <stdio.h>
30 #include <math.h>
31 #include <time.h>
32
33 #include "astro.h"
34
35
36 #define FALSE 0
37 #define TRUE 1
38
39 /* Astronomical constants */
40
41 #define epoch 2444238.5 /* 1980 January 0.0 */
42
43 /* Constants defining the Sun's apparent orbit */
44
45 #define elonge 278.833540 /* Ecliptic longitude of the Sun
46 at epoch 1980.0 */
47 #define elongp 282.596403 /* Ecliptic longitude of the Sun at
48 perigee */
49 #define eccent 0.016718 /* Eccentricity of Earth's orbit */
50 #define sunsmax 1.495985e8 /* Semi-major axis of Earth's orbit, km */
51 #define sunangsiz 0.533128 /* Sun's angular size, degrees, at
52 semi-major axis distance */
53
54 /* Elements of the Moon's orbit, epoch 1980.0 */
55
56 #define mmlong 64.975464 /* Moon's mean lonigitude at the epoch */
57 #define mmlongp 349.383063 /* Mean longitude of the perigee at the
58 epoch */
59 #define mlnode 151.950429 /* Mean longitude of the node at the
60 epoch */
61 #define minc 5.145396 /* Inclination of the Moon's orbit */
62 #define mecc 0.054900 /* Eccentricity of the Moon's orbit */
63 #define mangsiz 0.5181 /* Moon's angular size at distance a
64 from Earth */
65 #define msmax 384401.0 /* Semi-major axis of Moon's orbit in km */
66 #define mparallax 0.9507 /* Parallax at distance a from Earth */
67 #define synmonth 29.53058868 /* Synodic month (new Moon to new Moon) */
68 #define lunatbase 2423436.0 /* Base date for E. W. Brown's numbered
69 series of lunations (1923 January 16) */
70
71 /* Properties of the Earth */
72
73 #define earthrad 6378.16 /* Radius of Earth in kilometres */
74
75
76 #define PI 3.14159265358979323846 /* Assume not near black hole nor in
77 Tennessee */
78
79 /* Handy mathematical functions */
80
81 #define sgn(x) (((x) < 0) ? -1 : ((x) > 0 ? 1 : 0)) /* Extract sign */
82 #define abs(x) ((x) < 0 ? (-(x)) : (x)) /* Absolute val */
83 #define fixangle(a) ((a) - 360.0 * (floor((a) / 360.0))) /* Fix angle */
84 #define torad(d) ((d) * (PI / 180.0)) /* Deg->Rad */
85 #define todeg(d) ((d) * (180.0 / PI)) /* Rad->Deg */
86 #define dsin(x) (sin(torad((x)))) /* Sin from deg */
87 #define dcos(x) (cos(torad((x)))) /* Cos from deg */
88
89
90 /*
91 * UNIX_TO_JULIAN -- Convert internal Unix date/time to astronomical
92 * Julian time (i.e. Julian date plus day fraction, expressed as
93 * a double).
94 */
95 double
unix_to_julian(time_t t)96 unix_to_julian( time_t t )
97 {
98 return (double) t / 86400.0 + 2440587.4999996666666666666;
99 }
100
101
102 /*
103 * JYEAR -- Convert Julian date to year, month, day, which are
104 * returned via integer pointers to integers.
105 */
106 static void
jyear(double td,int * yy,int * mm,int * dd)107 jyear( double td, int* yy, int* mm, int* dd )
108 {
109 double j, d, y, m;
110
111 td += 0.5; /* Astronomical to civil */
112 j = floor(td);
113 j = j - 1721119.0;
114 y = floor(((4 * j) - 1) / 146097.0);
115 j = (j * 4.0) - (1.0 + (146097.0 * y));
116 d = floor(j / 4.0);
117 j = floor(((4.0 * d) + 3.0) / 1461.0);
118 d = ((4.0 * d) + 3.0) - (1461.0 * j);
119 d = floor((d + 4.0) / 4.0);
120 m = floor(((5.0 * d) - 3) / 153.0);
121 d = (5.0 * d) - (3.0 + (153.0 * m));
122 d = floor((d + 5.0) / 5.0);
123 y = (100.0 * y) + j;
124 if (m < 10.0)
125 m = m + 3;
126 else {
127 m = m - 9;
128 y = y + 1;
129 }
130 *yy = y;
131 *mm = m;
132 *dd = d;
133 }
134
135
136 /*
137 * MEANPHASE -- Calculates time of the mean new Moon for a given
138 * base date. This argument K to this function is
139 * the precomputed synodic month index, given by:
140 *
141 * K = (year - 1900) * 12.3685
142 *
143 * where year is expressed as a year and fractional year.
144 */
145 static double
meanphase(double sdate,double k)146 meanphase( double sdate, double k )
147 {
148 double t, t2, t3, nt1;
149
150 /* Time in Julian centuries from 1900 January 0.5 */
151 t = (sdate - 2415020.0) / 36525;
152 t2 = t * t; /* Square for frequent use */
153 t3 = t2 * t; /* Cube for frequent use */
154
155 nt1 = 2415020.75933 + synmonth * k
156 + 0.0001178 * t2
157 - 0.000000155 * t3
158 + 0.00033 * dsin(166.56 + 132.87 * t - 0.009173 * t2);
159
160 return nt1;
161 }
162
163
164 /*
165 * TRUEPHASE -- Given a K value used to determine the
166 * mean phase of the new moon, and a phase
167 * selector (0.0, 0.25, 0.5, 0.75), obtain
168 * the true, corrected phase time.
169 */
170 static double
truephase(double k,double pha)171 truephase( double k, double pha )
172 {
173 double t, t2, t3, pt, m, mprime, f;
174 int apcor = FALSE;
175
176 k += pha; /* Add phase to new moon time */
177 t = k / 1236.85; /* Time in Julian centuries from
178 1900 January 0.5 */
179 t2 = t * t; /* Square for frequent use */
180 t3 = t2 * t; /* Cube for frequent use */
181 pt = 2415020.75933 /* Mean time of phase */
182 + synmonth * k
183 + 0.0001178 * t2
184 - 0.000000155 * t3
185 + 0.00033 * dsin(166.56 + 132.87 * t - 0.009173 * t2);
186
187 m = 359.2242 /* Sun's mean anomaly */
188 + 29.10535608 * k
189 - 0.0000333 * t2
190 - 0.00000347 * t3;
191 mprime = 306.0253 /* Moon's mean anomaly */
192 + 385.81691806 * k
193 + 0.0107306 * t2
194 + 0.00001236 * t3;
195 f = 21.2964 /* Moon's argument of latitude */
196 + 390.67050646 * k
197 - 0.0016528 * t2
198 - 0.00000239 * t3;
199 if ((pha < 0.01) || (abs(pha - 0.5) < 0.01)) {
200
201 /* Corrections for New and Full Moon */
202
203 pt += (0.1734 - 0.000393 * t) * dsin(m)
204 + 0.0021 * dsin(2 * m)
205 - 0.4068 * dsin(mprime)
206 + 0.0161 * dsin(2 * mprime)
207 - 0.0004 * dsin(3 * mprime)
208 + 0.0104 * dsin(2 * f)
209 - 0.0051 * dsin(m + mprime)
210 - 0.0074 * dsin(m - mprime)
211 + 0.0004 * dsin(2 * f + m)
212 - 0.0004 * dsin(2 * f - m)
213 - 0.0006 * dsin(2 * f + mprime)
214 + 0.0010 * dsin(2 * f - mprime)
215 + 0.0005 * dsin(m + 2 * mprime);
216 apcor = TRUE;
217 } else if ((abs(pha - 0.25) < 0.01 || (abs(pha - 0.75) < 0.01))) {
218 pt += (0.1721 - 0.0004 * t) * dsin(m)
219 + 0.0021 * dsin(2 * m)
220 - 0.6280 * dsin(mprime)
221 + 0.0089 * dsin(2 * mprime)
222 - 0.0004 * dsin(3 * mprime)
223 + 0.0079 * dsin(2 * f)
224 - 0.0119 * dsin(m + mprime)
225 - 0.0047 * dsin(m - mprime)
226 + 0.0003 * dsin(2 * f + m)
227 - 0.0004 * dsin(2 * f - m)
228 - 0.0006 * dsin(2 * f + mprime)
229 + 0.0021 * dsin(2 * f - mprime)
230 + 0.0003 * dsin(m + 2 * mprime)
231 + 0.0004 * dsin(m - 2 * mprime)
232 - 0.0003 * dsin(2 * m + mprime);
233 if (pha < 0.5)
234 /* First quarter correction */
235 pt += 0.0028 - 0.0004 * dcos(m) + 0.0003 * dcos(mprime);
236 else
237 /* Last quarter correction */
238 pt += -0.0028 + 0.0004 * dcos(m) - 0.0003 * dcos(mprime);
239 apcor = TRUE;
240 }
241 if (!apcor) {
242 (void)fprintf (stderr,
243 "TRUEPHASE called with invalid phase selector.\n");
244 abort();
245 }
246 return pt;
247 }
248
249
250 /*
251 * PHASEHUNT5 -- Find time of phases of the moon which surround
252 * the current date. Five phases are found, starting
253 * and ending with the new moons which bound the
254 * current lunation.
255 */
256 void
phasehunt5(double sdate,double phases[5])257 phasehunt5( double sdate, double phases[5] )
258 {
259 double adate, k1, k2, nt1, nt2;
260 int yy, mm, dd;
261
262 adate = sdate - 45;
263
264 jyear(adate, &yy, &mm, &dd);
265 k1 = floor((yy + ((mm - 1) * (1.0 / 12.0)) - 1900) * 12.3685);
266
267 adate = nt1 = meanphase(adate, k1);
268 while (TRUE) {
269 adate += synmonth;
270 k2 = k1 + 1;
271 nt2 = meanphase(adate, k2);
272 if (nt1 <= sdate && nt2 > sdate)
273 break;
274 nt1 = nt2;
275 k1 = k2;
276 }
277 phases [0] = truephase (k1, 0.0);
278 phases [1] = truephase (k1, 0.25);
279 phases [2] = truephase (k1, 0.5);
280 phases [3] = truephase (k1, 0.75);
281 phases [4] = truephase (k2, 0.0);
282 }
283
284
285 /*
286 * PHASEHUNT2 -- Find time of phases of the moon which surround
287 * the current date. Two phases are found.
288 */
289 void
phasehunt2(double sdate,double phases[2],double which[2])290 phasehunt2( double sdate, double phases[2], double which[2] )
291 {
292 double phases5[5];
293
294 phasehunt5( sdate, phases5 );
295 phases[0] = phases5[0];
296 which[0] = 0.0;
297 phases[1] = phases5[1];
298 which[1] = 0.25;
299 if ( phases[1] <= sdate ) {
300 phases[0] = phases[1];
301 which[0] = which[1];
302 phases[1] = phases5[2];
303 which[1] = 0.5;
304 if ( phases[1] <= sdate ) {
305 phases[0] = phases[1];
306 which[0] = which[1];
307 phases[1] = phases5[3];
308 which[1] = 0.75;
309 if ( phases[1] <= sdate ) {
310 phases[0] = phases[1];
311 which[0] = which[1];
312 phases[1] = phases5[4];
313 which[1] = 0.0;
314 }
315 }
316 }
317 }
318
319
320 /*
321 * KEPLER -- Solve the equation of Kepler.
322 */
323 static double
kepler(double m,double ecc)324 kepler( double m, double ecc )
325 {
326 double e, delta;
327 #define EPSILON 1E-6
328
329 e = m = torad(m);
330 do {
331 delta = e - ecc * sin(e) - m;
332 e -= delta / (1 - ecc * cos(e));
333 } while (abs (delta) > EPSILON);
334 return e;
335 }
336
337
338 /*
339 * PHASE -- Calculate phase of moon as a fraction:
340 *
341 * The argument is the time for which the phase is requested,
342 * expressed as a Julian date and fraction. Returns the terminator
343 * phase angle as a percentage of a full circle (i.e., 0 to 1),
344 * and stores into pointer arguments the illuminated fraction of
345 * the Moon's disc, the Moon's age in days and fraction, the
346 * distance of the Moon from the centre of the Earth, and the
347 * angular diameter subtended by the Moon as seen by an observer
348 * at the centre of the Earth.
349 *
350 * pphase: Illuminated fraction
351 * mage: Age of moon in days
352 * dist: Distance in kilometres
353 * angdia: Angular diameter in degrees
354 * sudist: Distance to Sun
355 * suangdia: Sun's angular diameter
356 */
357 double
phase(double pdate,double * pphase,double * mage,double * dist,double * angdia,double * sudist,double * suangdia)358 phase( double pdate, double* pphase, double* mage, double* dist, double* angdia, double* sudist, double* suangdia )
359 {
360
361 double Day, N, M, Ec, Lambdasun, ml, MM, MN, Ev, Ae, A3, MmP,
362 mEc, A4, lP, V, lPP, NP, y, x, Lambdamoon, BetaM,
363 MoonAge, MoonPhase,
364 MoonDist, MoonDFrac, MoonAng, MoonPar,
365 F, SunDist, SunAng;
366
367 /* Calculation of the Sun's position */
368
369 Day = pdate - epoch; /* Date within epoch */
370 N = fixangle((360 / 365.2422) * Day); /* Mean anomaly of the Sun */
371 M = fixangle(N + elonge - elongp); /* Convert from perigee
372 co-ordinates to epoch 1980.0 */
373 Ec = kepler(M, eccent); /* Solve equation of Kepler */
374 Ec = sqrt((1 + eccent) / (1 - eccent)) * tan(Ec / 2);
375 Ec = 2 * todeg(atan(Ec)); /* True anomaly */
376 Lambdasun = fixangle(Ec + elongp); /* Sun's geocentric ecliptic
377 longitude */
378 /* Orbital distance factor */
379 F = ((1 + eccent * cos(torad(Ec))) / (1 - eccent * eccent));
380 SunDist = sunsmax / F; /* Distance to Sun in km */
381 SunAng = F * sunangsiz; /* Sun's angular size in degrees */
382
383
384 /* Calculation of the Moon's position */
385
386 /* Moon's mean longitude */
387 ml = fixangle(13.1763966 * Day + mmlong);
388
389 /* Moon's mean anomaly */
390 MM = fixangle(ml - 0.1114041 * Day - mmlongp);
391
392 /* Moon's ascending node mean longitude */
393 MN = fixangle(mlnode - 0.0529539 * Day);
394
395 /* Evection */
396 Ev = 1.2739 * sin(torad(2 * (ml - Lambdasun) - MM));
397
398 /* Annual equation */
399 Ae = 0.1858 * sin(torad(M));
400
401 /* Correction term */
402 A3 = 0.37 * sin(torad(M));
403
404 /* Corrected anomaly */
405 MmP = MM + Ev - Ae - A3;
406
407 /* Correction for the equation of the centre */
408 mEc = 6.2886 * sin(torad(MmP));
409
410 /* Another correction term */
411 A4 = 0.214 * sin(torad(2 * MmP));
412
413 /* Corrected longitude */
414 lP = ml + Ev + mEc - Ae + A4;
415
416 /* Variation */
417 V = 0.6583 * sin(torad(2 * (lP - Lambdasun)));
418
419 /* True longitude */
420 lPP = lP + V;
421
422 /* Corrected longitude of the node */
423 NP = MN - 0.16 * sin(torad(M));
424
425 /* Y inclination coordinate */
426 y = sin(torad(lPP - NP)) * cos(torad(minc));
427
428 /* X inclination coordinate */
429 x = cos(torad(lPP - NP));
430
431 /* Ecliptic longitude */
432 Lambdamoon = todeg(atan2(y, x));
433 Lambdamoon += NP;
434
435 /* Ecliptic latitude */
436 BetaM = todeg(asin(sin(torad(lPP - NP)) * sin(torad(minc))));
437
438 /* Calculation of the phase of the Moon */
439
440 /* Age of the Moon in degrees */
441 MoonAge = lPP - Lambdasun;
442
443 /* Phase of the Moon */
444 MoonPhase = (1 - cos(torad(MoonAge))) / 2;
445
446 /* Calculate distance of moon from the centre of the Earth */
447
448 MoonDist = (msmax * (1 - mecc * mecc)) /
449 (1 + mecc * cos(torad(MmP + mEc)));
450
451 /* Calculate Moon's angular diameter */
452
453 MoonDFrac = MoonDist / msmax;
454 MoonAng = mangsiz / MoonDFrac;
455
456 /* Calculate Moon's parallax */
457
458 MoonPar = mparallax / MoonDFrac;
459
460 *pphase = MoonPhase;
461 *mage = synmonth * (fixangle(MoonAge) / 360.0);
462 *dist = MoonDist;
463 *angdia = MoonAng;
464 *sudist = SunDist;
465 *suangdia = SunAng;
466 return fixangle(MoonAge) / 360.0;
467 }
468