/*-
 * SPDX-License-Identifier: BSD-3-Clause
 *
 * Copyright (c) 2019-2020 The DragonFly Project.  All rights reserved.
 *
 * This code is derived from software contributed to The DragonFly Project
 * by Aaron LI <aly@aaronly.me>
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in
 *    the documentation and/or other materials provided with the
 *    distribution.
 * 3. Neither the name of The DragonFly Project nor the names of its
 *    contributors may be used to endorse or promote products derived
 *    from this software without specific, prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
 * FOR A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE
 * COPYRIGHT HOLDERS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
 * INCIDENTAL, SPECIAL, EXEMPLARY OR CONSEQUENTIAL DAMAGES (INCLUDING,
 * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
 * AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
 * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
 * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 * SUCH DAMAGE.
 *
 * Reference:
 * Calendrical Calculations, The Ultimate Edition (4th Edition)
 * by Edward M. Reingold and Nachum Dershowitz
 * 2018, Cambridge University Press
 */

#include <math.h>
#include <stdbool.h>
#include <stddef.h>
#include <stdio.h>
#include <stdlib.h>

#include "basics.h"
#include "gregorian.h"
#include "moon.h"
#include "sun.h"
#include "utils.h"


/*
 * Time from new moon to new moon (i.e., a lunation).
 * Ref: Sec.(14.6), Eq.(14.44)
 */
const double mean_synodic_month = 29.530588861;


/*
 * Argument data 1 used by 'nth_new_moon()'.
 * Ref: Sec.(14.6), Table(14.3)
 */
static const struct nth_new_moon_arg1 {
        double  v;
        int     w;
        int     x;
        int     y;
        int     z;
} nth_new_moon_data1[] = {
        { -0.40720, 0,  0, 1,  0 },
        {  0.17241, 1,  1, 0,  0 },
        {  0.01608, 0,  0, 2,  0 },
        {  0.01039, 0,  0, 0,  2 },
        {  0.00739, 1, -1, 1,  0 },
        { -0.00514, 1,  1, 1,  0 },
        {  0.00208, 2,  2, 0,  0 },
        { -0.00111, 0,  0, 1, -2 },
        { -0.00057, 0,  0, 1,  2 },
        {  0.00056, 1,  1, 2,  0 },
        { -0.00042, 0,  0, 3,  0 },
        {  0.00042, 1,  1, 0,  2 },
        {  0.00038, 1,  1, 0, -2 },
        { -0.00024, 1, -1, 2,  0 },
        { -0.00007, 0,  2, 1,  0 },
        {  0.00004, 0,  0, 2, -2 },
        {  0.00004, 0,  3, 0,  0 },
        {  0.00003, 0,  1, 1, -2 },
        {  0.00003, 0,  0, 2,  2 },
        { -0.00003, 0,  1, 1,  2 },
        {  0.00003, 0, -1, 1,  2 },
        { -0.00002, 0, -1, 1, -2 },
        { -0.00002, 0,  1, 3,  0 },
        {  0.00002, 0,  0, 4,  0 },
};

/*
 * Argument data 2 used by 'nth_new_moon()'.
 * Ref: Sec.(14.6), Table(14.4)
 */
static const struct nth_new_moon_arg2 {
        double  i;
        double  j;
        double  l;
} nth_new_moon_data2[] = {
        { 251.88,  0.016321, 0.000165 },
        { 251.83, 26.651886, 0.000164 },
        { 349.42, 36.412478, 0.000126 },
        {  84.66, 18.206239, 0.000110 },
        { 141.74, 53.303771, 0.000062 },
        { 207.14,  2.453732, 0.000060 },
        { 154.84,  7.306860, 0.000056 },
        {  34.52, 27.261239, 0.000047 },
        { 207.19,  0.121824, 0.000042 },
        { 291.34,  1.844379, 0.000040 },
        { 161.72, 24.198154, 0.000037 },
        { 239.56, 25.513099, 0.000035 },
        { 331.55,  3.592518, 0.000023 },
};

/*
 * Calculate the moment of the $n-th new moon after the new moon of
 * January 11, 1 (Gregorian), which is the first new moon after RD 0.
 * This function is centered upon January, 2000.
 * Ref: Sec.(14.6), Eq.(14.45)
 */
double
nth_new_moon(int n)
{
        int n0 = 24724;  /* Months from RD 0 until j2000 */
        int k = n - n0;  /* Months since j2000 */
        double nm = 1236.85;  /* Months per century */
        double c = k / nm;  /* Julian centuries */
        double j2000 = 0.5 + gregorian_new_year(2000);

        double coef_ap[] = { 5.09766, mean_synodic_month * nm,
                             0.00015437, -0.000000150, 0.00000000073 };
        double approx = j2000 + poly(c, coef_ap, nitems(coef_ap));

        double extra = 0.000325 * sin_deg(299.77 + 132.8475848 * c -
                                          0.009173 * c*c);

        double coef_sa[] = { 2.5534, 29.10535670 * nm,
                             -0.0000014, -0.00000011 };
        double coef_la[] = { 201.5643, 385.81693528 * nm,
                             0.0107582, 0.00001238, -0.000000058 };
        double coef_ma[] = { 160.7108, 390.67050284 * nm,
                             -0.0016118, -0.00000227, 0.000000011 };
        double coef_om[] = { 124.7746, -1.56375588 * nm,
                             0.0020672, 0.00000215 };
        double solar_anomaly = poly(c, coef_sa, nitems(coef_sa));
        double lunar_anomaly = poly(c, coef_la, nitems(coef_la));
        double moon_argument = poly(c, coef_ma, nitems(coef_ma));
        double omega = poly(c, coef_om, nitems(coef_om));
        double E = 1.0 - 0.002516 * c - 0.0000074 * c*c;

        double sum_c = 0.0;
        const struct nth_new_moon_arg1 *arg1;
        for (size_t i = 0; i < nitems(nth_new_moon_data1); i++) {
                arg1 = &nth_new_moon_data1[i];
                sum_c += arg1->v * pow(E, arg1->w) * sin_deg(
                                arg1->x * solar_anomaly +
                                arg1->y * lunar_anomaly +
                                arg1->z * moon_argument);
        }
        double correction = -0.00017 * sin_deg(omega) + sum_c;

        double additional = 0.0;
        const struct nth_new_moon_arg2 *arg2;
        for (size_t i = 0; i < nitems(nth_new_moon_data2); i++) {
                arg2 = &nth_new_moon_data2[i];
                additional += arg2->l * sin_deg(arg2->i + arg2->j * k);
        }

        double dt = approx + correction + extra + additional;
        return universal_from_dynamical(dt);
}

/*
 * Mean longitude of moon at moment given in Julian centuries $c.
 * Ref: Sec.(14.6), Eq.(14.49)
 */
static double
lunar_longitude_mean(double c)
{
        double coef[] = { 218.3164477, 481267.88123421, -0.0015786,
                          1.0 / 538841.0, -1.0 / 65194000.0 };
        return mod_f(poly(c, coef, nitems(coef)), 360);
}

/*
 * Elongation of moon (angular distance from Sun) at moment given in Julian
 * centuries $c.
 * Ref: Sec.(14.6), Eq.(14.50)
 */
static double
lunar_elongation(double c)
{
        double coef[] = { 297.8501921, 445267.1114034, -0.0018819,
                          1.0 / 545868.0, -1.0 / 113065000.0 };
        return mod_f(poly(c, coef, nitems(coef)), 360);
}

/*
 * Mean anomaly of Sun (angular distance from perihelion) at moment given in
 * Julian centuries $c.
 * Ref: Sec.(14.6), Eq.(14.51)
 */
static double
solar_anomaly(double c)
{
        double coef[] = { 357.5291092, 35999.0502909, -0.0001536,
                          1.0 / 24490000.0 };
        return mod_f(poly(c, coef, nitems(coef)), 360);
}

/*
 * Mean anomaly of moon (angular distance from perigee) at moment given in
 * Julian centuries $c.
 * Ref: Sec.(14.6), Eq.(14.52)
 */
static double
lunar_anomaly(double c)
{
        double coef[] = { 134.9633964, 477198.8675055, 0.0087414,
                          1.0 / 69699.0, -1.0 / 14712000.0 };
        return mod_f(poly(c, coef, nitems(coef)), 360);
}

/*
 * Moon's argument of latitude (the distance from the moon's node) at moment
 * given in Julian centuries $c.
 * Ref: Sec.(14.6), Eq.(14.53)
 */
static double
moon_node(double c)
{
        double coef[] = { 93.2720950, 483202.0175233, -0.0036539,
                          -1.0 / 3526000.0, 1.0 / 863310000.0 };
        return mod_f(poly(c, coef, nitems(coef)), 360);
}

/*
 * Argument data used by 'lunar_longitude()'.
 * Ref: Sec.(14.6), Table(14.5)
 */
static const struct lunar_longitude_arg {
        int     v;
        int     w;
        int     x;
        int     y;
        int     z;
} lunar_longitude_data[] = {
        { 6288774, 0,  0,  1,  0 },
        { 1274027, 2,  0, -1,  0 },
        {  658314, 2,  0,  0,  0 },
        {  213618, 0,  0,  2,  0 },
        { -185116, 0,  1,  0,  0 },
        { -114332, 0,  0,  0,  2 },
        {   58793, 2,  0, -2,  0 },
        {   57066, 2, -1, -1,  0 },
        {   53322, 2,  0,  1,  0 },
        {   45758, 2, -1,  0,  0 },
        {  -40923, 0,  1, -1,  0 },
        {  -34720, 1,  0,  0,  0 },
        {  -30383, 0,  1,  1,  0 },
        {   15327, 2,  0,  0, -2 },
        {  -12528, 0,  0,  1,  2 },
        {   10980, 0,  0,  1, -2 },
        {   10675, 4,  0, -1,  0 },
        {   10034, 0,  0,  3,  0 },
        {    8548, 4,  0, -2,  0 },
        {   -7888, 2,  1, -1,  0 },
        {   -6766, 2,  1,  0,  0 },
        {   -5163, 1,  0, -1,  0 },
        {    4987, 1,  1,  0,  0 },
        {    4036, 2, -1,  1,  0 },
        {    3994, 2,  0,  2,  0 },
        {    3861, 4,  0,  0,  0 },
        {    3665, 2,  0, -3,  0 },
        {   -2689, 0,  1, -2,  0 },
        {   -2602, 2,  0, -1,  2 },
        {    2390, 2, -1, -2,  0 },
        {   -2348, 1,  0,  1,  0 },
        {    2236, 2, -2,  0,  0 },
        {   -2120, 0,  1,  2,  0 },
        {   -2069, 0,  2,  0,  0 },
        {    2048, 2, -2, -1,  0 },
        {   -1773, 2,  0,  1, -2 },
        {   -1595, 2,  0,  0,  2 },
        {    1215, 4, -1, -1,  0 },
        {   -1110, 0,  0,  2,  2 },
        {    -892, 3,  0, -1,  0 },
        {    -810, 2,  1,  1,  0 },
        {     759, 4, -1, -2,  0 },
        {    -713, 0,  2, -1,  0 },
        {    -700, 2,  2, -1,  0 },
        {     691, 2,  1, -2,  0 },
        {     596, 2, -1,  0, -2 },
        {     549, 4,  0,  1,  0 },
        {     537, 0,  0,  4,  0 },
        {     520, 4, -1,  0,  0 },
        {    -487, 1,  0, -2,  0 },
        {    -399, 2,  1,  0, -2 },
        {    -381, 0,  0,  2, -2 },
        {     351, 1,  1,  1,  0 },
        {    -340, 3,  0, -2,  0 },
        {     330, 4,  0, -3,  0 },
        {     327, 2, -1,  2,  0 },
        {    -323, 0,  2,  1,  0 },
        {     299, 1,  1, -1,  0 },
        {     294, 2,  0,  3,  0 },
};

/*
 * Calculate the geocentric longitude of moon (in degrees) at moment $t.
 * Ref: Sec.(14.6), Eq.(14.48)
 */
double
lunar_longitude(double t)
{
        double c = julian_centuries(t);
        double nu = nutation(t);

        double L_prime = lunar_longitude_mean(c);
        double D = lunar_elongation(c);
        double M = solar_anomaly(c);
        double M_prime = lunar_anomaly(c);
        double F = moon_node(c);
        double E = 1.0 - 0.002516 * c - 0.0000074 * c*c;

        double sum = 0.0;
        const struct lunar_longitude_arg *arg;
        for (size_t i = 0; i < nitems(lunar_longitude_data); i++) {
                arg = &lunar_longitude_data[i];
                sum += arg->v * pow(E, abs(arg->x)) * sin_deg(
                                arg->w * D +
                                arg->x * M +
                                arg->y * M_prime +
                                arg->z * F);
        }
        double correction = sum / 1e6;

        double venus = sin_deg(119.75 + 131.849 * c) * 3958 / 1e6;
        double jupiter = sin_deg(53.09 + 479264.29 * c) * 318 / 1e6;
        double flat_earth = sin_deg(L_prime - F) * 1962 / 1e6;

        return mod_f((L_prime + correction + venus + jupiter +
                      flat_earth + nu), 360);
}

/*
 * Argument data used by 'lunar_latitude()'.
 * Ref: Sec.(14.6), Table(14.6)
 */
static const struct lunar_latitude_arg {
        int     v;
        int     w;
        int     x;
        int     y;
        int     z;
} lunar_latitude_data[] = {
        { 5128122, 0,  0,  0,  1 },
        {  280602, 0,  0,  1,  1 },
        {  277693, 0,  0,  1, -1 },
        {  173237, 2,  0,  0, -1 },
        {   55413, 2,  0, -1,  1 },
        {   46271, 2,  0, -1, -1 },
        {   32573, 2,  0,  0,  1 },
        {   17198, 0,  0,  2,  1 },
        {    9266, 2,  0,  1, -1 },
        {    8822, 0,  0,  2, -1 },
        {    8216, 2, -1,  0, -1 },
        {    4324, 2,  0, -2, -1 },
        {    4200, 2,  0,  1,  1 },
        {   -3359, 2,  1,  0, -1 },
        {    2463, 2, -1, -1,  1 },
        {    2211, 2, -1,  0,  1 },
        {    2065, 2, -1, -1, -1 },
        {   -1870, 0,  1, -1, -1 },
        {    1828, 4,  0, -1, -1 },
        {   -1794, 0,  1,  0,  1 },
        {   -1749, 0,  0,  0,  3 },
        {   -1565, 0,  1, -1,  1 },
        {   -1491, 1,  0,  0,  1 },
        {   -1475, 0,  1,  1,  1 },
        {   -1410, 0,  1,  1, -1 },
        {   -1344, 0,  1,  0, -1 },
        {   -1335, 1,  0,  0, -1 },
        {    1107, 0,  0,  3,  1 },
        {    1021, 4,  0,  0, -1 },
        {     833, 4,  0, -1,  1 },
        {     777, 0,  0,  1, -3 },
        {     671, 4,  0, -2,  1 },
        {     607, 2,  0,  0, -3 },
        {     596, 2,  0,  2, -1 },
        {     491, 2, -1,  1, -1 },
        {    -451, 2,  0, -2,  1 },
        {     439, 0,  0,  3, -1 },
        {     422, 2,  0,  2,  1 },
        {     421, 2,  0, -3, -1 },
        {    -366, 2,  1, -1,  1 },
        {    -351, 2,  1,  0,  1 },
        {     331, 4,  0,  0,  1 },
        {     315, 2, -1,  1,  1 },
        {     302, 2, -2,  0, -1 },
        {    -283, 0,  0,  1,  3 },
        {    -229, 2,  1,  1, -1 },
        {     223, 1,  1,  0, -1 },
        {     223, 1,  1,  0,  1 },
        {    -220, 0,  1, -2, -1 },
        {    -220, 2,  1, -1, -1 },
        {    -185, 1,  0,  1,  1 },
        {     181, 2, -1, -2, -1 },
        {    -177, 0,  1,  2,  1 },
        {     176, 4,  0, -2, -1 },
        {     166, 4, -1, -1, -1 },
        {    -164, 1,  0,  1, -1 },
        {     132, 4,  0,  1, -1 },
        {    -119, 1,  0, -1, -1 },
        {     115, 4, -1,  0, -1 },
        {     107, 2, -2,  0,  1 },
};

/*
 * Calculate the geocentric latitude of moon (in degrees) at moment $t.
 * Lunar latitude ranges from about -6 to 6 degress.
 * Ref: Sec.(14.6), Eq.(14.63)
 */
double
lunar_latitude(double t)
{
        double c = julian_centuries(t);

        double L_prime = lunar_longitude_mean(c);
        double D = lunar_elongation(c);
        double M = solar_anomaly(c);
        double M_prime = lunar_anomaly(c);
        double F = moon_node(c);
        double E = 1.0 - 0.002516 * c - 0.0000074 * c*c;

        double sum = 0.0;
        const struct lunar_latitude_arg *arg;
        for (size_t i = 0; i < nitems(lunar_latitude_data); i++) {
                arg = &lunar_latitude_data[i];
                sum += arg->v * pow(E, abs(arg->x)) * sin_deg(
                                arg->w * D +
                                arg->x * M +
                                arg->y * M_prime +
                                arg->z * F);
        }
        double beta = sum / 1e6;

        double venus = (sin_deg(119.75 + 131.849 * c + F) +
                        sin_deg(119.75 + 131.849 * c - F)) * 175 / 1e6;
        double flat_earth = (-2235 * sin_deg(L_prime) +
                             127 * sin_deg(L_prime - M_prime) -
                             115 * sin_deg(L_prime + M_prime)) / 1e6;
        double extra = sin_deg(313.45 + 481266.484 * c) * 382 / 1e6;

        return (beta + venus + flat_earth + extra);
}

/*
 * Argument data used by 'lunar_distance()'.
 * Ref: Sec.(14.6), Table(14.7)
 */
static const struct lunar_distance_arg {
        int     v;
        int     w;
        int     x;
        int     y;
        int     z;
} lunar_distance_data[] = {
        { -20905355, 0,  0,  1,  0 },
        {  -3699111, 2,  0, -1,  0 },
        {  -2955968, 2,  0,  0,  0 },
        {   -569925, 0,  0,  2,  0 },
        {     48888, 0,  1,  0,  0 },
        {     -3149, 0,  0,  0,  2 },
        {    246158, 2,  0, -2,  0 },
        {   -152138, 2, -1, -1,  0 },
        {   -170733, 2,  0,  1,  0 },
        {   -204586, 2, -1,  0,  0 },
        {   -129620, 0,  1, -1,  0 },
        {    108743, 1,  0,  0,  0 },
        {    104755, 0,  1,  1,  0 },
        {     10321, 2,  0,  0, -2 },
        {         0, 0,  0,  1,  2 },
        {     79661, 0,  0,  1, -2 },
        {    -34782, 4,  0, -1,  0 },
        {    -23210, 0,  0,  3,  0 },
        {    -21636, 4,  0, -2,  0 },
        {     24208, 2,  1, -1,  0 },
        {     30824, 2,  1,  0,  0 },
        {     -8379, 1,  0, -1,  0 },
        {    -16675, 1,  1,  0,  0 },
        {    -12831, 2, -1,  1,  0 },
        {    -10445, 2,  0,  2,  0 },
        {    -11650, 4,  0,  0,  0 },
        {     14403, 2,  0, -3,  0 },
        {     -7003, 0,  1, -2,  0 },
        {         0, 2,  0, -1,  2 },
        {     10056, 2, -1, -2,  0 },
        {      6322, 1,  0,  1,  0 },
        {     -9884, 2, -2,  0,  0 },
        {      5751, 0,  1,  2,  0 },
        {         0, 0,  2,  0,  0 },
        {     -4950, 2, -2, -1,  0 },
        {      4130, 2,  0,  1, -2 },
        {         0, 2,  0,  0,  2 },
        {     -3958, 4, -1, -1,  0 },
        {         0, 0,  0,  2,  2 },
        {      3258, 3,  0, -1,  0 },
        {      2616, 2,  1,  1,  0 },
        {     -1897, 4, -1, -2,  0 },
        {     -2117, 0,  2, -1,  0 },
        {      2354, 2,  2, -1,  0 },
        {         0, 2,  1, -2,  0 },
        {         0, 2, -1,  0, -2 },
        {     -1423, 4,  0,  1,  0 },
        {     -1117, 0,  0,  4,  0 },
        {     -1571, 4, -1,  0,  0 },
        {     -1739, 1,  0, -2,  0 },
        {         0, 2,  1,  0, -2 },
        {     -4421, 0,  0,  2, -2 },
        {         0, 1,  1,  1,  0 },
        {         0, 3,  0, -2,  0 },
        {         0, 4,  0, -3,  0 },
        {         0, 2, -1,  2,  0 },
        {      1165, 0,  2,  1,  0 },
        {         0, 1,  1, -1,  0 },
        {         0, 2,  0,  3,  0 },
        {      8752, 2,  0, -1, -2 },
};

/*
 * Calculate the distance to moon (in meters) at moment $t.
 * Ref: Sec.(14.6), Eq.(14.65)
 */
double
lunar_distance(double t)
{
        double c = julian_centuries(t);

        double D = lunar_elongation(c);
        double M = solar_anomaly(c);
        double M_prime = lunar_anomaly(c);
        double F = moon_node(c);
        double E = 1.0 - 0.002516 * c - 0.0000074 * c*c;

        double correction = 0.0;
        const struct lunar_distance_arg *arg;
        for (size_t i = 0; i < nitems(lunar_distance_data); i++) {
                arg = &lunar_distance_data[i];
                correction += arg->v * pow(E, abs(arg->x)) * cos_deg(
                                arg->w * D +
                                arg->x * M +
                                arg->y * M_prime +
                                arg->z * F);
        }

        return 385000560.0 + correction;
}

/*
 * Calculate the altitude of moon (in degrees) above the horizon at
 * location ($latitude, $longitude) and moment $t, ignoring parallax
 * and refraction.
 * NOTE: This calculates the geocentric altitude viewed from the center
 * of the Earth.
 * Ref: Sec.(14.6), Eq.(14.64)
 */
double
lunar_altitude(double t, double latitude, double longitude)
{
        double lambda = lunar_longitude(t);
        double beta = lunar_latitude(t);
        double alpha = right_ascension(t, beta, lambda);
        double delta = declination(t, beta, lambda);
        double theta = sidereal_from_moment(t);
        double H = mod_f(theta + longitude - alpha, 360);

        double v = (sin_deg(latitude) * sin_deg(delta) +
                    cos_deg(latitude) * cos_deg(delta) * cos_deg(H));
        return mod3_f(arcsin_deg(v), -180, 180);
}

/*
 * Parallax of moon at moment $t and location ($latitude, $longitude).
 * Ref: Sec.(14.6), Eq.(14.66)
 */
static double
lunar_parallax(double t, double latitude, double longitude)
{
        double geo = lunar_altitude(t, latitude, longitude);
        double distance = lunar_distance(t);
        /* Equatorial horizontal parallax of the moon */
        double sin_pi = 6378140.0 / distance;
        return arcsin_deg(sin_pi * cos_deg(geo));
}

/*
 * Calculate the topocentric altitude of moon viewed from the Earth surface
 * at location ($latitude, $longitude) and at moment $t.
 * Ref: Sec.(14.6), Eq.(14.67)
 */
static double
lunar_altitude_topocentric(double t, double latitude, double longitude)
{
        return (lunar_altitude(t, latitude, longitude) -
                lunar_parallax(t, latitude, longitude));
}

/*
 * Calculate the observed altitude of the upper limb of moon at location
 * ($latitude, $longitude) and moment $t.
 * Ref: Sec.(14.7), Eq.(14.82)
 */
double
lunar_altitude_observed(double t, const struct location *loc)
{
        double moon_radius = 16.0 / 60.0;  /* 16 arcminutes */
        return (lunar_altitude_topocentric(t, loc->latitude, loc->longitude) +
                refraction(loc->elevation) + moon_radius);
}

/*
 * Calculate the phase of the moon, defined as the difference in
 * longitudes of the Sun and moon, at the given moment $t.
 * Ref: Sec.(14.6), Eq.(14.56)
 */
double
lunar_phase(double t)
{
        double phi = mod_f(lunar_longitude(t) - solar_longitude(t), 360);

        /*
         * To check whether the above result conflicts with the time of
         * new moon as calculated by the more precise 'nth_new_moon()'
         */
        double t0 = nth_new_moon(0);
        int n = (int)lround((t - t0) / mean_synodic_month);
        double tn = nth_new_moon(n);
        double phi2 = mod_f((t - tn) / mean_synodic_month, 1) * 360;

        /* prefer the approximation based on the 'nth_new_moon()' moment */
        return (fabs(phi - phi2) > 180) ? phi2 : phi;
}

/*
 * Calculate the moment of the next time at or after the given moment $t
 * when the phase of the moon is $phi degree.
 * Ref: Sec.(14.6), Eq.(14.58)
 */
double
lunar_phase_atafter(double phi, double t)
{
        double rate = mean_synodic_month / 360.0;
        double phase = lunar_phase(t);
        double tau = t + rate * mod_f(phi - phase, 360);

        /* estimate range (within 2 days) */
        double a = (t > tau - 2) ? t : tau - 2;
        double b = tau + 2;

        return invert_angular(lunar_phase, phi, a, b);
}

/*
 * Calculate the moment of the new moon before the given moment $t.
 * Ref: Sec.(14.6), Eq.(14.46)
 */
double
new_moon_before(double t)
{
        double t0 = nth_new_moon(0);
        double phi = lunar_phase(t);
        int n = (int)lround((t - t0) / mean_synodic_month - phi / 360.0);

        int k = n - 1;
        double t1 = nth_new_moon(k);
        while (t1 < t) {
                k++;
                t1 = nth_new_moon(k);
        }

        return nth_new_moon(k-1);
}

/*
 * Calculate the moment of the new moon at or after the given moment $t.
 * Ref: Sec.(14.6), Eq.(14.47)
 */
double
new_moon_atafter(double t)
{
        double t0 = nth_new_moon(0);
        double phi = lunar_phase(t);
        int n = (int)lround((t - t0) / mean_synodic_month - phi / 360.0);

        double t1 = nth_new_moon(n);
        while (t1 < t) {
                n++;
                t1 = nth_new_moon(n);
        }

        return t1;
}

/*
 * Calculate the moment of moonrise in standard time on fixed date $rd
 * at location $loc.
 * NOTE: Return an NaN if no moonrise.
 * Ref: Sec.(14.7), Eq.(14.83)
 */
double
moonrise(int rd, const struct location *loc)
{
        double t = (double)rd - loc->zone;  /* universal time */
        bool waning = lunar_phase(t) > 180.0;
        /* lunar altitude at midnight */
        double alt = lunar_altitude_observed(t, loc);
        double offset = alt / (4.0 * (90.0 - fabs(loc->latitude)));

        /* approximate rising time */
        double t_approx = t;
        if (waning) {
                t_approx -= offset;
                if (offset > 0)
                        t_approx += 1;
        } else {
                t_approx += 0.5 + offset;
        }

        /* binary search to determine the rising time */
        const double eps = 30.0 / 3600 / 24;  /* accuracy of 30 seconds */
        double a = t_approx - 6.0/24.0;
        double b = t_approx + 6.0/24.0;
        double t_rise;
        do {
                t_rise = (a + b) / 2.0;
                if (lunar_altitude_observed(t_rise, loc) > 0)
                        b = t_rise;
                else
                        a = t_rise;
        } while (fabs(a - b) >= eps);

        if (t_rise < t + 1) {
                t_rise += loc->zone;  /* standard time */
                /* may be just before to midnight */
                return (t_rise > (double)rd) ? t_rise : (double)rd;
        } else {
                /* no moonrise on this day */
                return NAN;
        }
}

/*
 * Calculate the moment of moonset in standard time on fixed date $rd
 * at location $loc.
 * NOTE: Return an NaN if no moonset.
 * Ref: Sec.(14.7), Eq.(14.84)
 */
double
moonset(int rd, const struct location *loc)
{
        double t = (double)rd - loc->zone;  /* universal time */
        bool waxing = lunar_phase(t) < 180.0;
        /* lunar altitude at midnight */
        double alt = lunar_altitude_observed(t, loc);
        double offset = alt / (4.0 * (90.0 - fabs(loc->latitude)));

        /* approximate setting time */
        double t_approx = t;
        if (waxing) {
                t_approx += offset;
                if (offset <= 0)
                        t_approx += 1;
        } else {
                t_approx += 0.5 - offset;
        }

        /* binary search to determine the setting time */
        const double eps = 30.0 / 3600 / 24;  /* accuracy of 30 seconds */
        double a = t_approx - 6.0/24.0;
        double b = t_approx + 6.0/24.0;
        double t_set;
        do {
                t_set = (a + b) / 2.0;
                if (lunar_altitude_observed(t_set, loc) < 0)
                        b = t_set;
                else
                        a = t_set;
        } while (fabs(a - b) >= eps);

        if (t_set < t + 1) {
                t_set += loc->zone;  /* standard time */
                /* may be just before to midnight */
                return (t_set > (double)rd) ? t_set : (double)rd;
        } else {
                /* no moonrise on this day */
                return NAN;
        }
}

/**************************************************************************/

/*
 * Print moon information at the given moment $t (in standard time)
 * and events in the year.
 */
void
show_moon_info(double t, const struct location *loc)
{
        char buf[64];
        int rd = (int)floor(t);
        double t_u = t - loc->zone;  /* universal time */

        /*
         * Lunar phase
         */

        double eps = 1e-5;
        double phi = lunar_phase(t_u);
        bool waxing = (phi <= 180);
        bool crescent = (phi <= 90 || phi > 270);
        const char *phase_name;
        if (fabs(phi) < eps || fabs(phi - 360) < eps)
                phase_name = "New Moon";
        else if (fabs(phi - 90) < eps)
                phase_name = "First Quarter";
        else if (fabs(phi - 180) < eps)
                phase_name = "Full Moon";
        else if (fabs(phi - 270) < eps)
                phase_name = "Last Quarter";
        else
                phase_name = NULL;

        if (phase_name) {
                snprintf(buf, sizeof(buf), "%s", phase_name);
        } else {
                snprintf(buf, sizeof(buf), "%s %s",
                         (waxing ? "Waxing" : "Waning"),
                         (crescent ? "Crescent" : "Gibbous"));
        }
        printf("Moon phase: %.2lf° (%s)\n", phi, buf);

        /*
         * Moon position
         */
        double lon = lunar_longitude(t_u);
        double alt = lunar_altitude_observed(t_u, loc);
        printf("Moon position: %.4lf° (longitude), %.4lf° (altitude)\n",
               lon, alt);

        /*
         * Moon rise and set
         */

        double moments[2] = { moonrise(rd, loc), moonset(rd, loc) };
        const char *names[2] = { "Moonrise", "Moonset" };
        if (!isnan(moments[0]) && !isnan(moments[1]) &&
            moments[0] > moments[1]) {
                double t_tmp = moments[0];
                moments[0] = moments[1];
                moments[1] = t_tmp;
                const char *p = names[0];
                names[0] = names[1];
                names[1] = p;
        }
        for (size_t i = 0; i < nitems(moments); i++) {
                if (isnan(moments[i]))
                        snprintf(buf, sizeof(buf), "(null)");
                else
                        format_time(buf, sizeof(buf), moments[i]);
                printf("%-8s: %s\n", names[i], buf);
        }

        /*
         * Moon phases in the year
         */

        int year = gregorian_year_from_fixed(rd);
        struct date date = { year, 1, 1 };
        double t_begin = fixed_from_gregorian(&date) - loc->zone;
        date.year++;
        double t_end = fixed_from_gregorian(&date) - loc->zone;

        printf("\nLunar events in year %d:\n", year);
        printf("%19s   %19s   %19s   %19s\n",
               "New Moon", "First Quarter", "Full Moon", "Last Quarter");

        double t_newmoon = t_begin;
        while ((t_newmoon = new_moon_atafter(t_newmoon)) < t_end) {
                t_newmoon += loc->zone;  /* to standard time */
                gregorian_from_fixed((int)floor(t_newmoon), &date);
                format_time(buf, sizeof(buf), t_newmoon);
                printf("%d-%02d-%02d %s",
                       date.year, date.month, date.day, buf);

                /*
                 * first quarter, full moon, last quarter
                 */
                double t_event = t_newmoon;
                int phi_events[] = { 90, 180, 270 };
                for (size_t i = 0; i < nitems(phi_events); i++) {
                        t_event = lunar_phase_atafter(phi_events[i], t_event);
                        t_event += loc->zone;
                        gregorian_from_fixed((int)floor(t_event), &date);
                        format_time(buf, sizeof(buf), t_event);
                        printf("   %d-%02d-%02d %s",
                               date.year, date.month, date.day, buf);
                }
                printf("\n");

                /* go to the next new moon */
                t_newmoon += 28;
        }
}