/* -*- Mode: C; tab-width: 8; indent-tabs-mode: t; c-basic-offset: 4 -*- */
/* weather-sun.c - Astronomy calculations for mateweather
*
* 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
* .
*/
/*
* Formulas from:
* "Practical Astronomy With Your Calculator" (3e), Peter Duffett-Smith
* Cambridge University Press 1988
* Unless otherwise noted, comments referencing "steps" are related to
* the algorithm presented in section 49 of above
*/
#ifdef HAVE_CONFIG_H
#include
#endif
#include
#include
#include
#define MATEWEATHER_I_KNOW_THIS_IS_UNSTABLE
#include "weather-priv.h"
#define ECCENTRICITY(d) (0.01671123 - (d)/36525.*0.00004392)
/*
* Ecliptic longitude of the sun at specified time (UT)
* The algoithm is described in section 47 of Duffett-Smith
* Return value is in radians
*/
gdouble
sunEclipLongitude(time_t t)
{
gdouble ndays, meanAnom, eccenAnom, delta, e, longitude;
/*
* Start with an estimate based on a fixed daily rate
*/
ndays = EPOCH_TO_J2000(t) / 86400.;
meanAnom = DEGREES_TO_RADIANS(MEAN_ECLIPTIC_LONGITUDE(ndays)
- PERIGEE_LONGITUDE(ndays));
/*
* Approximate solution of Kepler's equation:
* Find E which satisfies E - e sin(E) = M (mean anomaly)
*/
eccenAnom = meanAnom;
e = ECCENTRICITY(ndays);
while (1e-12 < fabs( delta = eccenAnom - e * sin(eccenAnom) - meanAnom))
{
eccenAnom -= delta / (1.- e * cos(eccenAnom));
}
/*
* Earth's longitude on the ecliptic
*/
longitude = fmod( DEGREES_TO_RADIANS (PERIGEE_LONGITUDE(ndays))
+ 2. * atan (sqrt ((1.+e)/(1.-e))
* tan (eccenAnom / 2.)),
2. * M_PI);
if (longitude < 0.) {
longitude += 2 * M_PI;
}
return longitude;
}
static gdouble
ecliptic_obliquity (gdouble time)
{
gdouble jc = EPOCH_TO_J2000 (time) / (36525. * 86400.);
gdouble eclip_secs = (84381.448
- (46.84024 * jc)
- (59.e-5 * jc * jc)
+ (1.813e-3 * jc * jc * jc));
return DEGREES_TO_RADIANS(eclip_secs / 3600.);
}
/*
* Convert ecliptic longitude and latitude (radians) to equitorial
* coordinates, expressed as right ascension (hours) and
* declination (radians)
*/
void
ecl2equ (gdouble time,
gdouble eclipLon, gdouble eclipLat,
gdouble *ra, gdouble *decl)
{
gdouble mEclipObliq = ecliptic_obliquity(time);
if (ra) {
*ra = RADIANS_TO_HOURS (atan2 ((sin (eclipLon) * cos (mEclipObliq)
- tan (eclipLat) * sin(mEclipObliq)),
cos (eclipLon)));
if (*ra < 0.)
*ra += 24.;
}
if (decl) {
*decl = asin (( sin (eclipLat) * cos (mEclipObliq))
+ cos (eclipLat) * sin (mEclipObliq) * sin(eclipLon));
}
}
/*
* Calculate rising and setting times for an object
* based on it equitorial coordinates (section 33 & 15)
* Returned "rise" and "set" values are sideral times in hours
*/
static void
gstObsv (gdouble ra, gdouble decl,
gdouble obsLat, gdouble obsLon,
gdouble *rise, gdouble *set)
{
double a = acos (-tan (obsLat) * tan (decl));
double b;
if (isnan (a) != 0) {
*set = *rise = a;
return;
}
a = RADIANS_TO_HOURS (a);
b = 24. - a + ra;
a += ra;
a -= RADIANS_TO_HOURS (obsLon);
b -= RADIANS_TO_HOURS (obsLon);
if ((a = fmod (a, 24.)) < 0)
a += 24.;
if ((b = fmod (b, 24.)) < 0)
b += 24.;
*set = a;
*rise = b;
}
static gdouble
t0 (time_t date)
{
gdouble t = ((gdouble)(EPOCH_TO_J2000 (date) / 86400)) / 36525.0;
gdouble t0 = fmod (6.697374558 + 2400.051366 * t + 2.5862e-5 * t * t, 24.);
if (t0 < 0.)
t0 += 24.;
return t0;
}
static gboolean
calc_sun2 (WeatherInfo *info, time_t t)
{
gdouble obsLat = info->location->latitude;
gdouble obsLon = info->location->longitude;
time_t gm_midn;
time_t lcl_midn;
gdouble gm_hoff, lambda;
gdouble ra1, ra2;
gdouble decl1, decl2;
gdouble decl_midn, decl_noon;
gdouble rise1, rise2;
gdouble set1, set2;
gdouble tt, t00;
gdouble x, u, dt;
/* Approximate preceding local midnight at observer's longitude */
obsLat = info->location->latitude;
obsLon = info->location->longitude;
gm_midn = t - (t % 86400);
gm_hoff = floor ((RADIANS_TO_DEGREES (obsLon) + 7.5) / 15.);
lcl_midn = gm_midn - 3600. * gm_hoff;
if (t - lcl_midn >= 86400)
lcl_midn += 86400;
else if (lcl_midn > t)
lcl_midn -= 86400;
lambda = sunEclipLongitude (lcl_midn);
/*
* Calculate equitorial coordinates of sun at previous
* and next local midnights
*/
ecl2equ (lcl_midn, lambda, 0., &ra1, &decl1);
ecl2equ (lcl_midn + 86400.,
lambda + DEGREES_TO_RADIANS(SOL_PROGRESSION), 0.,
&ra2, &decl2);
/*
* If the observer is within the Arctic or Antarctic Circles then
* the sun may be above or below the horizon for the full day.
*/
decl_midn = MIN(decl1,decl2);
decl_noon = (decl1+decl2)/2.;
info->midnightSun =
(obsLat > (M_PI/2.-decl_midn)) || (obsLat < (-M_PI/2.-decl_midn));
info->polarNight =
(obsLat > (M_PI/2.+decl_noon)) || (obsLat < (-M_PI/2.+decl_noon));
if (info->midnightSun || info->polarNight) {
info->sunriseValid = info->sunsetValid = FALSE;
return FALSE;
}
/*
* Convert to rise and set times based positions for the preceding
* and following local midnights.
*/
gstObsv (ra1, decl1, obsLat, obsLon - (gm_hoff * M_PI / 12.), &rise1, &set1);
gstObsv (ra2, decl2, obsLat, obsLon - (gm_hoff * M_PI / 12.), &rise2, &set2);
/* TODO: include calculations for regions near the poles. */
if (isnan(rise1) || isnan(rise2)) {
info->sunriseValid = info->sunsetValid = FALSE;
return FALSE;
}
if (rise2 < rise1) {
rise2 += 24.;
}
if (set2 < set1) {
set2 += 24.;
}
tt = t0(lcl_midn);
t00 = tt - (gm_hoff + RADIANS_TO_HOURS(obsLon)) * 1.002737909;
if (t00 < 0.)
t00 += 24.;
if (rise1 < t00) {
rise1 += 24.;
rise2 += 24.;
}
if (set1 < t00) {
set1 += 24.;
set2 += 24.;
}
/*
* Interpolate between the two to get a rise and set time
* based on the sun's position at local noon (step 8)
*/
rise1 = (24.07 * rise1 - t00 * (rise2 - rise1)) / (24.07 + rise1 - rise2);
set1 = (24.07 * set1 - t00 * (set2 - set1)) / (24.07 + set1 - set2);
/*
* Calculate an adjustment value to account for parallax,
* refraction and the Sun's finite diameter (steps 9,10)
*/
decl2 = (decl1 + decl2) / 2.;
x = DEGREES_TO_RADIANS(0.830725);
u = acos ( sin(obsLat) / cos(decl2) );
dt = RADIANS_TO_HOURS ( asin ( sin(x) / sin(u) ) / cos(decl2) );
/*
* Subtract the correction value from sunrise and add to sunset,
* then (step 11) convert sideral times to UT
*/
rise1 = (rise1 - dt - tt) * 0.9972695661;
if (rise1 < 0.)
rise1 += 24;
else if (rise1 >= 24.)
rise1 -= 24.;
info->sunriseValid = ((rise1 >= 0.) && (rise1 < 24.));
info->sunrise = (rise1 * 3600.) + lcl_midn;
set1 = (set1 + dt - tt) * 0.9972695661;
if (set1 < 0.)
set1 += 24;
else if (set1 >= 24.)
set1 -= 24.;
info->sunsetValid = ((set1 >= 0.) && (set1 < 24.));
info->sunset = (set1 * 3600.) + lcl_midn;
return (info->sunriseValid || info->sunsetValid);
}
/**
* calc_sun_time:
* @info: #WeatherInfo structure containing the observer's latitude
* and longitude in radians, fills in the sunrise and sunset times.
* @t: time_t
*
* Returns: gboolean indicating if the results are valid.
*/
gboolean
calc_sun_time (WeatherInfo *info, time_t t)
{
return info->location->latlon_valid && calc_sun2 (info, t);
}
/**
* calc_sun:
* @info: #WeatherInfo structure containing the observer's latitude
* and longitude in radians, fills in the sunrise and sunset times.
*
* Returns: gboolean indicating if the results are valid.
*/
gboolean
calc_sun (WeatherInfo *info)
{
return calc_sun_time(info, time(NULL));
}
/**
* weather_info_next_sun_event:
* @info: #WeatherInfo structure
*
* Returns: the interval, in seconds, until the next "sun event":
* - local midnight, when rise and set times are recomputed
* - next sunrise, when icon changes to daytime version
* - next sunset, when icon changes to nighttime version
*/
gint
weather_info_next_sun_event (WeatherInfo *info)
{
time_t now = time (NULL);
struct tm ltm;
time_t nxtEvent;
g_return_val_if_fail (info != NULL, -1);
if (!calc_sun (info))
return -1;
/* Determine when the next local midnight occurs */
(void) localtime_r (&now, <m);
ltm.tm_sec = 0;
ltm.tm_min = 0;
ltm.tm_hour = 0;
ltm.tm_mday++;
nxtEvent = mktime (<m);
if (info->sunsetValid &&
(info->sunset > now) && (info->sunset < nxtEvent))
nxtEvent = info->sunset;
if (info->sunriseValid &&
(info->sunrise > now) && (info->sunrise < nxtEvent))
nxtEvent = info->sunrise;
return (gint)(nxtEvent - now);
}