1 /* -*- Mode: C; tab-width: 8; indent-tabs-mode: t; c-basic-offset: 4 -*- */
2 /* weather-moon.c - Lunar calculations for gweather
3 *
4 * This program is free software; you can redistribute it and/or
5 * modify it under the terms of the GNU General Public License as
6 * published by the Free Software Foundation; either version 2 of the
7 * License, or (at your option) any later version.
8 *
9 * This program is distributed in the hope that it will be useful, but
10 * WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 * General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with this program; if not, see
16 * <https://www.gnu.org/licenses/>.
17 */
18
19 /*
20 * Formulas from:
21 * "Practical Astronomy With Your Calculator" (3e), Peter Duffett-Smith
22 * Cambridge University Press 1988
23 */
24
25 #ifdef HAVE_CONFIG_H
26 # include <config.h>
27 #endif
28
29 #ifdef __FreeBSD__
30 #include <sys/types.h>
31 #endif
32
33 #include <math.h>
34 #include <time.h>
35 #include <string.h>
36 #include <glib.h>
37
38 #include "gweather-private.h"
39
40 /*
41 * Elements of the Moon's orbit, epoch 2000 Jan 1.5
42 * http://ssd.jpl.nasa.gov/?sat_elem#earth
43 * The page only lists most values to 2 decimal places
44 */
45
46 #define LUNAR_MEAN_LONGITUDE 218.316
47 #define LUNAR_PERIGEE_MEAN_LONG 318.15
48 #define LUNAR_NODE_MEAN_LONG 125.08
49 #define LUNAR_PROGRESSION 13.176358
50 #define LUNAR_INCLINATION DEGREES_TO_RADIANS(5.145396)
51
52 /*
53 * calc_moon_internal:
54 * @info: WeatherInfo containing time_t of interest. The
55 * values moonphase, moonlatitude and moonValid are updated
56 * on success.
57 *
58 * Returns: gboolean indicating success or failure.
59 * moonphase is expressed as degrees where '0' is a new moon,
60 * '90' is first quarter, etc.
61 */
62 static gboolean
calc_moon_internal(time_t update,gdouble * moonphase,gdouble * moonlatitude)63 calc_moon_internal (time_t update, gdouble *moonphase, gdouble *moonlatitude)
64 {
65 time_t t;
66 gdouble ra_h;
67 gdouble decl_r;
68 gdouble ndays, sunMeanAnom_d;
69 gdouble moonLong_d;
70 gdouble moonMeanAnom_d, moonMeanAnom_r;
71 gdouble sunEclipLong_r;
72 gdouble ascNodeMeanLong_d;
73 gdouble corrLong_d, eviction_d;
74 gdouble sinSunMeanAnom;
75 gdouble Ae, A3, Ec, A4, lN_r;
76 gdouble lambda_r, beta_r;
77
78 /*
79 * The comments refer to the enumerated steps to calculate the
80 * position of the moon (section 65 of above reference)
81 */
82 t = update;
83 ndays = EPOCH_TO_J2000(t) / 86400.;
84 sunMeanAnom_d = fmod (MEAN_ECLIPTIC_LONGITUDE (ndays) - PERIGEE_LONGITUDE (ndays),
85 360.);
86 sunEclipLong_r = sunEclipLongitude (t);
87 moonLong_d = fmod (LUNAR_MEAN_LONGITUDE + (ndays * LUNAR_PROGRESSION),
88 360.);
89 /* 5: moon's mean anomaly */
90 moonMeanAnom_d = fmod ((moonLong_d - (0.1114041 * ndays)
91 - (LUNAR_PERIGEE_MEAN_LONG + LUNAR_NODE_MEAN_LONG)),
92 360.);
93 /* 6: ascending node mean longitude */
94 ascNodeMeanLong_d = fmod (LUNAR_NODE_MEAN_LONG - (0.0529539 * ndays),
95 360.);
96 eviction_d = 1.2739 /* 7: eviction */
97 * sin (DEGREES_TO_RADIANS (2.0 * (moonLong_d - RADIANS_TO_DEGREES (sunEclipLong_r))
98 - moonMeanAnom_d));
99 sinSunMeanAnom = sin (DEGREES_TO_RADIANS (sunMeanAnom_d));
100 Ae = 0.1858 * sinSunMeanAnom;
101 A3 = 0.37 * sinSunMeanAnom; /* 8: annual equation */
102 moonMeanAnom_d += eviction_d - Ae - A3; /* 9: "third correction" */
103 moonMeanAnom_r = DEGREES_TO_RADIANS (moonMeanAnom_d);
104 Ec = 6.2886 * sin (moonMeanAnom_r); /* 10: equation of center */
105 A4 = 0.214 * sin (2.0 * moonMeanAnom_r); /* 11: "yet another correction" */
106
107 /* Steps 12-14 give the true longitude after correcting for variation */
108 moonLong_d += eviction_d + Ec - Ae + A4
109 + (0.6583 * sin (2.0 * (moonMeanAnom_r - sunEclipLong_r)));
110
111 /* 15: corrected longitude of node */
112 corrLong_d = ascNodeMeanLong_d - 0.16 * sinSunMeanAnom;
113
114 /*
115 * Calculate ecliptic latitude (16-19) and longitude (20) of the moon,
116 * then convert to right ascension and declination.
117 */
118 lN_r = DEGREES_TO_RADIANS (moonLong_d - corrLong_d); /* l''-N' */
119 lambda_r = DEGREES_TO_RADIANS(corrLong_d)
120 + atan2 (sin (lN_r) * cos (LUNAR_INCLINATION), cos (lN_r));
121 beta_r = asin (sin (lN_r) * sin (LUNAR_INCLINATION));
122 ecl2equ (t, lambda_r, beta_r, &ra_h, &decl_r);
123
124 /*
125 * The phase is the angle from the sun's longitude to the moon's
126 */
127 *moonphase =
128 fmod (15.*ra_h - RADIANS_TO_DEGREES (sunEclipLongitude (update)),
129 360.);
130 if (*moonphase < 0)
131 *moonphase += 360;
132 *moonlatitude = RADIANS_TO_DEGREES (decl_r);
133
134 return TRUE;
135 }
136
137 void
_gweather_info_ensure_moon(GWeatherInfo * info)138 _gweather_info_ensure_moon (GWeatherInfo *info)
139 {
140 if (!info->location.latlon_valid)
141 return;
142
143 if (!info->moonValid)
144 info->moonValid = calc_moon_internal (info->current_time,
145 &info->moonphase,
146 &info->moonlatitude);
147 }
148
149 /*
150 * calc_moon_phases:
151 * @info: WeatherInfo containing the time_t of interest
152 * @phases: An array of four time_t values that will hold the returned values.
153 * The values are estimates of the time of the next new, quarter, full and
154 * three-quarter moons.
155 *
156 * Returns: gboolean indicating success or failure
157 */
158 static gboolean
calc_moon_phases(GWeatherInfo * info,time_t * phases)159 calc_moon_phases (GWeatherInfo *info, time_t *phases)
160 {
161 time_t tmp_update;
162 gdouble tmp_moonphase;
163 gdouble tmp_moonlatitude;
164 time_t *ptime;
165 int idx;
166 gdouble advance;
167 int iter;
168 time_t dtime;
169
170 _gweather_info_ensure_moon (info);
171
172 ptime = phases;
173
174 for (idx = 0; idx < 4; idx++) {
175 tmp_update = info->current_time;
176 tmp_moonphase = info->moonphase;
177
178 /*
179 * First estimate on how far the moon needs to advance
180 * to get to the required phase
181 */
182 advance = (idx * 90.) - info->moonphase;
183 if (advance < 0.)
184 advance += 360.;
185
186 for (iter = 0; iter < 10; iter++) {
187 /* Convert angle change (degrees) to dtime (seconds) */
188 dtime = advance / LUNAR_PROGRESSION * 86400.;
189 if ((dtime > -10) && (dtime < 10))
190 break;
191 tmp_update += dtime;
192 (void)calc_moon_internal (tmp_update, &tmp_moonphase, &tmp_moonlatitude);
193
194 if (idx == 0 && tmp_moonphase > 180.) {
195 advance = 360. - tmp_moonphase;
196 } else {
197 advance = (idx * 90.) - tmp_moonphase;
198 }
199 }
200 *ptime++ = tmp_update;
201 }
202
203 return TRUE;
204 }
205
206 /**
207 * gweather_info_get_upcoming_moonphases:
208 * @info: a #GWeatherInfo containing the time_t of interest
209 * @phases: (out) (array fixed-size=4) (element-type gulong): An array of four
210 * time_t values that will hold the returned values.
211 * The values are estimates of the time of the next new, quarter, full and
212 * three-quarter moons.
213 *
214 * Returns: gboolean indicating success or failure
215 */
216 gboolean
gweather_info_get_upcoming_moonphases(GWeatherInfo * info,time_t * phases)217 gweather_info_get_upcoming_moonphases (GWeatherInfo *info, time_t *phases)
218 {
219 g_return_val_if_fail (GWEATHER_IS_INFO (info), FALSE);
220 g_return_val_if_fail (phases != NULL, FALSE);
221
222 return calc_moon_phases(info, phases);
223 }
224