#include <iostream>
#include <cstring>
#include <cmath>
#include "dislin.h"

using namespace std;

//classe Planete
class Planete {
  public:
    //prototypes constructeurs
    Planete();     //par defaut
    Planete(char *name, double r, double t, double ecc); 
  
    //destructeur
    ~Planete();
  
    void input();    //saisie
    void calcOrbite(double t, double dt);    //calcule de l'orbite
    void position() const;     //affichage de la position de la planete
  
  private:
    void init();   //initialisation  
  
    char m_nom[20];    //nom
    double m_a;     //axe major
    double m_T;     //periode
    double m_V;     //vitesse orbitale
    double m_e;     //eccentricite
  
    double m_velIni;
    double m_posIni;
    double m_GM;     //G_N * M_Soleil

    
    double m_x0, m_y0;     //ancienne position
    double m_x, m_y;     // position actuell
    double m_t;     //temps
};

//constructeur par defaut
Planete::Planete() {
  m_nom[0] = '\0';
  m_a = 0.;
  m_T = 0.;
  m_e = 0.;
  m_velIni = 0.;
  m_posIni = 0.;
  m_GM = 0.;
  cout << "Planete construite sans initialisation !" << endl;
}    

//constructeur avec parametres   
Planete::Planete(char *name, double r, double t, double ecc=0.) {
  strcpy(m_nom,name);
  m_a = r;
  m_T = t;
  m_e = ecc;
  m_velIni = 0.;
  m_posIni = 0.;
  m_GM = 0.; 

  init();

  cout << "dans le constructeur de la planete " << m_nom << endl;
}

//destructeur
Planete::~Planete() {
  cout << "La planete " << m_nom << " a ete detruite !" << endl;
}    

//saisi de données de la planete
void Planete::input() {
  cout << "Nom de la palnete   : ";   cin >> m_nom;
  cout << "Distance (10**6 km) : ";   cin >> m_a;
  cout << "Periode (jours)     : ";   cin >> m_T;
  cout << "Eccentricite        : ";   cin >> m_e;
  init();
}

void Planete::init() {
  //conversion des distances en unites astronomique, temps en annees
  double distance = m_a/149.;
  double periode = m_T/365.26;
  //facteur GM = G_N * M_Soleil
  m_GM = pow(2.*M_PI*distance/periode,2.)*distance;
  
  //vitesse/position initiales
  //position initiale en x = a (y=0), vitesse initiale en y
  m_velIni = sqrt((m_GM/distance)*(1.+m_e)/(1.-m_e));    //perihelion
  m_posIni = distance*(1.-m_e);
  m_x0 = m_posIni;
  m_y0 = 0.;
}

//calcul de l'orbite avec la methode de Euler - Cromer
//methode equivalente a celle de Runge de deuxieme ordre
//sauf que on ne redefinit pas la velocite initiale
//par rapport a la methode d'Euler
//on calcule toujours d'abord la vitesse puis la position
void Planete::calcOrbite(double t, double dT) {
  int steps = int(t/dT);
  
  //allocation dynamique des vecteurs position, vitesse, acceleration et temps
  double *xpos = new double[steps];
  double *ypos = new double[steps];
  double *xvel = new double[steps];
  double *yvel = new double[steps];
  double *time = new double[steps];
  
  //initialisation
  xpos[0] = m_posIni;
  ypos[0] = 0.;
  xvel[0] = 0.;
  yvel[0] = m_velIni;
  
  double normR, xacc, yacc;
  //methode d'Euler - Cromer
  for (int i=1; i<steps; i++) {
    normR = sqrt(xpos[i-1]*xpos[i-1] + ypos[i-1]*ypos[i-1]);
    xacc = -m_GM*xpos[i-1] / (normR*normR*normR);
    yacc = -m_GM*ypos[i-1] / (normR*normR*normR);
    
    xvel[i] = xvel[i-1] + xacc*dT;   //vitesse 
    yvel[i] = yvel[i-1] + yacc*dT;
    
    xpos[i] = xpos[i-1] + xvel[i]*dT;   //position
    ypos[i] = ypos[i-1] + yvel[i]*dT;
  }
  m_x = xpos[steps-1];
  m_y = ypos[steps-1];
  m_t = t;
  
  //dessin de l'orbite
  curve(xpos,ypos,steps);
  //dessin de la planete
  incmrk(-1);
  marker(21);
  double x[1] = {m_x};
  double y[1] = {m_y};
  curve(x,y,1);
  incmrk(0);

  delete [] xpos, ypos, xvel, yvel;
}

void Planete::position() const {
  cout << "Planete " << m_nom << endl;
  cout << "ancienne position (x,y) :             " << m_x0 << "\t" << m_y0 << endl;
  cout << "nouvelle position (x,y) apres " << m_t << " ans : " << m_x << "\t" << m_y << endl;
} 

int main() {
  //creation des objets Planete
  Planete terre("Terre", 149.0, 365.26, 0.02);  
  Planete venus("Venus", 108.2, 224.7, 0.01);
  Planete mars("Mars", 227.94, 686.98, 0.09);
  Planete jupiter("Jupiter", 778.330, 4332.71, 0.05);
  Planete mercur("Mercur", 57.91, 87.97, 0.21);
  Planete saturne;
  saturne.input();   

  //parametres pour le calcul de l'orbite
  double time, dT;
  cout << "Entrez le temps en ans : ";
  cin >> time;
  cout << "et le pas d'integration : ";
  cin >> dT;
  
  //initialization de DISLIN
  metafl("XWIN");
  page(3000, 3000);  //pour dessiner une image proportionee
  disini();
  name("axe-X","X");
  name("axe-Y","Y");
  color("white");
  titlin("Orbites des planetes",1);
  graf(-15., 15., -15., 5., -15., 15., -15., 5.);
  title();
  
  //dessin du soleil
  color("yellow");
  incmrk(-1);
  marker(21);
  double x[1]={0.};
  double y[1]={0.};
  curve(x,y,1);
  incmrk(0);
  color("blue");
  chncrv("color");
  
  //calcul de l'orbite
  double time, dT;
  cout << "Entrez le temps en ans : ";
  cin >> time;
  cout << "et le pas d'integration : ";
  cin >> dT;
  terre.calcOrbite(time,dT);
  terre.position();
  venus.calcOrbite(time,dT);
  venus.position();
  mars.calcOrbite(time,dT);
  mars.position();
  jupiter.calcOrbite(time,dT);
  jupiter.position();
  mercure.calcOrbite(time,dT);
  mercure.position();
  saturne.calcOrbite(time,dT);
  saturne.position();
  
  disfin();   //fin de DISLIN

  return 0;
}