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
218
219
220
221
222
223
224
225
226
227
|
/* -*- Mode: C; tab-width: 8; indent-tabs-mode: t; c-basic-offset: 4 -*- */
/* weather-moon.c - Lunar calculations for gweather
*
* This program is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License as
* published by the Free Software Foundation; either version 2 of the
* License, or (at your option) any later version.
*
* This program is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, see
* <http://www.gnu.org/licenses/>.
*/
/*
* Formulas from:
* "Practical Astronomy With Your Calculator" (3e), Peter Duffett-Smith
* Cambridge University Press 1988
*/
#ifdef HAVE_CONFIG_H
# include <config.h>
#endif
#ifdef __FreeBSD__
#include <sys/types.h>
#endif
#include <math.h>
#include <time.h>
#include <string.h>
#include <glib.h>
#define GWEATHER_I_KNOW_THIS_IS_UNSTABLE
#include "weather-priv.h"
/*
* Elements of the Moon's orbit, epoch 2000 Jan 1.5
* http://ssd.jpl.nasa.gov/?sat_elem#earth
* The page only lists most values to 2 decimal places
*/
#define LUNAR_MEAN_LONGITUDE 218.316
#define LUNAR_PERIGEE_MEAN_LONG 318.15
#define LUNAR_NODE_MEAN_LONG 125.08
#define LUNAR_PROGRESSION 13.176358
#define LUNAR_INCLINATION DEGREES_TO_RADIANS(5.145396)
/*
* calc_moon_internal:
* @info: WeatherInfo containing time_t of interest. The
* values moonphase, moonlatitude and moonValid are updated
* on success.
*
* Returns: gboolean indicating success or failure.
* moonphase is expressed as degrees where '0' is a new moon,
* '90' is first quarter, etc.
*/
static gboolean
calc_moon_internal (time_t update, gdouble *moonphase, gdouble *moonlatitude)
{
time_t t;
gdouble ra_h;
gdouble decl_r;
gdouble ndays, sunMeanAnom_d;
gdouble moonLong_d;
gdouble moonMeanAnom_d, moonMeanAnom_r;
gdouble sunEclipLong_r;
gdouble ascNodeMeanLong_d;
gdouble corrLong_d, eviction_d;
gdouble sinSunMeanAnom;
gdouble Ae, A3, Ec, A4, lN_r;
gdouble lambda_r, beta_r;
/*
* The comments refer to the enumerated steps to calculate the
* position of the moon (section 65 of above reference)
*/
t = update;
ndays = EPOCH_TO_J2000(t) / 86400.;
sunMeanAnom_d = fmod (MEAN_ECLIPTIC_LONGITUDE (ndays) - PERIGEE_LONGITUDE (ndays),
360.);
sunEclipLong_r = sunEclipLongitude (t);
moonLong_d = fmod (LUNAR_MEAN_LONGITUDE + (ndays * LUNAR_PROGRESSION),
360.);
/* 5: moon's mean anomaly */
moonMeanAnom_d = fmod ((moonLong_d - (0.1114041 * ndays)
- (LUNAR_PERIGEE_MEAN_LONG + LUNAR_NODE_MEAN_LONG)),
360.);
/* 6: ascending node mean longitude */
ascNodeMeanLong_d = fmod (LUNAR_NODE_MEAN_LONG - (0.0529539 * ndays),
360.);
eviction_d = 1.2739 /* 7: eviction */
* sin (DEGREES_TO_RADIANS (2.0 * (moonLong_d - RADIANS_TO_DEGREES (sunEclipLong_r))
- moonMeanAnom_d));
sinSunMeanAnom = sin (DEGREES_TO_RADIANS (sunMeanAnom_d));
Ae = 0.1858 * sinSunMeanAnom;
A3 = 0.37 * sinSunMeanAnom; /* 8: annual equation */
moonMeanAnom_d += eviction_d - Ae - A3; /* 9: "third correction" */
moonMeanAnom_r = DEGREES_TO_RADIANS (moonMeanAnom_d);
Ec = 6.2886 * sin (moonMeanAnom_r); /* 10: equation of center */
A4 = 0.214 * sin (2.0 * moonMeanAnom_r); /* 11: "yet another correction" */
/* Steps 12-14 give the true longitude after correcting for variation */
moonLong_d += eviction_d + Ec - Ae + A4
+ (0.6583 * sin (2.0 * (moonMeanAnom_r - sunEclipLong_r)));
/* 15: corrected longitude of node */
corrLong_d = ascNodeMeanLong_d - 0.16 * sinSunMeanAnom;
/*
* Calculate ecliptic latitude (16-19) and longitude (20) of the moon,
* then convert to right ascension and declination.
*/
lN_r = DEGREES_TO_RADIANS (moonLong_d - corrLong_d); /* l''-N' */
lambda_r = DEGREES_TO_RADIANS(corrLong_d)
+ atan2 (sin (lN_r) * cos (LUNAR_INCLINATION), cos (lN_r));
beta_r = asin (sin (lN_r) * sin (LUNAR_INCLINATION));
ecl2equ (t, lambda_r, beta_r, &ra_h, &decl_r);
/*
* The phase is the angle from the sun's longitude to the moon's
*/
*moonphase =
fmod (15.*ra_h - RADIANS_TO_DEGREES (sunEclipLongitude (update)),
360.);
if (*moonphase < 0)
*moonphase += 360;
*moonlatitude = RADIANS_TO_DEGREES (decl_r);
return TRUE;
}
void
_gweather_info_ensure_moon (GWeatherInfo *info)
{
GWeatherInfoPrivate *priv;
priv = info->priv;
if (!priv->moonValid)
priv->moonValid = calc_moon_internal (priv->current_time,
&priv->moonphase,
&priv->moonlatitude);
}
/*
* calc_moon_phases:
* @info: WeatherInfo containing the time_t of interest
* @phases: An array of four time_t values that will hold the returned values.
* The values are estimates of the time of the next new, quarter, full and
* three-quarter moons.
*
* Returns: gboolean indicating success or failure
*/
static gboolean
calc_moon_phases (GWeatherInfo *info, time_t *phases)
{
GWeatherInfoPrivate *priv;
time_t tmp_update;
gdouble tmp_moonphase;
gdouble tmp_moonlatitude;
time_t *ptime;
int idx;
gdouble advance;
int iter;
time_t dtime;
_gweather_info_ensure_moon (info);
priv = info->priv;
ptime = phases;
for (idx = 0; idx < 4; idx++) {
tmp_update = priv->current_time;
tmp_moonphase = priv->moonphase;
/*
* First estimate on how far the moon needs to advance
* to get to the required phase
*/
advance = (idx * 90.) - priv->moonphase;
if (advance < 0.)
advance += 360.;
for (iter = 0; iter < 10; iter++) {
/* Convert angle change (degrees) to dtime (seconds) */
dtime = advance / LUNAR_PROGRESSION * 86400.;
if ((dtime > -10) && (dtime < 10))
break;
tmp_update += dtime;
(void)calc_moon_internal (tmp_update, &tmp_moonphase, &tmp_moonlatitude);
if (idx == 0 && tmp_moonphase > 180.) {
advance = 360. - tmp_moonphase;
} else {
advance = (idx * 90.) - tmp_moonphase;
}
}
*ptime++ = tmp_update;
}
return TRUE;
}
/**
* gweather_info_get_upcoming_moonphases:
* @info: a #GWeatherInfo containing the time_t of interest
* @phases: (out) (array fixed-size=4) (element-type gulong): An array of four
* time_t values that will hold the returned values.
* The values are estimates of the time of the next new, quarter, full and
* three-quarter moons.
*
* Returns: gboolean indicating success or failure
*/
gboolean
gweather_info_get_upcoming_moonphases (GWeatherInfo *info, time_t *phases)
{
g_return_val_if_fail (GWEATHER_IS_INFO (info), FALSE);
g_return_val_if_fail (phases != NULL, FALSE);
return calc_moon_phases(info, phases);
}
|
