/* -*- Mode: C; tab-width: 8; indent-tabs-mode: t; c-basic-offset: 4 -*- */

/* weather-sun.c - Astronomy 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

 * Unless otherwise noted, comments referencing "steps" are related to

 * the algorithm presented in section 49 of above

 */

#ifdef HAVE_CONFIG_H

#include <config.h>

#endif

#include <math.h>

#include <time.h>

#include <glib.h>

#include "gweather-private.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_sun (GWeatherInfo *info, time_t t)

{

    GWeatherInfoPrivate *priv;

    gdouble obsLat;

    gdouble obsLon;

    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 */

    priv = info->priv;

    obsLat = priv->location.latitude;

    obsLon = priv->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.;

    priv->midnightSun =

	(obsLat > (M_PI/2.-decl_midn)) || (obsLat < (-M_PI/2.-decl_midn));

    priv->polarNight =

	(obsLat > (M_PI/2.+decl_noon)) || (obsLat < (-M_PI/2.+decl_noon));

    if (priv->midnightSun || priv->polarNight) {

	priv->sunriseValid = priv->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)) {

	priv->sunriseValid = priv->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.;

    priv->sunriseValid = ((rise1 >= 0.) && (rise1 < 24.));

    priv->sunrise = (rise1 * 3600.) + lcl_midn;

    set1  = (set1 + dt - tt) * 0.9972695661;

    if (set1 < 0.)

	set1 += 24;

    else if (set1 >= 24.)

	set1 -= 24.;

    priv->sunsetValid = ((set1 >= 0.) && (set1 < 24.));

    priv->sunset = (set1 * 3600.) + lcl_midn;

    return (priv->sunriseValid || priv->sunsetValid);

}

void

_gweather_info_ensure_sun (GWeatherInfo *info)

{

    GWeatherInfoPrivate *priv;

    priv = info->priv;

    if (!priv->sunriseValid && !priv->sunsetValid)

	calc_sun (info, priv->current_time);

}

/**

 * 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

gweather_info_next_sun_event (GWeatherInfo *info)

{

    time_t    now = time (NULL);

    struct tm ltm;

    time_t    nxtEvent;

    GWeatherInfoPrivate *priv;

    priv = info->priv;

    g_return_val_if_fail (info != NULL, -1);

    _gweather_info_ensure_sun (info);

    /* 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 (priv->sunsetValid &&

	(priv->sunset > now) && (priv->sunset < nxtEvent))

	nxtEvent = priv->sunset;

    if (priv->sunriseValid &&

	(priv->sunrise > now) && (priv->sunrise < nxtEvent))

	nxtEvent = priv->sunrise;

    return (gint)(nxtEvent - now);

}