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

#include "scol_maths.h"

/* Resolution d'une equation du second degre : ax2 + bx + c = 0 */
static void scol_eq2nd (double *d, double a, double b, double c)
{
    double D;

    /* discriminant de l'équation */
    D = pow (b, 2.0) - (4 * a * c);
    if (D < 0)  /* solutions complexes */
    {
        d = NULL;
        return;
    }

    d[0] = (-b + sqrt (D)) / (2 * a);
    d[1] = (-b - sqrt (D)) / (2 * a);
    return;
}




/**
 * \brief Phi, nombre d'or
 *
 * Phi = (1 + 5 ^ 0.5) / 2
 *
 * \Return : F : Phi
 */
int SCOLmathsPhi (mmachine m)
{
    double phi;
    float phif;

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

    phi = (1 + pow (5, 0.5)) / 2;

    FSET (phif, (float) phi);
    MMpush (m, phif);
    return 0;
}


int SCOLmathsVolumeSphere (mmachine m)
{
    int mradius, mheight, mflag;
    int flag = SCOLMATHS_VOLUME_SPHERE;
    float radius, height, volume, out;

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

    mflag = MTOI (MMpull (m));
    mheight = MTOI (MMpull (m));
    mradius = MTOI (MMpull (m));

    if ((flag > 0) || (flag < 5))
        flag = mflag;

    switch (flag)
    {
        case (SCOLMATHS_VOLUME_SPHERE) :

            radius = FGET (mradius);
            volume = (4.0 / 3.0) * M_PI * pow (radius, 3.0);
            break;

        case (SCOLMATHS_VOLUME_SPHERICAL_CAP) :

            radius = FGET (mradius);
            height = FGET (mheight);
            volume = ((M_PI / 3.0) * pow (height, 2.0)) * ((3.0 * radius) - height);
            break;

        case (SCOLMATHS_VOLUME_SPHERE_CONE_INTERSECT) :

            height = FGET (mheight);
            volume = (M_PI / 6) * pow (height, 3.0);
            break;

        case (SCOLMATHS_VOLUME_SPHERE_CYLINDER) :

            radius = FGET (mradius);
            height = FGET (mheight);
            volume = ((2.0 * M_PI) / 3.0) * height * pow (radius, 2.0);
            break;

        default :

            volume = NIL;
    }

    FSET (out, volume);
    MMpush (m, out);
    return 0;
}



/**
 * \brief Resolution equation du second degre : ax2 + bx +c = 0
 *
 * Retourne le couple de résultats si et seulement si ce sont des valeurs relles
 * les nombres complexes n'étant pas geres.
 *
 * Stack :
 * stage 0 : F : c
 * stage 1 : F : b
 * stage 2 : F : a
 * Return : [F F]
 */
int SCOLmathsEquation2nd (mmachine m)
{
    int ma, mb, mc;
    float fx1, fx2;
    double a = 0.0, b = 0.0, c = 0.0;
    double r[2];

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

    mc =  (MMpull (m));
    mb =  (MMpull (m));
    ma =  (MMpull (m));

    if (ma != NIL)
        a = (double) FGET (ma);
    if (mb != NIL)
        b = (double) FGET (mb);
    if (mc != NIL)
        c = (double) FGET (mc);

    scol_eq2nd (r, a, b, c);

    if (r == NULL)
    {
        MMechostr (0, "SCOLmathsEquation2nd warning : no solutions in R, solutions in C\n");
        MMpush (m, NIL);
        return 0;
    }

    FSET (fx1, (float) r[0]);
    FSET (fx2, (float) r[1]);
    MMpush (m, fx1);
    MMpush (m, fx2);
    MMpush (m, ITOM (2));
    MBdeftab (m);

    return 0;
}

/**
 * \brief Cubic function : ax3 + bx2 + cx + d = 0
 *
 * Cardano method (Methode de Cardan - Tartagliani)
 *
 */
int SCOLmathsEquation3rd (mmachine m)
{
    int ma, mb, mc, md;
    float fx1, fx2, fx3;
    double a = 0.0, b = 0.0, c = 0.0, d = 0.0;
    double x1, x2, x3, p, q, D;

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

    md = MMpull (m);
    mc = MMpull (m);
    mb = MMpull (m);
    ma = MMpull (m);

    if (ma != NIL)
        a = (double) FGET (ma);
    if (mb != NIL)
        b = (double) FGET (mb);
    if (mc != NIL)
        c = (double) FGET (mc);
    if (md != NIL)
        d = (double) FGET (md);

    p = (-(pow (b, 2.0)) / (3.0 * pow (a, 2.0))) + (c / a);
    q = (b / (27.0 * a)) * (((2.0* pow (b, 2.0)) / pow (a, 2.0)) - ((9.0 * c) / a)) + d / a;

    /* discriminant */
    D = pow (q, 2.0) + ((4.0 * pow (p, 3.0)) / 27.0);
    /* TO DO */
    MMpush (m, NIL);
    return 0;
}
/* // ??? pas Cardan !!!
ÉQUATION DU TROISIÈME DEGRÉ : ax³ + bx² + cx + d = 0

SI b = c = 0 alors ax³+d = 0 alors l'algorithme général ne fonctionne pas pour ce cas particulier.

          o Les 3 solutions particulières :
            x1=(-d/a)^(1/3)
            x2= -x1/2 + i*x1*3^0.5/2
            x3= -x1/2 - i*x1*3^0.5/2
            avec i² = -1

Algorithme général de résolution :

Calculs intermédiaires :

          o a0= d/a
            a1= c/a
            a2= b/a
            a3= a2/3
            p= a1-a3*a2
            q= a0-a1*a3+2*a3^3

Calcul du discriminant : delta = (q/2)^2 + (p/3)^3

Si delta > 0 alors :

          o Calcul intermédiaire :
            w= ( -q/2 + delta^0.5 )^(1/3)
          o 3 Solutions :
            x1= w - p/3/w - a3
            x2= - (a2 + x1)/2 + i*(3^0.5 /2 * (w+p/3/w) )
            x3= - (a2 + x1)/2 - i*(3^0.5 /2 * (w+p/3/w) )
            avec i² = -1

Si delta = 0 alors :

          o 2 Solutions et 1 double :
            x1= 3*q/p - a3
            x2= -3*q/2/p - a3
            x3= x2

Si delta < 0 alors :

          o Calculs intermédiaires :
            u= 2*(-p/3)^0.5
            v= -q/2/(-p/3)^(3/2)
            t= arccos(v)/3
          o 3 Solutions :
            x1= u*cos(t)-a3
            x2= u*cos(t+2*pi/3)-a3
            x3= u*cos(t+4*pi/3)-a3
*/




/**
 * \brief Chute libre : retourne différentes valeurs selon le flag passé : fun [F F F I F] F
 *
 * Equations horaires d'un objet en chute soumis à la seule force du champ
 * de pesanteur (à l'exclusion donc de la poussée d'Archimede et des forces de
 * frottements).
 * Cette fonction permet ainsi, en fonction de paramètres initiaux tels que la
 * valeur de l'accélération de pesanteur, la vitesse initiale ou la hauteur initiale
 * de connaître une vitesse ou une altitude à un temps t, la vitesse à une altitude
 * donnée, etc.
 * Ce peut être utile en 3D pour le calcul de la trajectoire d'un objet.
 *
 * Equation horaire de la vitesse en fonction du temps :
 * V(t) = -gt + V0
 *
 * Equation de l'altitude en fonction du temps :
 * Z(t) = -1/2gt² + V0t + Z0
 *
 * où Vo et Z0 sont respectivement la vitesse et la hauteur initiales.
 *
 * Rappelons qu'une trajectoire en chute libre est une trajectoire en ligne droite,
 * éventuellement en double-sens (descente seule ou montée puis descente), dans un
 * plan vertical.
 * Rappelons également que cette fonction retourne une valeur conforme aux spécifications
 * suivantes :
 * - mouvement vertical, le vecteur unité étant dirigé vers le haut;
 * - le temps est en seconde;
 * - l'altitude en mètre;
 * - la vitesse en m.s-1 (ce sont en fait les unités internationales)
 * ce sera donc au developpeur Scol de faire les conversions si nécessaire.
 * La valeur linéaire d'accélaration est à 9.8 par défaut (valeur sur Terre).
 * La vitesse et la hauteur initiale sont à 0 par défaut.
 *
 * Selon ce qui est demandé en retour (flag), un paramètre est requis.
 *
 * Usage Scol :
 *
 * fun [F F F I F] F
 * g : F : acceleration de  pesanteur. Par défaut : 9.8 si à nil.
 * v0 : F : vitesse initiale (la valeur positive ou négative impliquera ou non une montée avant la descente)
 * z0 : F : la hauteur initiale
 * flag : I : peut valoir les valeurs suivantes :
 *          - SCOLMATHS_MVT_GETSPEED_TOTIME : avoir la vitesse en fonction du temps
 *          - SCOLMATHS_MVT_GETTIME_TOSPEED : avoir le temps écoulé en fonction de la vitesse
 *          - SCOLMATHS_MVT_GETHEIGHT_TOTIME : avoir la hauteur en fonction du temps
 *          - SCOLMATHS_MVT_GETHEIGHT_TOSPEED : avoir la hauteur en fonction de la vitesse
 *          - SCOLMATHS_MVT_GETSPEED_TOHEIGHT : avoir la vitesse en fonction de la hauteur
 *          - SCOLMATHS_MVT_GETTIME_TOHEIGHT : avoir le temps écoulé en fonction de la hauteur
 * param : F : paramètre qui correspond à la valeur demandée par le flag ("en fonction de ...")
 *          Par exemple, si flag vaut SCOLMATHS_MVT_GETSPEED_TOTIME, param vaudra le temps demandé (2.0 pour 2 secondes écoulées)
 * Retourne : F : la valeur désirée en fonction du flag
 */
int SCOLmathsChuteLibre (mmachine m)
{
    int mg, mv0, mz0, mflag, mparam;
    double G = 9.8, v0 = 0.0, z0 = 0.0, param = 0.0;
    double r;
    double *d;
    float f;

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

    mparam = MMpull (m);
    mflag = MTOI (MMpull (m));
    mz0 = MMpull (m);
    mv0 = MMpull (m);
    mg = MMpull (m);

    if (mg != NIL)
        G = FGET (mg);
    if (mv0 != NIL)
        v0 = FGET (mv0);
    if(mz0 != NIL)
        z0 = FGET (mz0);
    if(mparam != NIL)
        param = FGET (mparam);

    /* donne la vitesse au temps t */
    if (mflag == SCOLMATHS_MVT_GETSPEED_TOTIME)
        r = (-1.0) * (G * param) + v0;

    /* donne le temps pour atteindre une vitesse */
    else if (mflag == SCOLMATHS_MVT_GETTIME_TOSPEED)
        r = (v0 - param) / G;

    /* altitude à un temps t */
    else if (mflag == SCOLMATHS_MVT_GETHEIGHT_TOTIME)
        r = (-1 * 0.5 * G * pow (param, 2.0)) + v0 * param + z0;

    /* altitude atteinte pour une vitesse donnée */
    else if (mflag == SCOLMATHS_MVT_GETHEIGHT_TOSPEED)
        r = (-0.5 * (pow (v0 - param, 2.0) / G)) + (v0 * ((v0 - param) / G)) + z0;

    /* vitesse à une altitude donnée */
    else if (mflag == SCOLMATHS_MVT_GETSPEED_TOHEIGHT)
    {
        scol_eq2nd (d, (-0.5) * G, v0, z0 - param);
        r = (-1) * G * d[0] + v0;
        MMechostr (0, "SCOLmathsChuteLibre t0 : %f\n", (-1) * G * d[0] + v0);
        MMechostr (0, "SCOLmathsChuteLibre t1 : %f\n", (-1) * G * d[1] + v0);
    }
    /* tempsmis pour atteindre une altitude donnée */
    else if (mflag == SCOLMATHS_MVT_GETTIME_TOHEIGHT)
    {
        scol_eq2nd (d, (-0.5) * G, v0, z0 - param);
        r = d[0];
        MMechostr (0, "SCOLmathsChuteLibre t0 : %f\n", d[0]);
        MMechostr (0, "SCOLmathsChuteLibre t1 : %f\n", d[1]);
    }
    else
        r = NIL;

    FSET (f, r);
    MMpush (m, f);
    return 0;
}





/* SCOL API */

/* functions name */
char* science_maths_name[SCIENCE_MATHS_PKG_NB]=
{
    /* flags */
    "SCOLMATHS_VOLUME_SPHERE", "SCOLMATHS_VOLUME_SPHERICAL_CAP", "SCOLMATHS_VOLUME_SPHERE_CONE_INTERSECT",
    "SCOLMATHS_VOLUME_SPHERE_CYLINDER",

    "SCOLMATHS_MVT_GETSPEED_TOTIME", "SCOLMATHS_MVT_GETTIME_TOSPEED", "SCOLMATHS_MVT_GETHEIGHT_TOTIME",
    "SCOLMATHS_MVT_GETHEIGHT_TOSPEED", "SCOLMATHS_MVT_GETSPEED_TOHEIGHT", "SCOLMATHS_MVT_GETTIME_TOHEIGHT",

    /* functions */
    "PHIf",
    "_scienceMathsVolumeSphere",
    "_scienceMathsEquation2nd",
    "_scienceMathsChuteLibre"
};

/* internals functions */
int (*science_maths_fun[SCIENCE_MATHS_PKG_NB])(mmachine m)=
{
    /* flags */
    (bullshit) (1*2), (bullshit) (2*2), (bullshit) (3*2),
    (bullshit) (4*2),

    (bullshit) (1*2), (bullshit) (2*2), (bullshit) (3*2),
    (bullshit) (4*2), (bullshit) (5*2), (bullshit) (6*2),

    /* functions */
    SCOLmathsPhi,
    SCOLmathsVolumeSphere,
    SCOLmathsEquation2nd,
    SCOLmathsChuteLibre
};

/* arguments (number) */
int science_maths_narg[SCIENCE_MATHS_PKG_NB]=
{
    /* flags */
    TYPVAR, TYPVAR, TYPVAR,
    TYPVAR,

    TYPVAR, TYPVAR, TYPVAR,
    TYPVAR, TYPVAR, TYPVAR,

    /* functions */
    0,          /* PHIf */
    3,          /* _scienceMathsVolumeSphere */
    3,           /* _scienceMathsEquation2nd */
    5           /* _scienceMathsChuteLibre */
};

/* functions type */
char* science_maths_type[SCIENCE_MATHS_PKG_NB]=
{
    /* flags */
    "I", "I", "I",
    "I",

    "I", "I", "I",
    "I", "I", "I",

    /* functions */
    "fun [ ] F",                /* PHIf */
    "fun [F F I] F",             /* _scienceMathsVolumeSphere */
    "fun [F F F] [F F]",         /* _scienceMathsEquation2nd */
    "fun [F F F I F] F"         /* _scienceMathsChuteLibre */
};

int SCOLmathsLoadAPI (mmachine m)
{
    int k;

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

    k = PKhardpak (m, "Science_Maths", SCIENCE_MATHS_PKG_NB, science_maths_name, science_maths_fun, science_maths_narg, science_maths_type);
    return k;
}