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

typedef struct vect{double x,y;} VECT;
double G;	// Constante des forces gravifiques
double c;	// Vitesse de la lumière
double M;	// Masse du Soleil

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 u1,u2;	// Norme de vitesse
double rv1,rv2;	// r•v
double t;	// Temps
double T;	// Intervalle de temps
int N ; 	// Nombre d'itération par maille
int A;		// (d2>d1) lors de la maille précédente, phase d'éloignement 

G=6.6738e-11;	// Constante des forces de gravité
M=1.98866e30;	// Masse du Soleil
c=299792458;	// Constante de la vitesse de la lumière

T=1;		// Intervalle de temps
N=2;		// Nombre d'itération 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("Masse du Soleil :                  %.8g\n", M);
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 :                           Modèle en relativité restreinte\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);
u1 = sqrt(v1.x*v1.x + v1.y*v1.y);
rv1 = r1.x*v1.x + r1.y*v1.y;
a1.x = ((rv1*v1.x/(c*c) - r1.x)/d1)*G*M/(d1*d1);
a1.y = ((rv1*v1.y/(c*c) - r1.y)/d1)*G*M/(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);

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

	u2 = sqrt(v2.x*v2.x + v2.y*v2.y);

        rv2 = r2.x*v2.x + r2.y*v2.y;
        a2.x = ((rv2*v2.x/(c*c) - r2.x)/d2)*G*M/(d2*d2);
        a2.y = ((rv2*v2.y/(c*c) - r2.y)/d2)*G*M/(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;
    v1=v2;
    u1=u2;
    rv1=rv2;
    a1=a2;
    t=t+T;
}
}