/* ----------------------------------------*/
/* Xavier Gnata -- Laure Gonnord, oct 2004 */
/* gccgraph -lm trois_corps.c -o trois_corps */
/* ----------------------------------------*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h> /*pour utiliser pow*/
#include "graphdsu.h"

/*fonctions f_i*/
double f (double GM0, double GM1, double GM2, double X[3], double Y[3],
   double Vx[3], double Vy[3], int i)
{
  /*écrire les fonctions ici*/
}


/* La méthode de RK ordre 4 integation equas diffs*/
void rk (double GM0, double GM1, double GM2, double X[3], double Y[3],
    double Vx[3], double Vy[3], double h)
{
  int i, i_rk;
  double TempX[3], TempY[3], TempVx[3], TempVy[3];
  double k1[12], k2[12], k3[12], k4[12];

  for (i_rk = 0; i_rk < 12; i_rk++)
    k1[i_rk] = h * f (GM0, GM1, GM2, X, Y, Vx, Vy, i_rk);

  for (i = 0; i < 3; i++)
    {
      TempX[i] = X[i] + k1[i] / 2.;
      TempY[i] = Y[i] + k1[i + 3] / 2.;
      TempVx[i] = Vx[i] + k1[i + 6] / 2.;
      TempVy[i] = Vy[i] + k1[i + 9] / 2.;
    }
  for (i_rk = 0; i_rk < 12; i_rk++)
    k2[i_rk] = h * f (GM0, GM1, GM2, TempX, TempY, TempVx, TempVy, i_rk);

  for (i = 0; i < 3; i++)
    {
      TempX[i] = X[i] + k2[i] / 2.;
      TempY[i] = Y[i] + k2[i + 3] / 2.;
      TempVx[i] = Vx[i] + k2[i + 6] / 2.;
      TempVy[i] = Vy[i] + k2[i + 9] / 2.;
    }
  for (i_rk = 0; i_rk < 12; i_rk++)
    k3[i_rk] = h * f (GM0, GM1, GM2, TempX, TempY, TempVx, TempVy, i_rk);

  for (i = 0; i < 3; i++)
    {
      TempX[i] = X[i] + k3[i];
      TempY[i] = Y[i] + k3[i + 3];
      TempVx[i] = Vx[i] + k3[i + 6];
      TempVy[i] = Vy[i] + k3[i + 9];
    }
  for (i_rk = 0; i_rk < 12; i_rk++)
    k4[i_rk] = h * f (GM0, GM1, GM2, TempX, TempY, TempVx, TempVy, i_rk);

  for (i = 0; i < 3; i++)
    {
      X[i] = X[i] + (k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i]) / 6.;
      Y[i] =
	Y[i] + (k1[i + 3] + 2 * k2[i + 3] + 2 * k3[i + 3] + k4[i + 3]) / 6.;
      Vx[i] =
	Vx[i] + (k1[i + 6] + 2 * k2[i + 6] + 2 * k3[i + 6] + k4[i + 6]) / 6.;
      Vy[i] =
	Vy[i] + (k1[i + 9] + 2 * k2[i + 9] + 2 * k3[i + 9] + k4[i + 9]) / 6.;
    }
}

int main (void)
{

  /*diverses déclarations*/

  /*conditions initiales*/

  /*demande de la précision et du nb de pas de calculs*/

  /* Affichage utilisant la lib graphdsu */

  /*Calcul et affichage, sauvegarde du dernier pas calculé*/

  /*fermeture de la fenêtre*/

  return 0;
}