/******************************************************************************!
 * \file sun.c
 * \author Sebastien Beaugrand
 * \sa http://beaugrand.chez.com/
 * \copyright CeCILL 2.1 Free Software license
 * \note Source: Astronomical algorithms - Jean Meeus - 1991
 ******************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "algos.h"

/******************************************************************************!
 * \fn main
 ******************************************************************************/
int main(int argc, char* argv[])
{
    int year;
    int month;
    int day;
    double lat;
    double lon;
    double h0;

    if (argc != 4 && argc != 5 && argc != 6) {
        printf("Usage: %s <YYYY-MM-DD> <latitude> <longitude> [h0]\n",
               argv[0]);
        return EXIT_FAILURE;
    }
    sscanf(argv[1], "%d-%d-%d", &year, &month, &day);
    sscanf(argv[2], "%lf", &lat);
    sscanf(argv[3], "%lf", &lon);
    // Geometric altitude of the center of the body
    // at the time of apparent rising or setting
    // 34/60 = 0.5667
    // 0.833 in https://www.esrl.noaa.gov/gmd/grad/solcalc
    if (argc == 5) {
        sscanf(argv[4], "%lf", &h0);
    } else {
        h0 = -deg(0, 34, 0);
    }
    double gmtoff = gmtOffset(year, month, day);

    double jd = julianDay(year, month, day);
    double t = julianTime(jd);
    double eps;
    double dec;
    double ra;
    double nutationInLongitude;
    sunApparentRightAscensionAndDeclination(jd, &ra, &dec,
                                            &nutationInLongitude, &eps);
    double theta0 = reduceAngle(sideralTime(jd, t));
    theta0 = apparentSideralTime(theta0, nutationInLongitude, eps);

    double H0 =
        DEG(acos((SIN(h0) - SIN(lat) * SIN(dec)) / (COS(lat) * COS(dec))));
    double m0 = (ra - lon - theta0) / 360;
    if (m0 < 0) {
        m0 += 1;
    } else if (m0 >= 1) {
        m0 -= 1;
    }
    double m1 = m0 - H0 / 360;
    if (m1 < 0) {
        m1 += 1;
    } else if (m1 >= 1) {
        m1 -= 1;
    }
    double m2 = m0 + H0 / 360;
    if (m2 < 0) {
        m2 += 1;
    } else if (m2 >= 1) {
        m2 -= 1;
    }
#   ifdef WITH_INTERPOLATION
    double a1;
    double a2 = ra;
    double a3;
    double d1;
    double d2 = dec;
    double d3;
    sunApparentRightAscensionAndDeclination(jd - 1, &a1, &d1,
                                            &nutationInLongitude, &eps);
    sunApparentRightAscensionAndDeclination(jd + 1, &a3, &d3,
                                            &nutationInLongitude, &eps);
    double t1 = reduceAngle(theta0 + 360.985647 * m1);
    double t2 = reduceAngle(theta0 + 360.985647 * m0);
    double t3 = reduceAngle(theta0 + 360.985647 * m2);
    // deltaT
    // https://www.iers.org/IERS/EN/DataProducts/EarthOrientationData/eop.html
    double deltaT = 32.184 + 37;  // 2018
    double n0 = m0 + deltaT / 86400;
    double n1 = m1 + deltaT / 86400;
    double n2 = m2 + deltaT / 86400;
    double a;
    double b;
    double c;
    a = a2 - a1;
    b = a3 - a2;
    c = b - a;
    a1 = a2 + n1 * (a + b + n1 * c) / 2;
    a2 = a2 + n0 * (a + b + n0 * c) / 2;
    a3 = a2 + n2 * (a + b + n2 * c) / 2;
    a = d2 - d1;
    b = d3 - d2;
    c = b - a;
    d1 = d2 + n1 * (a + b + n1 * c) / 2;
    d3 = d2 + n2 * (a + b + n2 * c) / 2;
    double H1 = t1 + lon - a1;
    double H2 = t2 + lon - a2;
    double H3 = t3 + lon - a3;
    double h1 = DEG(asin(SIN(lat) * SIN(d1) + COS(lat) * COS(d1) * COS(H1)));
    double h3 = DEG(asin(SIN(lat) * SIN(d3) + COS(lat) * COS(d3) * COS(H3)));
    m0 += -H2 / 360;
    m1 += (h1 - h0) / (360 * COS(d1) * COS(lat) * SIN(H1));
    m2 += (h3 - h0) / (360 * COS(d3) * COS(lat) * SIN(H3));
#   endif
    double hr0 = m0 * 24 + gmtoff;
    double hr1 = m1 * 24 + gmtoff;
    double hr2 = m2 * 24 + gmtoff;
    int mi0 = round((hr0 - ((int) hr0)) * 60);
    int mi1 = round((hr1 - ((int) hr1)) * 60);
    int mi2 = round((hr2 - ((int) hr2)) * 60);
    if (mi0 == 60) {
        mi0 = 0;
        hr0 += 1;
    }
    if (mi1 == 60) {
        mi1 = 0;
        hr1 += 1;
    }
    if (mi2 == 60) {
        mi2 = 0;
        hr2 += 1;
    }
    /*  */ if (argc == 6 && *argv[5] == 'h') {
        printf("%02d", (int) hr0);
    } else if (argc == 6 && *argv[5] == 'm') {
        printf("%02d", mi0);
    } else {
        printf("%02d:%02d - %02d:%02d\n", (int) hr1, mi1, (int) hr2, mi2);
    }

    return EXIT_SUCCESS;
}