Mouvements d'astres

Maths sur la toile

Les parties mathématiques sont créées en MathML gràce à

Il suffit d'insérer dans la page < script type="text/javascript" src="son_emplacement\ASCIIMathML.js" ></script > pour que les parties entre \$ ...\$ ou \`...\` soient interprètées comme des commandes mathématiques et présentées comme telles.

⇑Objectifs

Je souhaitais observer l'évolution de systèmes stellaires exotiques. Typiquement, des étoiles doubles avec des planètes autour. D'où plusieurs problèmes à résoudre :

⇑Mise en équations

⇑Equations de la gravitation

Pour des systèmes formés d'une étoile et d'une planète, on sait résoudre de façon exacte les équations du mouvement. Ici, avec de multiples astres, il faudra en faire une résolution approchée.
Mais tout d'abord, il faut déterminer les relations physique qui vont régir les mouvements.

la dérivée de la position est la vistesse
la dérivée de la vitesse est l'accélération qui dépend des positions relatives des astres.
La norme de la force exercée par l'astre `j` sur l'astre `i` est $\frac {G.m_i.m_j} {{A_i A_j}^2}$ dans la direction de `\vec {A_i A_j}` donc, en se ramenant à un vecteur de norme 1 :
le vecteur force est : `\frac {G.m_i.m_j} {A_i A_j^3}\vec {A_i A_j}`
La composante de l'accélération due à $A_j$ est donc : `\frac {G.m_j} {A_i A_j^3}\vec {A_i A_j}`
Le vecteur accélération de l'astre $i$ est donc `sum_{j != i}\frac {G.m_j} {A_i A_j^3} \vec {A_i A_j} `

⇑Résolution numérique des relations

Il existe différentes méthodes de résolution de ce type d'équation (dite différentielle) liant des quantitéq et leurs dérivées.

la plus simple : méthode d'Euler

la plus courante (et beaucoup plus précise) : méthode de Runge-Kutta d'ordre 4 (dite RK4)

⇑Euler

On considère ici l'équation différentielle `y'(t)=f(y(t))` ou plus généralement, `y'(t)=f(t,y(t))`.
Dans le cas des astres, la quantité `y` sera une colonne entassant position `p` et vitesse `v` : `Y(t)=((p(t)),(v(t)))` avec `Y'(t)=((v(t)),(a(t)))=f((p(t)),(v(t)))`
La variable `t` étant dans ce cas le temps (qui n'apparait pas directement dans la relation, il est caché dans les variables `p`et `v`).

la méthode d'Euler consiste à calculer une solution approchée pas `t_i` à pas `t_{i+1}=t_i+h` (longueur du pas : `h`), en approchant entre chaque pas la fonction `f` par sa valeur en `(y(t_i))
`y'` étant constant égale à `f(y(t_i))`, on a alors `y(t_i+h)=y(t_i)+h.f(y(t_i))`
de façon plus imagée, si la vitesse `v` est constante, en un temps `h` la position se sera déplacée de `h.v`.

on calcule donc de proche en proche par `y_{suivant}=y_{précédent}+h.f(y_{précédent})` ou en notation mathématiques `y_{i+1}=y_{i}+h.f(y_{i})`

⇑RK4

La méthode est ici plus subtile. On tâte le terrain (la dérivée) autour du point de départ et par une subtile moyenne pondérée, on détermine le point suivant.

`k0=f(y_i)` est la dérivée au point de départ `y_i`
`y_i+h/2*k0` est la position estimée un instant `h/2` plus tard et `k1=f(y_i+h/2*k0)` la dérivée en ce point
`k2=f(y_i+h/2*k1)` est également la dérivée estimée, mais si, cette fois, la dérivée au point de départ avait été celle estimée précédemment.
`k3= f(y_i+h*k2)` de même aprsè un temps `h` en supposant la dérivée égale à la valeur précédente.
`(k0+2k1+2k2+k3)/6` est la moyenne subtilement pondérée entre ces 4 dérivées et
`y_{i+1}=y_i+h*(k0+2k1+2k2+k3)/6` la position estimée un temps `h` plus tard.

⇑précision

Les résoltutions numériques apporchées conduisenet à deux types d'erreurs :

de méthode : la dérivée n'est pas constante égale à `f(y(t_i))`, et donc `y_{i+1}=y(t_i)+h.f(y(t_i))}` n'est pas la valeur exacte de `y(t_i+h)`
Plus le pas `h` est grand, est plus grossière sera l'approximation de la dérivée par constante.
Pour un pas donnée, l'erreur évolue proportionnellement à `1/h` et pour RK4, elle évolue comme `1/h^4`
En prenant une valeur de `h` dix fois plus petite, la précision pour Euler sera 10 fois meilleur et 10000 fois meilleur pour RK4.
de calcul : les calculs se font avec la précision de la machine, d'où des erreurs d'arrondi, en générale sur la dernière décimale.
Si les calculs se font avec 9 décimales (type float), au bout d'un million de calculs, l'erreur accumulée peut atteindre le troisième chiffre significatif.
...et cela se fait rapidement sur un processeur à 2Ghz (2 milliard d'opérations par seconde) )

Pour étudier le mouvement sur une durée donnée, on doit donc choisir entre

un pas plus petit qui diminue l'erreur de méthode mais augemente le nombre de calculs et donc les erreurs d'arrondis et

un pas plus grand qui auguemente l'erreur de méthode mais minimise les calculs (et le temps de calculs) et les erreurs d'arrondi.

⇑le code brut




#include <SDL.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>


#define HAUTEUR 800
#define LARGEUR 800
#define NB_ASTRES_MAX 10
#define PAS 0.001
#define RAFRAICHISSEMENT 20
#define G 0.5
#define NB_PAS 1


typedef struct vecteur
{
    float x,y,z;
}  vecteur;

typedef struct derivees
{
    vecteur vitesse;
    vecteur acceleration;
}  derivee;


void affecte(vecteur u,vecteur *v)
{
    v->x=u.x;
    v->y=u.y;
    v->z=u.z;

}
void diff(vecteur v1, vecteur v2, vecteur * v3)
{
    v3->x=v1.x-v2.x;
    v3->y=v1.y-v2.y;
    v3->z=v1.z-v2.z;
    //printf("difference %f %f %f:\n",v3->x,v3->y,v3->z);

}

void somme(vecteur v1, vecteur v2, vecteur * v3)
{
    v3->x=v1.x+v2.x;
    v3->y=v1.y+v2.y;
    v3->z=v1.z+v2.z;

}

// fait le produit sur place
void produit(float alpha,vecteur src,vecteur *dst)
{
    dst->x=alpha*src.x;
    dst->y=alpha*src.y;
    dst->z=alpha*src.z;
};


float prodScal(vecteur u, vecteur v)
{
    return u.x*v.x+u.y*v.y+u.z*v.z;
}

float norme(vecteur v)
{
    return sqrt(prodScal(v,v));
}

float distCarre(vecteur v1, vecteur v2)
{
    vecteur v3;
    diff( v1,v2,&v3);
    return norme(v3);

}

float dist(vecteur v1, vecteur v2)
{
    return sqrt(distCarre(v1,v2));
}

float cosalpha=sqrt(3)/2, sinalpha=0.5;

void afficheVecteur(vecteur u)
{
    printf("%f %f %f ",u.x,u.y,u.z);
}


// rotation de (u,v) d'un angle alpha
void rotation(float cosalpha,float sinalpha,vecteur *u,vecteur *v)
{
    vecteur tampon,tampon1,tampon2;
    printf("norme de u %f et de v %f avant \n",norme(*u),norme(*v));
    produit(cosalpha,*u,&tampon1);
    produit(sinalpha,*v,&tampon2);
    somme(tampon1,tampon2,&tampon);
    produit(-sinalpha,*u,&tampon1);
    produit(cosalpha,*v,&tampon2);
    somme(tampon1,tampon2,v);
    affecte(tampon,u);
    printf("norme de u %f et de v %f apres \n",norme(*u),norme(*v));

}


typedef struct astre
{
    vecteur position;
    vecteur vitesse;
    float masse;
    Uint32 couleur;
}  astre;


int initialiseLesAstres(astre * astres)
{
    FILE *fichier=NULL;
    fichier = fopen("lesastres.txt", "r");
    if (!fichier)
    {
        printf("ouverture de lesastres.txt impossible\n");
    }
    printf("ouverture de lesastres.txt  correcte\n");
    int n=NB_ASTRES_MAX;
    fscanf(fichier, "%d", &n);
    printf("nombre d'astres : %d\n", n);
    if (n>NB_ASTRES_MAX)
    {
        n=NB_ASTRES_MAX;
    }
    printf("nombre d'astres rectifie 2 : %d\n", n);

    int i=0;
    for (i=0;i<n;i++)
    {
        fscanf(fichier, "%f %f %f %f %f %f %f ",
            &astres[i].position.x  ,&astres[i].position.y,&astres[i].position.z,
            &astres[i].vitesse.x,&astres[i].vitesse.y,&astres[i].vitesse.z,
            &astres[i].masse);

    };
    fclose(fichier);

    return n;
}


// la couleur est proportionnelle à la masse du rouge au turquoise

void initialiseCouleurs(SDL_Surface *ecran,int n,astre *astres)
{
    int i=0;
    // détermine la masse maximale
    float m=0;
    for (i=0;i<n;i++)
    {
        if(astres[i].masse>m ){ m=astres[i].masse ;}
    }

    //printf("masse maximale : %f\n",m);
    // affecte les couleurs porportionellement
    Uint8 r=0,b=0,v=0;
    for (i=0;i<n;i++)
    {
        r=Uint8((1-astres[i].masse/m)*255);
        b=Uint8(astres[i].masse/m*255);
        v=b;
        printf("couleur de %d : (rouge %d vert %d bleu %d)\n",i,r,v,b);
        //Uint32 couleur=SDL_MapRGB(ecran->format,r,v,b);
        //printf("OK ?\n");

        astres[i].couleur=SDL_MapRGB(ecran->format,r,v,b);
        printf("couleur de %d : Uint32 %d)\n",i,astres[i].couleur);
    }
}

void fctBarycentre(int n, astre *astres, astre *barycentre)
{
    int i;
    float mTotale=0,m,x=0,y=0,z=0,vx=0,vy=0,vz=0;
    for (i=0;i<n;i++)
    {
        m=astres[i].masse;
        mTotale+=m;
        x+=astres[i].position.x*m;
        y+=astres[i].position.y*m;
        z+=astres[i].position.z*m;
        vx+=astres[i].vitesse.x*m;
        vy+=astres[i].vitesse.y*m;
        vz+=astres[i].vitesse.z*m;
    };
    barycentre->position.x=x/mTotale;
    barycentre->position.y=y/mTotale;
    barycentre->position.z=z/mTotale;
    barycentre->vitesse.x=vx/mTotale;
    barycentre->vitesse.z=vy/mTotale;
    barycentre->vitesse.z=vz/mTotale;
}

// projection orthogonale sur (centre,u,v)

void projection(vecteur centre,vecteur u, vecteur v, astre astreSrc, vecteur *dst)
{
    vecteur CA; // centre->astre
    diff(astreSrc.position,centre,&CA);
    // coordonnées sur u et v
    dst->x=prodScal(u,CA);
    dst->y=prodScal(v,CA);
}

// facteur d'echelle pour que tout rentre dans l'écran
float fctEchelle(int n,astre * astres, vecteur centre,vecteur u,vecteur v)
{
    float xMin=0,yMin=0,xMax=0,yMax=0,coef=1;
    int i=0;
    vecteur tampon;

    for (i=0;i<n;i++)
    {
        projection(centre,u,v,astres[i],&tampon);
        if (tampon.x<xMin) {xMin=tampon.x;}
        if (tampon.x>xMax) {xMax=tampon.x;}
        if (tampon.y<yMin) {yMin=tampon.y;}
        if (tampon.y>yMax) {yMax=tampon.y;}
    };

    // détermine le plus grand écart relatif
    coef=fabs(xMin)/LARGEUR;
    if(fabs(xMax)/LARGEUR>coef) {coef=fabs(xMax)/LARGEUR;}
    if(fabs(yMin)/HAUTEUR>coef) {coef=fabs(yMin)/HAUTEUR;}
    if(fabs(yMax)/HAUTEUR>coef) {coef=fabs(yMax)/HAUTEUR;}

    if (coef==0)
    {
        return 1;
    }
    else
    {
        return 0.5/coef;
    }
}



//colle les astres sur "ecran" après projection et mise à l'échelle
void afficheAstres(SDL_Surface *ecran,vecteur centre,vecteur u,vecteur v,float echelle, int n, astre * astres)
////////////////////
///////////////////
// il faut projeté avant de détermine rle facteur d'échelle !
{
    // noirci le ciel
    // pour chaque astre
    vecteur dst;
    int i=0;
    for (i=0;i<n;i++)
    {
        // calcule le projeté
        projection(centre,u,v,astres[i], &dst);
        // met à l'échelle
        produit(echelle,dst, &dst);

        //printf( "astre %d projete sur : %f %f \n",i,dst.x,dst.y);
        dst.x+=LARGEUR*0.5;
        dst.y+=HAUTEUR*0.5;
        SDL_Rect rect;
        rect.x=dst.x;
        rect.y=dst.y;
        rect.w=2;
        rect.h=2;
        SDL_FillRect(ecran,&rect,astres[i].couleur);

    }
    SDL_Flip(ecran);
}

// détermine le numéro de l'astre cliqué et -1 s'il n'y en a pas
int astreClique(SDL_Surface *ecran,int clicX,int clicY,int n,astre *astres,vecteur centre,vecteur u, vecteur v,float echelle)
{

    vecteur dst;
    int i=0;
    for (i=0;i<n;i++)
    {
        // calcule le projeté
        projection(centre,u,v,astres[i], &dst);
        // met à l'échelle
        produit(echelle,dst, &dst);

        dst.x+=LARGEUR*0.5;
        dst.y+=HAUTEUR*0.5;
        if (fabs(dst.x-clicX)<4 &&  fabs(dst.y-clicY)<4)
        {

            SDL_Rect rect;
            rect.x=dst.x-2;
            rect.y=dst.y-2;
            rect.w=4;
            rect.h=4;
            SDL_FillRect(ecran,&rect,SDL_MapRGB(ecran->format,255,255,255));
            SDL_Flip(ecran);
            return i;
        }
    }
    return -1;
}



int initialiseSDLVideo(SDL_Surface *ecran){
    if ( SDL_Init( SDL_INIT_VIDEO ) < 0 )
    {
        printf( "erreur d'iinitialisation video: %s\n", SDL_GetError() );
        return 1;
    }

    //s'assure de la fermeture propre
    //atexit(SDL_Quit);

    // creé la fenetre
    ecran = SDL_SetVideoMode(LARGEUR, HAUTEUR, 32,
                                           SDL_HWSURFACE|SDL_DOUBLEBUF);
    if ( !ecran )
    {
        printf("impossible de crerer la fenetre: %s\n", SDL_GetError());
        return 2;
    }
    return 0;
};


void attraction(int n,astre *astres,int numAstre,vecteur *force)
{
    force->x=0;
    force->y=0;
    force->z=0;

    float d=0,f=0;
    int i=0;
    vecteur v;
    for (i=0;i<n;i++)
    {
        if (i!=numAstre)
        {
            diff(astres[i].position,astres[numAstre].position,&v);
            d=norme(v);
            //printf("distance %d %d : %f \n",i,numAstre,d);
            if (d>0)
            {
                f=G*astres[i].masse/(d*d*d);
                produit(f,v,&v);
                somme(v,*force,force);
            }
        }
    }
}

void suivantRK4(int n,astre *astres)
{
    derivees k0[NB_ASTRES_MAX];
    derivees k1[NB_ASTRES_MAX];
    derivees k2[NB_ASTRES_MAX];
    derivees k3[NB_ASTRES_MAX];
    astre astresTampon[NB_ASTRES_MAX];

    vecteur vTampon,vSomme;

    //vecteur acceleration;
    int i=0;

    // méthode de Rounge-Kutta d'ordre 4:
    // pour la relation y'(t)=f(y(t)) et pour un pas h
    // position'=vitesse
    // vitesse'=accélération=attraction (position)
    // ici position'=vitesse et vitesse'=attraction(...)
    // on définit successivement
    // en partant de l'état y0 vers l'état y1 :
    // k0=f(y0)
    // k1=f(y0+h/2*k0)
    // k2=f(y0+h/2*k1)
    // k3= f(y0+h*k2)
    // y1=y0+h*(k0+2k1+2k2+k3)/6


    // k0=f(y0)
    for (i=0;i<n;i++)
    {
        k0[i].vitesse.x=astres[i].vitesse.x;
        k0[i].vitesse.y=astres[i].vitesse.y;
        k0[i].vitesse.z=astres[i].vitesse.z;
        attraction(n,astres,i,&(k0[i].acceleration));
    }

    // y0+h/2*k0 dans astresTampon
    for (i=0;i<n;i++)
    {
        //calcule y0+h/2*k0
        //calcule h/2*k0

        //positions
        produit(PAS/2,k0[i].vitesse,&vTampon);
        somme(astres[i].position,vTampon,&(astresTampon[i].position));

        //masse
        astresTampon[i].masse=astres[i].masse;

        //vitesses
        produit(PAS/2,k0[i].acceleration,&vTampon);
        somme(astres[i].vitesse,vTampon,&( k1[i].vitesse));
    }

    // acceleration
    for (i=0;i<n;i++)
    {
        attraction(n,astresTampon,i,&(k1[i].acceleration));
    }

    // k2=f(y0+h/2*k1)
    for (i=0;i<n;i++)
    {
        //calcule y0+h/2*k0
        //calcule h/2*k0

        //positions
        produit(PAS/2,k1[i].vitesse,&vTampon);
        somme(astres[i].position,vTampon,&(astresTampon[i].position));

        //masse
        astresTampon[i].masse=astres[i].masse;

        //vitesses
        produit(PAS/2,k1[i].acceleration,&vTampon);
        somme(astres[i].vitesse,vTampon,&( k2[i].vitesse));
    }

    // acceleration
    for (i=0;i<n;i++)
    {
        attraction(n,astresTampon,i,&(k2[i].acceleration));
    }

    // k3=f(y0+h*k2)
    //printf("k3 \n");
    for (i=0;i<n;i++)
    {
        //calcule y0+h*k0
        //calcule h*k0

        //positions
        produit(PAS,k2[i].vitesse,&vTampon);
        somme(astres[i].position,vTampon,&(astresTampon[i].position));

        //masse
        astresTampon[i].masse=astres[i].masse;

        //vitesses
        produit(PAS/2,k2[i].acceleration,&vTampon);
        somme(astres[i].vitesse,vTampon,&( k3[i].vitesse));
    }

    // acceleration
    for (i=0;i<n;i++)
    {
        attraction(n,astresTampon,i,&(k3[i].acceleration));
    }

    // y1=y0+h*(k0+2k1+2k2+k3)/6
    for (i=0;i<n;i++)
    {
        // k0+2k1+2k2+k3
        produit(2,k1[i].vitesse,&vTampon);
        somme(k0[i].vitesse,vTampon,&vSomme);
        produit(2,k2[i].vitesse,&vTampon);
        somme(vSomme,vTampon,&vSomme);
        somme(vSomme,k3[i].vitesse,&vSomme);

        // h/6*(k0+2k1+2k2+k3)
        produit(PAS/6,vSomme,&vSomme);

        //y1=y0+h/6*(k0+2k1+2k2+k3)
        somme(astres[i].position,vSomme,&(astres[i].position));


        produit(2,k1[i].acceleration,&vTampon);
        somme(k0[i].acceleration,vTampon,&vSomme);
        produit(2,k2[i].acceleration,&vTampon);
        somme(vSomme,vTampon,&vSomme);
        somme(vSomme,k3[i].acceleration,&vSomme);
        produit(PAS/6,vSomme,&vSomme);
        somme(astres[i].vitesse,vSomme,&(astres[i].vitesse));


    }

}

void suivantEuler(int n,astre *astres)
{
    astre astresSuivants[NB_ASTRES_MAX];

    vecteur acceleration,tampon;
    int i=0;

    for (i=0;i<n;i++)
    {
        produit(PAS,astres[i].vitesse,&tampon);
        somme(tampon,astres[i].position,&(astresSuivants[i].position));
        attraction(n,astres,i,&acceleration);
        produit(PAS,acceleration,&acceleration);
        somme(astres[i].vitesse,acceleration,&(astresSuivants[i].vitesse));

    }
    for (i=0;i<n;i++)
    {
        astres[i].position.x=astresSuivants[i].position.x;
        astres[i].position.y=astresSuivants[i].position.y;
        astres[i].position.z=astresSuivants[i].position.z;
        astres[i].vitesse.x=astresSuivants[i].vitesse.x;
        astres[i].vitesse.y=astresSuivants[i].vitesse.y;
        astres[i].vitesse.z=astresSuivants[i].vitesse.z;
    }

}


    // calcule la nouvelle position suivant le pas de la méthode RK4
    // avec nbPas par RAFRAICHISSEMENT (1/20ème de seconde )
    // signale s'il n'y arrive pas
    // et modifie l'affichage

    // le centre est repéré par un numéro  : -1 pour le barycenter, et le numéro de l'astre sinon.

void evolution(SDL_Surface *ecran,int n,astre *astres, int *numCentre,vecteur *u, vecteur *v,vecteur *w,float *echelle,int *nbPas)
{
    vecteur centre;
    astre barycentre;
    bool trace=true;

    if (*numCentre==-1)
    {
        fctBarycentre(n,astres,&barycentre);
        affecte(barycentre.position,¢re);
    }
    else
    {
        affecte(astres[*numCentre].position,¢re);
    }
    printf("centre position.x %f %f %f\n",centre.x,centre.y,centre.z);

    int nb=0,i=0; // nombre de pas effectués

    SDL_Event event; bool poursuit=true,pause=false,efface=false;
    Uint32 temps=SDL_GetTicks();
    int numAstre=0;
    int clicX=0, clicY=0;

    // avant chaque calcul de positoin, vérifie les événements et
    // si le temps est passé, le signale.
    while (poursuit)
        {
            // gere les interractions clavier
            if( SDL_PollEvent(&event))
            {
                switch(event.type)
                {
                    case SDL_QUIT : // Si c'est un évènement de type "Quitter"
                        poursuit=false;
                    break;
                    case SDL_MOUSEBUTTONDOWN:
                        clicX=event.button.x, clicY=event.button.y;
                        numAstre=astreClique(ecran,clicX,clicY,n,astres, centre,*u,*v,*echelle);

                        // changement de centre
                        if (event.button.button==SDL_BUTTON_LEFT)
                        {
                            *numCentre=numAstre;
                            SDL_FillRect(ecran,NULL,SDL_MapRGB(ecran->format, 0,0,0));
                        }

                        // efface l'étoile sélectionnée (couleur noir, poids nul)
                        if (event.button.button==SDL_BUTTON_RIGHT)
                        {
                            efface=false;
                            if (numAstre>-1)
                            {
                                astres[numAstre].masse=0;
                                astres[numAstre].couleur=0;
                            }
                        }
                    break;

                    case SDL_KEYDOWN : // une touche clavier est enfoncée
                        switch(event.key.keysym.sym)
                        {
                            case SDLK_q: // a en clavier azerty : ajuste et recentre la fenetre

                                *echelle=fctEchelle(n,astres,centre,*u,*v)*0.9;
                                for (i=0;i<n;i++)
                                {
                                    diff(astres[i].position,centre,&(astres[i].position));
                                }
                                printf("recentrage sur le centre\n");
                                SDL_FillRect(ecran,NULL,SDL_MapRGB(ecran->format, 0,0,0));

                            break;

                            case SDLK_t: // active/désactive la trace
                                trace=!trace;
                            break;

                            case SDLK_KP_PLUS: // zoom avant
                                SDL_FillRect(ecran,NULL,SDL_MapRGB(ecran->format, 0,0,0));
                                *echelle=*echelle*2;
                                printf("zoom avant\n");
                            break;


                            case SDLK_KP_MINUS: // zoom arrière
                                //SDL_FillRect(ecran,NULL,SDL_MapRGB(ecran->format, 0,0,0));
                                *echelle=*echelle/=2;
                                printf("zoom arrière \n");
                            break;

                             case SDLK_DOWN: // rotation autour de u
                                //SDL_FillRect(ecran,NULL,SDL_MapRGB(ecran->format, 0,0,0));
                                rotation(cosalpha,sinalpha,v,w);
                            break;
                            case SDLK_KP2: // rotation autour de u
                                //SDL_FillRect(ecran,NULL,SDL_MapRGB(ecran->format, 0,0,0));
                                rotation(cosalpha,sinalpha,v,w);
                            break;


                            case SDLK_UP : // rotation autour de u
                                //SDL_FillRect(ecran,NULL,SDL_MapRGB(ecran->format, 0,0,0));
                                rotation(cosalpha,-sinalpha,v,w);
                            break;
                            case  SDLK_KP8: // rotation autour de u
                                //SDL_FillRect(ecran,NULL,SDL_MapRGB(ecran->format, 0,0,0));
                                rotation(cosalpha,-sinalpha,v,w);
                            break;

                            case SDLK_RIGHT: // rotation autour de v
                                //SDL_FillRect(ecran,NULL,SDL_MapRGB(ecran->format, 0,0,0));
                                rotation(cosalpha,sinalpha,u,w);
                            break;
                            case SDLK_KP6: // rotation autour de v
                                //SDL_FillRect(ecran,NULL,SDL_MapRGB(ecran->format, 0,0,0));
                                rotation(cosalpha,sinalpha,u,w);
                            break;

                            case SDLK_LEFT: // rotation autour de v
                                //SDL_FillRect(ecran,NULL,SDL_MapRGB(ecran->format, 0,0,0));
                                rotation(cosalpha,-sinalpha,u,w);
                            break;
                            case SDLK_KP4: // rotation autour de v
                                //SDL_FillRect(ecran,NULL,SDL_MapRGB(ecran->format, 0,0,0));
                                rotation(cosalpha,-sinalpha,u,w);
                            break;



                            case SDLK_p: // pause attend le retour d'un p
                                pause=!pause;
                            break;

                            case SDLK_r: // plus rapide
                                *nbPas*=2;
                            break;


                            case SDLK_l: // plus lent
                                *nbPas/=2; *nbPas+=1;
                            break;


                            case SDLK_ESCAPE: // recentre et ajuste la fenetre (clavier qwerty)
                                poursuit=false;
                            break;

                            default:
                            break;
                                printf("autre touche\n");
                            break;

                        }
                    break;
                }
            }
            if (!poursuit)
            {
                break;
            }
            // vérifie que le nombre depas n'est pas excessif
            if (!pause)
            {
                if (SDL_GetTicks()>temps+RAFRAICHISSEMENT)
                {
                    printf("nombre de pas par rafraichissement réalisé %d pour %d\n",nb,*nbPas);
                    temps=SDL_GetTicks();
                }
                // si le nombre de pas est effectué, attend le temps d'affichage
                if (nb>*nbPas)
                {
                    nb=0;
                    //printf("attend la fin du dé&lai avant affichage\n");
                    if (*numCentre==-1)
                    {
                        fctBarycentre(n,astres,&barycentre);
                        affecte(barycentre.position,¢re);
                    }
                    else
                    {
                        affecte(astres[*numCentre].position,¢re);
                    }
                    if (trace)
                    {
                        SDL_FillRect(ecran,NULL,SDL_MapRGB(ecran->format, 0,0,0));
                    }
                    afficheAstres(ecran,centre,*u,*v,*echelle,n,astres);
                    while(SDL_GetTicks()<temps+20){};
                    temps=SDL_GetTicks();
                }
                else
                {
                    suivantRK4(n,astres);
                    nb++;
                }
            }
        }
}





int main ( int argc, char** argv )
{
    // initialisation de l'écran
    SDL_Surface *ecran;
    ecran = SDL_SetVideoMode(LARGEUR, HAUTEUR, 32,
                SDL_HWSURFACE|SDL_DOUBLEBUF);
    if ( !ecran )
    {
        printf("impossible de crerer la fenetre: %s\n", SDL_GetError());
        return 2;
    }
    atexit(SDL_Quit);
    printf("ecran cree \n");


    astre astres[NB_ASTRES_MAX];
    printf("creation des astres tres correcte\n");


    int n=NB_ASTRES_MAX;
    n=initialiseLesAstres(astres);
    printf("initialisation des astres correcte\n");


    initialiseCouleurs(ecran, n,astres);
    printf("initialisation des couleurs correcte\n");


    astre barycentre;
    fctBarycentre(n,astres, &barycentre);

    float coef=1;
    vecteur u={1,0,0};
    vecteur v={0,1,0};
    vecteur w={0,0,1};

    coef=fctEchelle(n,astres, barycentre.position,u,v)*0.9;
    printf("coef %f \n",coef);


    afficheAstres(ecran,barycentre.position,u,v,coef,n,astres);


    int nbPas=NB_PAS;int numCentre=-1;
    evolution(ecran,n,astres, &numCentre,&u,&v,&w,&coef,&nbPas);


    SDL_FreeSurface(ecran);

    // all is well ;)
    printf("Sortie correcte\n");
    return 0;

}