/* -*- 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 <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	135.27
#define LUNAR_NODE_MEAN_LONG	125.08
#define LUNAR_PROGRESSION	13.176358
#define LUNAR_INCLINATION	DEGREES_TO_RADIANS(5.145396)



/**
 * calc_moon:
 * @info:
 *
 * Sets "moonphase" expressed in degrees
 *   0 = New moon;   90 = First quarter
 * 180 = Full moon, 270 = Last quarter
 */

gboolean
calc_moon (WeatherInfo *info)
{
    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 = info->update;
    ndays = EPOCH_TO_J2000(t) / 86400.;
    sunMeanAnom_d = fmod ((ndays * SOL_PROGRESSION
			   + MEAN_ECLIPTIC_LONGITUDE (ndays)
			   - PERIGEE_LONGITUDE (ndays)),
			  360.);
    sunEclipLong_r = sunEclipLongitude (t);

    moonLong_d = fmod (LUNAR_MEAN_LONGITUDE + (ndays * LUNAR_PROGRESSION),
		       360.);
    moonMeanAnom_d = fmod ((moonLong_d - (0.1114041 * ndays)
			    - LUNAR_PERIGEE_MEAN_LONG),
			   360.);                          /*  5 */
    ascNodeMeanLong_d = fmod (LUNAR_NODE_MEAN_LONG - (0.0529539 * ndays),
			      360.);                       /*  6 */
    eviction_d = 1.2739
        * sin (DEGREES_TO_RADIANS (2.0 * (moonLong_d - RADIANS_TO_DEGREES (sunEclipLong_r))
				   - moonMeanAnom_d));          /*  7 */
    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" */

    /* 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: correct long of node */
    corrLong_d = 0.16 * sin (DEGREES_TO_RADIANS (sunMeanAnom_d));

                                               /* ( l'' - N' ) */
    lN_r = DEGREES_TO_RADIANS (moonLong_d + ascNodeMeanLong_d - corrLong_d);

    lambda_r = DEGREES_TO_RADIANS (corrLong_d) /* 16-19: ecliptic latitude */
        + atan2 (sin (lN_r) * cos (LUNAR_INCLINATION),  cos (lN_r));
                                               /* 20:   ecliptic longitude */
    beta_r = asin (sin (lN_r) * sin (LUNAR_INCLINATION));

    ecl2equ (t, lambda_r, beta_r, &ra_h, &decl_r);

    /*
     * The age or phase of the moon is based on the angle between the
     * longitude of the sun and moon
     */
    info->moonphase =
        fmod (15.*ra_h - RADIANS_TO_DEGREES (sunEclipLongitude (info->update)),
	      360.);
    if (info->moonphase < 0)
        info->moonphase += 360;
    info->moonlatitude = RADIANS_TO_DEGREES (decl_r);
    info->moonValid = TRUE;

    return TRUE;
}
