← Back to branch summary

~ubuntu-branches/ubuntu/lucid/libgweather/lucid

« back to all changes in this revision

Viewing changes to libgweather/weather-moon.c

  • Committer: Bazaar Package Importer
  • Author(s): Sebastien Bacher
  • Date: 2010-01-07 14:05:06 UTC
  • mfrom: (1.1.20 upstream)
  • Revision ID: james.westby@ubuntu.com-20100107140506-mz4gowhy70do2k4z
Tags: 2.29.4-0ubuntu1
* New upstream version
* debian/patches/lp291853_obsolete_timezones.patch,
  debian/patches/lp404616_default_station_for_montreal.patch:
  - the changes are in the new version
* debian/patches/99_autoreconf.patch:
  - new version update

Show diffs side-by-side

added added

removed removed

Lines of Context:
libgweather/weather-moon.c
 
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
 * <http://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
#define GWEATHER_I_KNOW_THIS_IS_UNSTABLE
 
39
#include "weather-priv.h"
 
40
 
 
41
/*
 
42
 * Elements of the Moon's orbit, epoch 2000 Jan 1.5
 
43
 * http://ssd.jpl.nasa.gov/?sat_elem#earth
 
44
 * The page only lists most values to 2 decimal places
 
45
 */
 
46
 
 
47
#define LUNAR_MEAN_LONGITUDE    218.316
 
48
#define LUNAR_PERIGEE_MEAN_LONG 318.15
 
49
#define LUNAR_NODE_MEAN_LONG    125.08
 
50
#define LUNAR_PROGRESSION       13.176358
 
51
#define LUNAR_INCLINATION       DEGREES_TO_RADIANS(5.145396)
 
52
 
 
53
/**
 
54
 * calc_moon:
 
55
 * @info:  WeatherInfo containing time_t of interest.  The
 
56
 *    values moonphase, moonlatitude and moonValid are updated
 
57
 *    on success.
 
58
 *
 
59
 * Returns: gboolean indicating success or failure.
 
60
 *    moonphase is expressed as degrees where '0' is a new moon,
 
61
 *    '90' is first quarter, etc.
 
62
 */
 
63
 
 
64
gboolean
 
65
calc_moon (WeatherInfo *info)
 
66
{
 
67
    time_t  t;
 
68
    gdouble ra_h;
 
69
    gdouble decl_r;
 
70
    gdouble ndays, sunMeanAnom_d;
 
71
    gdouble moonLong_d;
 
72
    gdouble moonMeanAnom_d, moonMeanAnom_r;
 
73
    gdouble sunEclipLong_r;
 
74
    gdouble ascNodeMeanLong_d;
 
75
    gdouble corrLong_d, eviction_d;
 
76
    gdouble sinSunMeanAnom;
 
77
    gdouble Ae, A3, Ec, A4, lN_r;
 
78
    gdouble lambda_r, beta_r;
 
79
 
 
80
    /*
 
81
     * The comments refer to the enumerated steps to calculate the
 
82
     * position of the moon (section 65 of above reference)
 
83
     */
 
84
    t = info->update;
 
85
    ndays = EPOCH_TO_J2000(t) / 86400.;
 
86
    sunMeanAnom_d = fmod (MEAN_ECLIPTIC_LONGITUDE (ndays) - PERIGEE_LONGITUDE (ndays),
 
87
                          360.);
 
88
    sunEclipLong_r = sunEclipLongitude (t);
 
89
    moonLong_d = fmod (LUNAR_MEAN_LONGITUDE + (ndays * LUNAR_PROGRESSION),
 
90
                       360.);
 
91
                                               /*  5: moon's mean anomaly */
 
92
    moonMeanAnom_d = fmod ((moonLong_d - (0.1114041 * ndays)
 
93
                            - (LUNAR_PERIGEE_MEAN_LONG + LUNAR_NODE_MEAN_LONG)),
 
94
                           360.);
 
95
                                               /*  6: ascending node mean longitude */
 
96
    ascNodeMeanLong_d = fmod (LUNAR_NODE_MEAN_LONG - (0.0529539 * ndays),
 
97
                              360.);
 
98
    eviction_d = 1.2739                        /*  7: eviction */
 
99
        * sin (DEGREES_TO_RADIANS (2.0 * (moonLong_d - RADIANS_TO_DEGREES (sunEclipLong_r))
 
100
                                   - moonMeanAnom_d));
 
101
    sinSunMeanAnom = sin (DEGREES_TO_RADIANS (sunMeanAnom_d));
 
102
    Ae = 0.1858 * sinSunMeanAnom;
 
103
    A3 = 0.37   * sinSunMeanAnom;              /*  8: annual equation    */
 
104
    moonMeanAnom_d += eviction_d - Ae - A3;    /*  9: "third correction" */
 
105
    moonMeanAnom_r = DEGREES_TO_RADIANS (moonMeanAnom_d);
 
106
    Ec = 6.2886 * sin (moonMeanAnom_r);        /* 10: equation of center */
 
107
    A4 = 0.214 * sin (2.0 * moonMeanAnom_r);   /* 11: "yet another correction" */
 
108
 
 
109
    /* Steps 12-14 give the true longitude after correcting for variation */
 
110
    moonLong_d += eviction_d + Ec - Ae + A4
 
111
        + (0.6583 * sin (2.0 * (moonMeanAnom_r - sunEclipLong_r)));
 
112
 
 
113
                                        /* 15: corrected longitude of node */
 
114
    corrLong_d = ascNodeMeanLong_d - 0.16 * sinSunMeanAnom;
 
115
 
 
116
    /*
 
117
     * Calculate ecliptic latitude (16-19) and longitude (20) of the moon,
 
118
     * then convert to right ascension and declination.
 
119
     */
 
120
    lN_r = DEGREES_TO_RADIANS (moonLong_d - corrLong_d);   /* l''-N' */
 
121
    lambda_r = DEGREES_TO_RADIANS(corrLong_d)
 
122
        + atan2 (sin (lN_r) * cos (LUNAR_INCLINATION),  cos (lN_r));
 
123
    beta_r = asin (sin (lN_r) * sin (LUNAR_INCLINATION));
 
124
    ecl2equ (t, lambda_r, beta_r, &ra_h, &decl_r);
 
125
 
 
126
    /*
 
127
     * The phase is the angle from the sun's longitude to the moon's 
 
128
     */
 
129
    info->moonphase =
 
130
        fmod (15.*ra_h - RADIANS_TO_DEGREES (sunEclipLongitude (info->update)),
 
131
              360.);
 
132
    if (info->moonphase < 0)
 
133
        info->moonphase += 360;
 
134
    info->moonlatitude = RADIANS_TO_DEGREES (decl_r);
 
135
    info->moonValid = TRUE;
 
136
 
 
137
    return TRUE;
 
138
}
 
139
 
 
140
 
 
141
/**
 
142
 * calc_moon_phases:
 
143
 * @info:   WeatherInfo containing the time_t of interest
 
144
 * @phases: An array of four time_t values that will hold the returned values.
 
145
 *    The values are estimates of the time of the next new, quarter, full and
 
146
 *    three-quarter moons.
 
147
 *
 
148
 * Returns: gboolean indicating success or failure
 
149
 */
 
150
 
 
151
gboolean
 
152
calc_moon_phases (WeatherInfo *info, time_t *phases)
 
153
{
 
154
    WeatherInfo temp;
 
155
    time_t      *ptime;
 
156
    int         idx;
 
157
    gdouble     advance;
 
158
    int         iter;
 
159
    time_t      dtime;
 
160
 
 
161
    g_return_val_if_fail (info != NULL &&
 
162
                          (info->moonValid || calc_moon (info)),
 
163
                          FALSE);
 
164
 
 
165
    ptime = phases;
 
166
    memset(&temp, 0, sizeof(WeatherInfo));
 
167
 
 
168
    for (idx = 0; idx < 4; idx++) {
 
169
        temp.update = info->update;
 
170
        temp.moonphase = info->moonphase;
 
171
 
 
172
        /*
 
173
         * First estimate on how far the moon needs to advance
 
174
         * to get to the required phase
 
175
         */
 
176
        advance = (idx * 90.) - info->moonphase;
 
177
        if (advance < 0.)
 
178
            advance += 360.;
 
179
 
 
180
        for (iter = 0; iter < 10; iter++) {
 
181
            /* Convert angle change (degrees) to dtime (seconds) */
 
182
            dtime = advance / LUNAR_PROGRESSION * 86400.;
 
183
            if ((dtime > -10) && (dtime < 10))
 
184
                break;
 
185
            temp.update += dtime;
 
186
            (void)calc_moon (&temp);
 
187
 
 
188
            if (idx == 0 && temp.moonphase > 180.) {
 
189
                advance = 360. - temp.moonphase;
 
190
            } else {
 
191
                advance = (idx * 90.) - temp.moonphase;
 
192
            }
 
193
        }
 
194
        *ptime++ = temp.update;
 
195
    }
 
196
 
 
197
    return TRUE;
 
198
}