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

typedef struct vect{double x,y;} VECT;
double G;	// Constante des forces gravifiques
double M;	// Masse du Soleil
double r0;	// Rayon du Soleil
double c;       // Vitesse de la lumière
double K;       // Constante servant au calcul de la masse du champ


double Aphelie(double t, VECT r1,VECT v1,VECT a1,VECT p){
    VECT r,v,a;
    double rv,vv,ra;	
    double d,h;
    int i;
    static double ta1,ta2; 	// Aphélie
    static double aa1,aa2; 	// Aphélie
    h=0;
    for(i=0;i<6;i++){
        r.x = r1.x + v1.x*h + a1.x*h*h/2 + p.x*h*h*h/6;
        r.y = r1.y + v1.y*h + a1.y*h*h/2 + p.y*h*h*h/6;
        v.x = v1.x + a1.x*h + p.x*h*h/2;
        v.y = v1.y + a1.y*h + p.y*h*h/2;
        d = sqrt(r.x*r.x + r.y*r.y); 
        a.x = - G*M*(r.x / d) / (d * d);
        a.y = - G*M*(r.y / d) / (d * d);
        rv = r.x*v.x + r.y*v.y;
        vv = v.x*v.x + v.y*v.y;
        ra = r.x*a.x + r.y*a.y;
        h = h - rv/(vv+ra);
    }
    ta2 = t+h;
    aa2 = acos(r.x/d)*(180/3.14159265358979324);
    printf("%6.3f µ°/Ms    %10.7f Ms\n",((aa2-aa1)*1e6 )/( (ta2-ta1)*1e-6),(ta2-ta1)*1e-6);
    ta1 = ta2;
    aa1 = aa2;    
    return h ;
}

int main() {  
VECT r1,r2;	// Position
VECT v1,v2;	// Vitesse
VECT a1,a2;	// Acceleration
VECT p;		// Pente
double d1,d2;	// Distance
double M1,M2;	// Masse du Soleil et de son champ de gravité
double t;	// Temps
double T;	// Intervalle de temps
int N ; 	// Nombre d'itération par maille
int A;		// A=(d2>d1) lors de la maille précédente

G=6.6738e-11; 	// Constante des forces de gravité
M=1.98865e30;	// Masse du Soleil
r0 = 6.96342e8; // Rayon du Soleil
c = 2.99792458e8; // Vitesse de la lumière
K = G*M*M/(2*c*c);// Constante servant au calcul de la masse du champ

T=1;		// Intervalle de temps
N=2;		// Nombre d'itérations par maille

r1.x=6.9817079e10;
r1.y=0;
v1.x=0;
v1.y=3.886e4;

printf("------------------------------------------------------------------------------\n");
printf("Constante des forces de gravité :  %.8g\n", G);
printf("Vitesse de la lumière :            %.8g\n", c);
printf("Masse du Soleil :                  %.8g\n", M);
printf("rayon du Soleil :                  %.8g\n", r0);
printf("Constante de masse de champ :      %.8g\n", K);
printf("Position de Mercure :              (%.8g, %.8g)\n", r1.x, r1.y);
printf("Vitesse de Mercure :               (%.8g, %.8g)\n", v1.x, v1.y);
printf("Modèle :                           Newton + masse du champ\n");
printf("Intervalle de temps d'une maille : %.8g\n", T);
printf("Nombre d'itération par maille :    %d\n", N);
printf("------------------------------------------------------------------------------\n");
printf(" Précession      Période\n");

t=0;
A=0;

d1 = sqrt(r1.x*r1.x + r1.y*r1.y);
M1 = M+(1/r0-1/d1)*K;
a1.x = - G*M1*r1.x/(d1*d1*d1);
a1.y = - G*M1*r1.y/(d1*d1*d1);
p.x = 0;
p.y = 0;

while(1){
    int i;
    for(i=0;i<N;i++){
        r2.x = r1.x + (v1.x + (a1.x/2 + p.x*T/6)*T)*T;
        r2.y = r1.y + (v1.y + (a1.y/2 + p.y*T/6)*T)*T;
	d2 = sqrt(r2.x*r2.x + r2.y*r2.y);
        M2 = M+(1/r0-1/d2)*K;

	v2.x = v1.x + (a1.x + p.x*T/2)*T;
	v2.y = v1.y + (a1.y + p.y*T/2)*T;

        a2.x = - G*M2*r2.x/(d2*d2*d2);
        a2.y = - G*M2*r2.y/(d2*d2*d2);

        p.x = (a2.x - a1.x)/T;
        p.y = (a2.y - a1.y)/T;
    }

    if (d2>d1 && A==0) A=1; 
    if (d2<d1 && A==1) {A=0; Aphelie(t,r1,v1,a1,p);}

    r1=r2;
    d1=d2;
    M1=M2;
    v1=v2;
    a1=a2;

    t=t+T;
}
}