← Back to branch summary

~ubuntu-branches/ubuntu/hardy/gnome-applets/hardy-updates

« back to all changes in this revision

Viewing changes to libgweather/weather-sun.c

  • Committer: Bazaar Package Importer
  • Author(s): Loic Minier
  • Date: 2008-01-15 15:16:45 UTC
  • mfrom: (1.1.25 upstream)
  • Revision ID: james.westby@ubuntu.com-20080115151645-7adqdys4ovi43xqs
Tags: 2.21.4-0ubuntu1
* New upstream development release.
  - Upstream news:
    . libgweather has been moved off to it's own module.
    . Accessx Status Applet:
       - Make bug buddy work for this applet (Kjartan Maraas).
    . Character Picker:
       - Documentation fixes.
    . CPU Frequency Monitor:
       - Fix memory leaks (kripkensteiner@gmail.com).
    . Keyboard Indicator:
       - The layout can now be printed from the preview (Sergey Udaltsov).
       - Require libgnomekbd 2.21.4.1 or later (Sergey Udaltsov).
    . Mixer:
       - Accelerators have been removed from the menu since they were.
         virtually never usable and now you won't be tempted to type Ctrl-T.
         in an unfortunate context. The normal under-score style menu.
         sortcuts are still available (Ted Gould).
    . Null Applet:
       - Make sure our example doesn't crash (Ray Strode).
    . Sticky Notes:
       - Better alignment in the UI (Christian Rose, Andrew Burton).
    . General:
       - Don't install documentation from obsolete applets (Kjartan Maraas).
    . Translations:
      ar, ca, es, et, eu, ga, he, nb, nn, oc, vi, sk, sv
    . Documentation Translations:
      ca, sv
  - Add a libgweather-dev (>= 2.21.1) build-dep.
  - Bump up libgnomekbdui-dev build-dep to >= 2.21.4.1.
  - Drop patches 05_gweather-dist-fixes, 06_gweather-dist-files,
    08_gweather_locations: the files which these patches used to dist and
    fix are now in libgweather.
  - Update autoreconf patch, 98_autoreconf.
  - Drop now useless gnome-applets-dev package.
  - Update gnome-applets' install file.
  - Drop DEB_DH_MAKESHLIBS_ARGS_gnome-applets.

Show diffs side-by-side

added added

removed removed

Lines of Context:
libgweather/weather-sun.c
1
 
/* $Id: weather-sun.c 10286 2007-07-12 08:52:56Z callum $ */
2
 
 
3
 
/*
4
 
 * Astronomy calculations for gweather
5
 
 *
6
 
 * Formulas from:
7
 
 * "Practical Astronomy With Your Calculator" (3e), Peter Duffett-Smith
8
 
 * Cambridge University Press 1988
9
 
 *
10
 
 * Planetary Mean Orbit parameters from http://ssd.jpl.nasa.gov/elem_planets.html,
11
 
 * converting longitudes from heliocentric to geocentric coordinates
12
 
 * epoch J2000 (2000 Jan 1 12:00:00 UTC)
13
 
 */
14
 
 
15
 
#ifdef HAVE_CONFIG_H
16
 
#  include <config.h>
17
 
#endif
18
 
 
19
 
#ifdef __FreeBSD__
20
 
#include <sys/types.h>
21
 
#endif
22
 
 
23
 
#include <math.h>
24
 
#include <time.h>
25
 
#include <glib.h>
26
 
 
27
 
#include "weather-priv.h"
28
 
 
29
 
#define EPOCH_TO_J2000(t)       (t-946728000)
30
 
#define MEAN_ECLIPTIC_LONGITUDE 280.46435
31
 
#define PERIGEE_LONGITUDE       282.94719
32
 
#define ECCENTRICITY            0.01671002
33
 
#define MEAN_ECLIPTIC_OBLIQUITY 
34
 
#define SOL_PROGRESSION         (360./365.242191)
35
 
 
36
 
/*
37
 
 * Convert ecliptic longitude (radians) to equitorial coordinates,
38
 
 * expressed as right ascension (hours) and declination (radians)
39
 
 * Assume ecliptic latitude is 0.0
40
 
 */
41
 
static void ecl2equ (gdouble time, gdouble eclipLon,
42
 
                     gdouble *ra, gdouble *decl)
43
 
{
44
 
    gdouble jc = EPOCH_TO_J2000(time) / (36525. * 86400.);
45
 
    gdouble mEclipObliq = DEGREES_TO_RADIANS(23.439291667
46
 
                                             - .013004166 * jc
47
 
                                             - 1.666667e-7 * jc * jc
48
 
                                             + 5.027777e-7 * jc * jc * jc);
49
 
    *ra = RADIANS_TO_HOURS(atan2 ( sin(eclipLon) * cos(mEclipObliq), cos(eclipLon)));
50
 
    if (*ra < 0.)
51
 
        *ra += 24.;
52
 
    *decl = asin ( sin(mEclipObliq) * sin(eclipLon) );
53
 
}
54
 
 
55
 
/*
56
 
 * Calculate rising and setting times of an object
57
 
 * based on it equitorial coordinates
58
 
 */
59
 
static void gstObsv (gdouble ra, gdouble decl,
60
 
                     gdouble obsLat, gdouble obsLon,
61
 
                     gdouble *rise, gdouble *set)
62
 
{
63
 
    double a = acos(-tan(obsLat) * tan(decl));
64
 
    double b;
65
 
 
66
 
    if (isnan(a) != 0) {
67
 
        *set = *rise = a;
68
 
        return;
69
 
    }
70
 
    a = RADIANS_TO_HOURS(a);
71
 
    b = 24. - a + ra;
72
 
    a += ra;
73
 
    a -= RADIANS_TO_HOURS(obsLon);
74
 
    b -= RADIANS_TO_HOURS(obsLon);
75
 
    if ((a = fmod(a, 24.)) < 0)
76
 
      a += 24.;
77
 
    if ((b = fmod(b, 24.)) < 0)
78
 
      b += 24.;
79
 
 
80
 
    *set = a;
81
 
    *rise = b;
82
 
}
83
 
 
84
 
 
85
 
static gdouble t0(time_t date)
86
 
{
87
 
    register gdouble t = ((gdouble)(EPOCH_TO_J2000(date) / 86400)) / 36525.0;
88
 
    gdouble t0 = fmod( 6.697374558 + 2400.051366 * t + 2.5862e-5 * t * t, 24.);
89
 
    if (t0 < 0.)
90
 
        t0 += 24.;
91
 
    return t0;
92
 
}
93
 
 
94
 
 
95
 
 
96
 
static gboolean calc_sun2 (time_t t, gdouble obsLat, gdouble obsLon,
97
 
          time_t *rise, time_t *set)
98
 
{
99
 
    time_t gm_midn;
100
 
    time_t lcl_midn;
101
 
    gdouble gm_hoff, ndays, meanAnom, eccenAnom, delta, lambda;
102
 
    gdouble ra1, ra2, decl1, decl2;
103
 
    gdouble rise1, rise2, set1, set2;
104
 
    gdouble tt, t00;
105
 
    gdouble x, u, dt;
106
 
 
107
 
    /* Approximate preceding local midnight at observer's longitude */
108
 
    gm_midn = t - (t % 86400);
109
 
    gm_hoff = floor((RADIANS_TO_DEGREES(obsLon) + 7.5) / 15.);
110
 
    lcl_midn = gm_midn - 3600. * gm_hoff;
111
 
    if (t - lcl_midn >= 86400)
112
 
        lcl_midn += 86400;
113
 
    else if (lcl_midn > t)
114
 
        lcl_midn -= 86400;
115
 
    
116
 
    /* 
117
 
     * Ecliptic longitude of the sun at midnight local time
118
 
     * Start with an estimate based on a fixed daily rate
119
 
     */
120
 
    ndays = EPOCH_TO_J2000(lcl_midn) / 86400.;
121
 
    meanAnom = DEGREES_TO_RADIANS(ndays * SOL_PROGRESSION
122
 
                                          + MEAN_ECLIPTIC_LONGITUDE - PERIGEE_LONGITUDE);
123
 
    
124
 
    /*
125
 
     * Approximate solution of Kepler's equation:
126
 
     * Find E which satisfies  E - e sin(E) = M (mean anomaly)
127
 
     */
128
 
    eccenAnom = meanAnom;
129
 
    
130
 
    while (1e-12 < fabs( delta = eccenAnom - ECCENTRICITY * sin(eccenAnom) - meanAnom))
131
 
    {
132
 
        eccenAnom -= delta / (1.- ECCENTRICITY * cos(eccenAnom));
133
 
    }
134
 
 
135
 
    /* Sun's longitude on the ecliptic */
136
 
    lambda = fmod( DEGREES_TO_RADIANS ( PERIGEE_LONGITUDE )
137
 
                   + 2. * atan( sqrt((1.+ECCENTRICITY)/(1.-ECCENTRICITY))
138
 
                                * tan ( eccenAnom / 2. ) ),
139
 
                   2. * M_PI);
140
 
    
141
 
    /* Calculate equitorial coordinates of sun at previous and next local midnights */
142
 
 
143
 
    ecl2equ (lcl_midn, lambda, &ra1, &decl1);
144
 
    ecl2equ (lcl_midn + 86400., lambda + DEGREES_TO_RADIANS(SOL_PROGRESSION), &ra2, &decl2);
145
 
    
146
 
    /* Convert to rise and set times assuming the earth were to stay in its position
147
 
     * in its orbit around the sun.  This will soon be followed by interpolating
148
 
     * between the two
149
 
     */
150
 
 
151
 
    gstObsv (ra1, decl1, obsLat, obsLon - (gm_hoff * M_PI / 12.), &rise1, &set1);
152
 
    gstObsv (ra2, decl2, obsLat, obsLon - (gm_hoff * M_PI / 12.), &rise2, &set2);
153
 
    
154
 
    /* TODO: include calculations for regions near the poles. */
155
 
    if (isnan(rise1) || isnan(rise2))
156
 
        return FALSE;
157
 
 
158
 
    if (rise2 < rise1) {
159
 
        rise2 += 24.;
160
 
    }
161
 
    if (set2 < set1) {
162
 
        set2 += 24.;
163
 
    }    
164
 
 
165
 
    tt = t0(lcl_midn);
166
 
    t00 = tt - (gm_hoff + RADIANS_TO_HOURS(obsLon)) * 1.002737909;
167
 
 
168
 
    if (t00 < 0.)
169
 
        t00 += 24.;
170
 
 
171
 
    if (rise1 < t00) {
172
 
        rise1 += 24.;
173
 
        rise2 += 24.;
174
 
    }
175
 
    if (set1 < t00) {
176
 
        set1  += 24.;
177
 
        set2  += 24.;
178
 
    }
179
 
    
180
 
    /*
181
 
     * Interpolate between the two
182
 
     */
183
 
    rise1 = (24.07 * rise1 - t00 * (rise2 - rise1)) / (24.07 + rise1 - rise2);
184
 
    set1 = (24.07 * set1 - t00 * (set2 - set1)) / (24.07 + set1 - set2);
185
 
 
186
 
    decl2 = (decl1 + decl2) / 2.;
187
 
    x = DEGREES_TO_RADIANS(0.830725);
188
 
    u = acos ( sin(obsLat) / cos(decl2) );
189
 
    dt =  RADIANS_TO_HOURS ( asin ( sin(x) / sin(u) )
190
 
                                     / cos(decl2) );
191
 
    
192
 
    rise1 = (rise1 - dt - tt) * 0.9972695661;
193
 
    set1  = (set1 + dt - tt) * 0.9972695661;
194
 
    if (rise1 < 0.)
195
 
      rise1 += 24;
196
 
    else if (rise1 >= 24.)
197
 
      rise1 -= 24.;
198
 
    if (set1 < 0.)
199
 
      set1 += 24;
200
 
    else if (set1 >= 24.)
201
 
      set1 -= 24.;
202
 
    
203
 
    *rise = (rise1 * 3600.) + lcl_midn;
204
 
    *set = (set1 * 3600.) + lcl_midn;
205
 
    return TRUE;
206
 
}
207
 
 
208
 
 
209
 
gboolean calc_sun (WeatherInfo *info)
210
 
{
211
 
  time_t now = time (NULL);
212
 
 
213
 
  return info->location->latlon_valid
214
 
    && calc_sun2(now, info->location->latitude, info->location->longitude,
215
 
           &info->sunrise, &info->sunset);
216
 
}
217
 
 
218
 
 
219
 
/*
220
 
 * Returns the interval, in seconds, until the next "sun event":
221
 
 *  - local midnight, when rise and set times are recomputed
222
 
 *  - next sunrise, when icon changes to daytime version
223
 
 *  - next sunset, when icon changes to nighttime version
224
 
 */
225
 
gint weather_info_next_sun_event(WeatherInfo *info)
226
 
{
227
 
  time_t    now = time(NULL);
228
 
  struct tm ltm;
229
 
  time_t    nxtEvent;
230
 
 
231
 
  if (!calc_sun(info))
232
 
    return -1;
233
 
 
234
 
  /* Determine when the next local midnight occurs */
235
 
  (void) localtime_r (&now, &ltm);
236
 
  ltm.tm_sec = 0;
237
 
  ltm.tm_min = 0;
238
 
  ltm.tm_hour = 0;
239
 
  ltm.tm_mday++;
240
 
  nxtEvent = mktime(&ltm);
241
 
 
242
 
  if (info->sunset > now && info->sunset < nxtEvent)
243
 
    nxtEvent = info->sunset;
244
 
  if (info->sunrise > now && info->sunrise < nxtEvent)
245
 
    nxtEvent = info->sunrise;
246
 
  return (gint)(nxtEvent - now);
247
 
}