//etude de la trajectorie d'une particule dans
//un champ electromagnetique
#include <iostream>
#include <cmath>
#include "dislin.h"

using namespace std;

int main() {
  //parametres
  const int steps = 60000; 
  const double M = 0.001;      //masse
  const double Q = 0.001;   //charge
  const double E = 100.;    //E
  const double B = 1.;    //B
  double xpos[steps], ypos[steps], xvel[steps], yvel[steps];
  //conditions initiales
  xpos[0] = 0.0;   ypos[0] = 0.0;
  xvel[0] = 0.0;   yvel[0] = 0.0;
    
  double dT = 0.0005;   //pas d'integration

  //methode de Runge
  double xacc, yacc;     //varialbes auxiliares
  double vdrift;  
  //redefinition de conditions initiales
  xacc = Q/M*yvel[0]*B;   //acceleration
  yacc = -Q/M*xvel[0]*B + Q/M*E;
  xvel[0] = xvel[0] - xacc*dT/2.;     //vitesse 
  yvel[0] = yvel[0] - yacc*dT/2.;
  
  for (int i=1; i<steps; i++) {
    xacc = Q/M*yvel[i-1]*B;   //acceleration
    yacc = -Q/M*xvel[i-1]*B + Q/M*E;
    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;
    
    //vitesse de derive
    if (yvel[i-1]>0. && yvel[i]<0.) {
      vdrift = xpos[i] / (i*dT);
      cout << "La vitesse de derive est : " << vdrift << " m/s." << endl;
    }
  }

  //partie graphique avec DISLIN
  metafl("XWIN");   //XWIN ou PDF
  disini();     //initialisation de DISLIN

  //dessin de la trajectoire
  name("axe-X [m]", "X");
  name("axe-Y [m]", "Y");
  titlin("Champ ElectroMagnetique",1); 
  graf(0.,3100.,0.,500.,-20.,220.,-20.,40.);
  title();
  thkcrv(5);
  color("RED");
  curve(xpos,ypos,steps);
  disfin();    //fin de DISLIN
 
  return 0;
}