//etude des orbites de deux etoiles binaires
//methode de Runge
#include <iostream>
#include <cmath>
#include "dislin.h"

/*
fonction SLEEP(ms)
malheureusement la fonction SLEEP est differente sur Windows et Unix :-(
pour pouvoir utiliser le meme programme sur Windows et Unix,
on utilise le directives de precompilation pour choisir la bonne fonction
*/  
#ifdef __UNIX__
#include <unistd.h>
#define SLEEP(ms) usleep(ms*1000);
#else
#include <windows.h>
#define SLEEP(ms) Sleep(ms);
#endif

using namespace std;

int main() {
  const int steps = 30000;
  const double GN = 1.;
  //attention aux conditions initiales  
  //etoile A
  const double MA = 1.;  
  double xposA[steps], yposA[steps], xvelA[steps], yvelA[steps]; 
  xposA[0] = 0.3; yposA[0] = 0.;
  xvelA[0] = 0.; yvelA[0] = 1.0;
  //etoile B
  const double MB = MA;
  double xposB[steps], yposB[steps], xvelB[steps], yvelB[steps];   
  xposB[0] = -xposA[0]; yposB[0] = 0.;
  xvelB[0] = 0.; yvelB[0] = -yvelA[0];
  
  double dT = 0.005;
  
  //methode de Runge
  double dist, xaccA, yaccA, xaccB, yaccB;   //variables auxiliares
  //redefinition des vitesses initiales
  dist = sqrt(pow((xposA[0]-xposB[0]),2) + pow((yposA[0]-yposB[0]),2)); 
  //etoile A
  xaccA = -GN*MB / pow(dist,3) * (xposA[0]-xposB[0]);
  yaccA = -GN*MB / pow(dist,3) * (yposA[0]-yposB[0]);    
  xvelA[0] = xvelA[0] -xaccA*dT/2.;   //vitesse 
  yvelA[0] = yvelA[0] -yaccA*dT/2.;
  //etoile B  
  xaccB = -GN*MA / pow(dist,3) * (xposB[0]-xposA[0]);
  yaccB = -GN*MA / pow(dist,3) * (yposB[0]-yposA[0]);
  xvelB[0] = xvelB[0] -xaccB*dT/2.;   //vitesse 
  yvelB[0] = yvelB[0] -yaccB*dT/2.;
  for (int i=1; i<steps; i++) {
    //distance entre les etoiles      
    dist = 
      sqrt(pow((xposA[i-1]-xposB[i-1]),2)+pow((yposA[i-1]-yposB[i-1]),2)); 
    
    //etoile A
    xaccA = -GN*MB / pow(dist,3) * (xposA[i-1]-xposB[i-1]);
    yaccA = -GN*MB / pow(dist,3) * (yposA[i-1]-yposB[i-1]);
    
    xvelA[i] = xvelA[i-1] + xaccA*dT;   //vitesse 
    yvelA[i] = yvelA[i-1] + yaccA*dT;
    
    xposA[i] = xposA[i-1] + xvelA[i]*dT;   //position
    yposA[i] = yposA[i-1] + yvelA[i]*dT;
    
    //etoile B
    //A et B sont inverse dans le calcul de l'acceleration !!!
    xaccB = -GN*MA / pow(dist,3) * (xposB[i-1]-xposA[i-1]);
    yaccB = -GN*MA / pow(dist,3) * (yposB[i-1]-yposA[i-1]);
    
    xvelB[i] = xvelB[i-1] + xaccB*dT;   //vitesse 
    yvelB[i] = yvelB[i-1] + yaccB*dT;
    
    xposB[i] = xposB[i-1] + xvelB[i]*dT;   //position
    yposB[i] = yposB[i-1] + yvelB[i]*dT;
  }
    
  //graphique DISLIN
  metafl("XWIN");
  page(3000, 3000);
  disini();
  name("axe-X","X");
  name("axe-Y","Y");
  titlin("Animation d'un systeme binaire",1);

/*
  //dessin des orbites
  graf(-1.,1.,-1.,1.,-1.,1.,-1.,1.);
  thkcrv(2);
  color("red"); 
  curve(xposA, yposB, steps);
  color("yellow"); 
  curve(xposB, yposB, steps);
*/
  
  //animation DISLIN
  double xx[1], yy[1];   //variables auxilieres
  for (int j=1 ; j<steps; j++) {
    //pour accelerer l'animation on ne dessine que toutes les 100 iterations
    SLEEP(5);    //pour arreter l'execution pendant 10 ms
    if (j%10==0) {
      erase();   //on efface le tableau
      color("FORE");   //on remet un axe decrit en dessous en blanc
      graf(-0.5,0.5,-0.5,0.2,-0.5,0.5,-0.5,0.2);
      title();

      //dessiner les etoiles
      incmrk(-1);
      color("GREEN");
      marker(21);
      xx[0]=xposA[j];
      yy[0]=yposA[j];
      curve(xx,yy,1);
      
      incmrk(-1);
      color("RED");
      marker(21);
      xx[0]=xposB[j];
      yy[0]=yposB[j];
      curve(xx,yy,1);
      
      sendbf();
      endgrf();
    }
  }
  
  disfin();   //fin de DISLIN

  return 0;
}