//******************************************************************************
//* Programme Comete de simulation d'une trajectoire de comete autour du Soleil
//* Dominique Lefebvre  novembre 2010
//* www.tangenteX.com
//******************************************************************************
#include <cstdlib>
#include <iostream>
#include <stdio.h>
#include <math.h>
#include <dislin_d.h>

using namespace std;

// Constantes physiques de l'expérience
#define PI         3.141592
#define G          6.67e-11            // constante de gravitation en USI
#define M          1.99e30             // masse du Soleil en USI
#define UA         1.496e11            // unité astronomique e, USI
 
// paramètres de simulation
#define NBPAS      2000
#define N          4          // nombre d'équations du système différentiel

//************************************************************* 
// Variables globales
//*************************************************************
double K;

//******************************************************************************
// Définition des subroutines graphiques
//******************************************************************************
void InitGraphDislin(void)
{
    metafl("XWIN");
    disini();
    pagera();
    pagfll(255);
    color("black");
    winfnt("arial");  
}

void LibelleGraphique(void)
{
    wintit("Trajectoire d'une comete");
    height(30);           // hauteur des caractères 
    hname(30);            // hauteur des labels d'axe
    htitle(40);           // hauteur du titre
//    titlin("Trajectoire d'une comete",1);
}

//******************************************************************************
// Fonctions de calcul
//******************************************************************************
double Min(double x[], int n)
{
  double min;
  int i;
  
  min = 0.0;
  for (i=0;i<n;i++)
    if (x[i] < min) min = x[i];
  return min;
}

double Max(double x[], int n)
{
  double max;
  int i;
  
  max = 0.0;
  for (i=0;i<n;i++)
    if (x[i] > max) max = x[i];
  return max;
}
//******************************************************************************
// Description du système différentiel pour RK4
//******************************************************************************
void Derivee(double y[], double dydt[])
{
    double d;
    
    d = pow((y[0]*y[0]+y[1]*y[1]),1.5);
    dydt[0] = y[2];            // dx/dt
    dydt[1] = y[3];            // dy/dt
    dydt[2] = -K*y[0]/d;       // d²x/dt² = -K*x/(x²+y²)^3/2
    dydt[3] = -K*y[1]/d;       // d²y/dt² = -K*y/(x²+y²)^3/2   
}

//******************************************************************************
// Routine RK4
//******************************************************************************
void RK4(double yin[], double yout[], double h)
{
   int i;
   double dydx[N], dydxt[N], yt[N], k1[N], k2[N], k3[N], k4[N];
  
   Derivee(yin, dydx); 
   for (i = 0; i < N; i++) 
    {
      k1[i] = h*dydx[i];
      yt[i] = yin[i] + 0.5*k1[i];
    }
   Derivee(yt, dydxt); 
   for (i = 0; i < N; i++)
    {
      k2[i] = h*dydxt[i];
      yt[i] = yin[i] + 0.5*k2[i];
    }   
   Derivee(yt, dydxt); 
   for (i = 0; i < N; i++)
    {
      k3[i] = h*dydxt[i];
      yt[i] = yin[i] + k3[i];
    }
   Derivee(yt, dydxt); 
   for (i = 0; i < N; i++)
    {
      k4[i] = h*dydxt[i];
      yout[i] = yin[i] + k1[i]/6. + k2[i]/3. + k3[i]/3. + k4[i]/6.;
    }
  return;
}

  
//******************************************************************************
// Programme principal
//******************************************************************************
int main(int argc, char *argv[])
{
    double x[NBPAS+1], y[NBPAS+1], vx[NBPAS+1], vy[NBPAS+1];
    double x0, y0, r0, v0;
    double Yin[N], Yout[N], Ep[NBPAS+1],Ec[NBPAS+1],Et[NBPAS+1];    
    double m, deltat, Periode;
    double a, b, e;
    int i, NbPas;
    char cbuf[60];

    // En-tête
    cout << "Simulation de la trajectoire d'une comete autour du Soleil\n";
    
    // définition des paramètres de la simulation
    K = 4*PI*PI;    // GM en unités astronomiques
    m = 1.;         // masse de la comète (unitaire par convention)
    
    
    // Saisie des données de l'expérience
    cout << "Aphelie (en UA): " << flush;
    cin >> r0;
    cout << "Vitesse tangentielle a l'aphelie (en UA): " << flush;
    cin >> v0;
    cout << "Periode de la comete (en UA): " << flush;
    cin >> Periode;
    
    // Paramètre de calcul
    deltat = Periode/NBPAS;    

    // calcul de la trajectoire de la comète
    x[0] = r0; y[0] =0.;
    vx[0] = 0.; vy[0] = v0;
    for(i=0;i<NBPAS;i++)
      {
        // calcul de l'énergie potentielle; cinétique et totale
        Ep[i] = -K*m/sqrt(x[i]*x[i] + y[i]*y[i]);
        Ec[i] = 0.5*m*(vx[i]*vx[i] + vy[i]*vy[i]);
        Et[i] = Ep[i] + Ec[i];
        // calcul de la trajectoire
        Yin[0] = x[i];Yin[1] = y[i];Yin[2] = vx[i];Yin[3] = vy[i];
        RK4(Yin, Yout, deltat);
        x[i+1] = Yout[0];y[i+1] = Yout[1];vx[i+1] = Yout[2];vy[i+1] = Yout[3];   
      }        

    // calcul des paramètres de la trajectoire
    a = (r0 + fabs(Min(x,NBPAS)))/2;     // demi-grand axe de l'orbite
    b = Max(y,NBPAS);                    // demi-petit axe de l'orbite
    e = sqrt(1 - (b*b)/(a*a));           // excentricité de l'orbite
    printf("a: %f  b: %f  e: %f \n",a,b,e);   
    
    // Initialisation de la page graphique
    setpag("DA4L");      // définition de la page 2970*2100
    InitGraphDislin();
    axslen(2000,900);
    LibelleGraphique();

    // tracé de la trajectoire de la comète
    axspos(500, 1000);
    axstyp("CROSS");
    graf(-(2*a-r0)-0.1,r0+0.1,-(2*a-r0)-0.1,0.5,-b-0.1,b+0.1,-b-0.1,0.5);
    messag("Trajectoire de la comete",2000,150);
    rlstrt(x[0],y[0]);
    for(i=1;i<NBPAS;i++)rlconn(x[i],y[i]);
    endgrf();
        
    // tracé du diagramme d'énergie
    axspos(500, 2000);
    axstyp("RECT");
    graf(0,Periode,0,0.5,-150,150,-150,50);
    messag("Evolution de l'energie",600,1150);
    legini(cbuf,3,20);
    // tracé de l'énergie potentielle
    rlstrt(0,Ep [0]);
    color("blue");
    legpat(0,5,-1,-1,-1,1);
    for(i=0;i<NBPAS;i++) rlconn(i*deltat,Ep[i]);
    // tracé de l'énergie cinétique
    rlstrt(0,Ec [0]);
    color("green");
    legpat(0,5,-1,-1,-1,2);
    for(i=0;i<NBPAS;i++) rlconn(i*deltat,Ec[i]);
    // tracé de l'énergie totale
    rlstrt(0,Ec [0]+Ep[0]);
    color("red");
    legpat(0,5,-1,-1,-1,3);
    for(i=0;i<NBPAS;i++) rlconn(i*deltat,Et[i]);  
    // affichage de la légende
    color("black");
    leglin(cbuf,"Energie potentielle",1);
    leglin(cbuf,"Energie cinetique",2);
    leglin(cbuf,"Energie totale",3);
    legend(cbuf,7);         
    endgrf();
    
    // Fin du programme
    disfin();
    return EXIT_SUCCESS;
}