/*
This source file is part of Scol
For the latest info, see http://www.scolring.org

Copyright (c) 2010 Stephane Bisaro, aka Iri <iri@irizone.net>

This program is free software; you can redistribute it and/or modify it under
the terms of the GNU Lesser 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 Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License along with
this program; if not, write to the Free Software Foundation, Inc., 59 Temple
Place - Suite 330, Boston, MA 02111-1307, USA, or go to
http://www.gnu.org/copyleft/lesser.txt

For others informations, please contact us from http://www.scolring.org/
*/

/*
$ IRI : A part of this code comes from the application "Ephemeride". http://www.irizone.net
I developped it under GNU GPL v3 license but i authorize its use in Scol only.
*/

#include "date_scol.h"


/* Return if a year is bissextile or not */
int scol_date_bissextile (int year)
{
    return ((year % 400 == 0) || ((year % 4==0) && (year % 100 != 0)));
}

/* Return the number of days from a month
By example, if month is 4 (april), it returns 121 (or 122 if leap year) */
int scol_date_nbdays (int year, int month)
{
    int nbjours = 0;

    switch (month)
    {
        case 12: nbjours += 30;
        case 11: nbjours += 31;
        case 10: nbjours += 30;
        case 9: nbjours += 31;
        case 8: nbjours += 31;
        case 7: nbjours += 30;
        case 6: nbjours += 31;
        case 5: nbjours += 30;
        case 4: nbjours += 31;
        case 3: if (scol_date_bissextile (year)) nbjours += 29; else nbjours += 28;
        case 2: nbjours += 31;
        default : nbjours += 0;
    }
    return nbjours;
}

int scol_date_quantieme (int year, int month, int day)
{
    return day + scol_date_nbdays (year, month);
}

/* Sunrise and sunset */

void scol_date_sun_conversion (struct Sun_Entree **entree, double d, double mn, int flag)
{
    if (flag == 1)
        (*entree)->longitude = (-1.0) * d + mn / 60.0;
    else
        (*entree)->latitude = d + mn / 60.0;
    return;
}

void scol_date_sun_quantieme (struct Sun_Entree **entree)
{
    int nbjoursfrom1mars = -31;

    if (scol_date_bissextile ((*entree)->annee))
        nbjoursfrom1mars -= 29;
    else
        nbjoursfrom1mars -= 28;

    if ((*entree)->mois >= 3)
        nbjoursfrom1mars *= 1;
    else
        nbjoursfrom1mars += 365;

    (*entree)->fjour = (double) ((*entree)->jour + scol_date_nbdays ((*entree)->annee, (*entree)->mois) + nbjoursfrom1mars - 1 + (*entree)->uheure / 24);
    return;
}

double scol_date_sun_equation (struct Sun_Entree ** entree, double flag)
{
    /* heure du milieu de journée (TU) + equation du temps +/- la moitié de la durée d'ensoleillement du jour.
    lever : flag = -1;
    coucher : flag = +1 */
    return (*entree)->uheure + (((*entree)->soleil.eqtemps + (flag * (*entree)->soleil.anghoraire)) / Sun_E_hr);
}











/**
 * \brief Return if any year is a leap year or not
 *
 * Stack :
 * Stage 0 : I : any year : should be posterior or equal at 1582 (gregorian calendar)
 * return : I : 1 if the year is a leap year, 0 if not and nil if error
 */
int SCOLscienceDateLeapYear (mmachine m)
{
    int myear;

    MMechostr (MSKDEBUG, "SCOLscienceDateLeapYear : entering\n");

    myear = MTOI (MMpull (m));

    if (myear < 1582)
    {
        MMechostr (MSKDEBUG, "SCOLscienceDateLeapYear error : invalid year (must be >= 1582)\n");
        MMpush (m, NIL);
        return 0;
    }

    MMpush (m, ITOM (scol_date_bissextile (myear)));
    return 0;
}

/**
 * \brief Get the number of days from the 1st january
 *
 * Stack :
 * Stage 0 : I : the day (from 1 to 28, 29, 30 or 31)
 * Stage 1 : I : the month (from 1 to 12)
 * Stage 2 : I : the year (from 1582)
 * return : I : this number or nil if error
 */
int SCOLscienceDateDayNumber (mmachine m)
{
    int myear, mmonth, mday;

    MMechostr (MSKDEBUG, "SCOLscienceDateDayNumber : entering\n");

    mday = MTOI (MMpull (m));
    mmonth = MTOI (MMpull (m));
    myear = MTOI (MMpull (m));

    if (myear < 1582)
    {
        MMechostr (MSKDEBUG, "SCOLscienceDateLeapYear error : invalid year (must be >= 1582)\n");
        MMpush (m, NIL);
        return 0;
    }

    if ((mmonth < 1) || (mmonth > 12) ||(mday < 1) ||(mday > 31))
    {
        MMechostr (MSKDEBUG, "SCOLscienceDateLeapYear error : invalid datas\n");
        MMpush (m, NIL);
        return 0;
    }

    MMpush (m, ITOM (scol_date_quantieme (myear, mmonth, mday)));
    return 0;
}


/**
 * \brief Return the sunrise and sunset, anywhere
 *
 * Stack :
 * Stage 0 : [I I I] : date : day (1 - 31), month (1-12), year (>= 1582). Can be nil (if in this case, get the current date)
 * Stage 1 : [I I] : longitude : degree (0-180),minute (0-59). Positive to east, negative to west.
 * Stage 0 : [I I] : latitude : degree (0-180), minute (0 - 59). Positive to north, negative to south
 * Return [S S] : sunrise and sunset, at UTC, "--" if not sunrise / sunset, nil if error
 * example : ["7:21" "16:50"]
 */
int SCOLscienceDateSun (mmachine m)
{
    int mlong, mlat, mdate;
    int mlongd, mlongmn, mlatd, mlatmn, myear, mmonth, mday;
    double xrect, yrect, cs, P, mn, h;

    struct Sun_Entree *entree;
    time_t mt;
    struct tm *t;

    MMechostr (MSKDEBUG, "SCOLscienceDateSun : entering\n");

    mt = time (NULL);
    t = gmtime (&mt);

    entree = malloc (sizeof (struct Sun_Entree));
    if (entree == NULL)
    {
        MMechostr (MSKDEBUG, "SCOLscienceDateSun error : not enough memory\n");
        MMpush (m, NIL);
        return 0;
    }

    mdate = MTOP (MMpull (m));
    mlong = MTOP (MMpull (m));
    mlat = MTOP (MMpull (m));

    mlongd = MTOI (MMfetch (m, mlong, 0));
    mlongmn = MTOI (MMfetch (m, mlong, 1));

    mlatd = MTOI (MMfetch (m, mlat, 0));
    mlatmn = MTOI (MMfetch (m, mlat, 1));

    mday = MTOI (MMfetch (m, mdate, 0));
    mmonth = MTOI (MMfetch (m, mdate, 1));
    myear = MTOI (MMfetch (m, mdate, 2));

    /*

    Initialisation

    */


    /* Year */
    if (myear == NIL)
        entree->annee = 1900+t->tm_year;
    else if (myear < 1582)
    {
        MMechostr (MSKDEBUG, "SCOLscienceDateSun error : invalid year (must be >= 1582)\n");
        MMpush (m, NIL);
        free (entree);
        return 0;
    }
    else
        entree->annee = myear;

    /* Month and day */
    if (mmonth == NIL)
        entree->mois = 1.0+t->tm_mon;
    else if (mday == NIL)
        entree->jour = t->tm_mday;
    if ((mmonth < 1) || (mmonth > 12) ||(mday < 1) ||(mday > 31))
    {
        MMechostr (MSKDEBUG, "SCOLscienceDateSun error : invalid datas\n");
        MMpush (m, NIL);
        free (entree);
        return 0;
    }
    else
    {
        entree->mois = (double) mmonth;
        entree->jour = (double) mday;
    }

    /* Latitude and longitude */
    if ((mlatd == NIL) || (mlatmn == NIL))
        entree->latitude = 48.833;  /* Paris */
    else if ((mlatd > 90) || (mlatd < (-90)) || (mlatmn >= 60) || (mlatmn < 0))
    {
        MMechostr (MSKDEBUG, "SCOLscienceDateSun error : invalid latitude\n");
        MMpush (m, NIL);
        free (entree);
        return 0;
    }
    else
        scol_date_sun_conversion (&entree, (double) mlatd, (double) mlatmn, 0);

    if ((mlongd == NIL) || (mlongmn == NIL))
        entree->longitude = -2.333; /* Paris */
    else if ((mlongd > 90) || (mlongd < (-90)) || (mlongmn >= 60) || (mlongmn < 0))
    {
        MMechostr (MSKDEBUG, "SCOLscienceDateSun error : invalid longitude\n");
        MMpush (m, NIL);
        free (entree);
        return 0;
    }
    else
        scol_date_sun_conversion (&entree, (double) mlongd, (double) mlongmn, 1);

    /* coordonnées en radians */
    entree->longitude *= (M_PI / 180.0);
    /*(*entree)->latitude = (M_PI / 180.0) * (*entree)->latitude;*/
    entree->latitude *= (M_PI / 180.0);
    /* Heure du zénith en Temps Universel */
    entree->uheure = 12.0 + (entree->longitude / Sun_E_hr);
    /* Jour écoulé depuis le 1er mars */
    scol_date_sun_quantieme (&entree);


    /*

    Sun

    */

    /* 0 si le Soleil se lève ET se couche ce jour, 1 sinon */
    entree->soleil.solnon = 0;
    /* Anomalie moyenne */
    entree->soleil.moyanomalie = Sun_E_k * (entree->fjour - Sun_E_jm);
    /* Longitude moyenne */
    entree->soleil.moylongitude = Sun_E_k * (entree->fjour - Sun_E_jl);
    /* Longitude vraie */
    entree->soleil.vraielongitude = entree->soleil.moylongitude + (2.0 * Sun_E_e * sin (entree->soleil.moyanomalie)) +
                                        (1.25 * pow (Sun_E_e, 2.0) * sin ( 2.0 * entree->soleil.moyanomalie));

    /* Coordonnées rectangulaires 1 */
    entree->soleil.rectX = cos (entree->soleil.vraielongitude);
    entree->soleil.rectY = cos (Sun_E_ob) * sin (entree->soleil.vraielongitude);
    entree->soleil.rectZ = sin (Sun_E_ob) * sin (entree->soleil.vraielongitude);
    /* Coordonnées rectangulaires 2 */
    xrect = (cos (entree->soleil.moylongitude) * entree->soleil.rectX) + (sin (entree->soleil.moylongitude) * entree->soleil.rectY);
    yrect = ((-1.0) * sin (entree->soleil.moylongitude) * entree->soleil.rectX) + (cos (entree->soleil.moylongitude) * entree->soleil.rectY);
    entree->soleil.rectX = xrect;
    entree->soleil.rectY = yrect;
    /* Équation du temps */
    entree->soleil.eqtemps = atan2 (entree->soleil.rectY, entree->soleil.rectX);
    /* Déclinaison solaire */
    entree->soleil.declinaison = atan2 (entree->soleil.rectZ, (sqrt (1 - pow (entree->soleil.rectZ, 2.0))));
    /* Angle horaire */
    cs = (sin (Sun_E_ht) - (sin (entree->latitude) * sin (entree->soleil.declinaison))) / cos (entree->latitude) / cos (entree->soleil.declinaison);
    if (cs > 1.0)
        entree->soleil.solnon = 1; /* soleil ne se lève pas */
    if (cs < -1.0)
        entree->soleil.solnon = 1; /* soleil ne se couche pas */
    if ((!(equaldouble (cs, 0.0))))
        entree->soleil.anghoraire = M_PI / 2.0;
    else
        entree->soleil.anghoraire = atan (sqrt (1.0 - pow (cs, 2.0)) / cs); /* Don't use atan2 here ! */
    if (cs < 0.0)
        entree->soleil.anghoraire = entree->soleil.anghoraire + M_PI;


    /*

    Out

    */

    /* Sunrise */
    P = scol_date_sun_equation (&entree, -1.0);
    if (P < 0.0)
        P+= 24.0;
    else if (P > 24.0)
        P-= 24.0;
    h = floor (P);
    mn = floor (60.0 * (P - h));

    if (entree->soleil.solnon == 0)
        snprintf (entree->sortie.lever, 6, "%.0f:%.0f", h, mn);
    else
        snprintf (entree->sortie.lever, 6, "--");

    /* sunset */
    P = scol_date_sun_equation (&entree, 1.0);
    if (P < 0.0)
        P+= 24.0;
    else if (P > 24.0)
        P-= 24.0;
    h = floor (P);
    mn = floor (60.0 * (P - h));
    if (entree->soleil.solnon == 0)
        snprintf (entree->sortie.coucher, 6, "%.0f:%.0f", h, mn);
    else
        snprintf (entree->sortie.coucher, 6, "--");

    Mpushstrbloc (m, entree->sortie.lever);
    Mpushstrbloc (m, entree->sortie.coucher);
    MMpush (m, ITOM (2));
    MBdeftab (m);

    free (entree);

    return 0;
}




/* SCOL API */

/* functions name */
char* science_date_name[SCIENCE_DATE_PKG_NB]=
{
    "_scienceDateSun",
    "_scienceDateLeapYear",
    "_scienceDateDayNumber"
};

/* internals functions */
int (*science_date_fun[SCIENCE_DATE_PKG_NB])(mmachine m)=
{
    SCOLscienceDateSun,
    SCOLscienceDateLeapYear,
    SCOLscienceDateDayNumber
};

/* arguments(number) */
int science_date_narg[SCIENCE_DATE_PKG_NB]=
{
    3,           /* _scienceDateSun */
    1,           /* _scienceDateLeapYear */
    3           /* _scienceDateDayNumber */
};

/* functions type */
char* science_date_type[SCIENCE_DATE_PKG_NB]=
{
    "fun [[I I] [I I] [I I I]] [S S]",       /* _scienceDateSun */
    "fun [I] I",                            /* _scienceDateLeapYear */
    "fun [I I I] I"                         /* _scienceDateDayNumber */
};

int SCOLdateLoadAPI (mmachine m)
{
    int k;

    MMechostr(MSKDEBUG, "SCOLdateLoadAPI: entering\n");

    k = PKhardpak (m, "Science_Date", SCIENCE_DATE_PKG_NB, science_date_name, science_date_fun, science_date_narg, science_date_type);
    return k;
}