//oscillateur anharmonique
#include <iostream>
#include <cmath>
#include "dislin.h"

using namespace std;

int main() {
  cout << "Etude de l'oscillateur anharmonique" << endl;
  cout << "Entrez la constante k/m de l'oscillateur : ";
  double k = 1.;
  cin >> k;
  cout << "Entrez l'exposant alpha (alpha = 1 -> osc. harmonique) : ";
  double alpha = 1.;
  cin >> alpha;
  cout << "Entrez l'elongation initiale : ";
  double angle0 = 10.;
  cin >> angle0;

  //initialisation de DISLIN
  metafl("XWIN");   //XWIN ou PDF
  disini();
  name("temps (s)", "X");
  name("elongation (m)", "Y"); 
  titlin("Oscillateur anharmonique", 1);
  graf(0.,100.,0.,20.,-2*angle0,2*angle0,-2*angle0,angle0);
  title();
  
  const int STEPS = 10000;
  double dt = 0.01;
  double angle[STEPS], omega[STEPS], time[STEPS];
  angle[0] = angle0;   //elongation initial 
  omega[0] = 0.;    //vitesse initiale
  time[0] = 0.;
    
  //methode de Runge
  //conditions initiales
  double acc;
  acc = - angle[0]/abs(angle[0])*k * pow(abs(angle[0]),alpha);
  omega[0] = omega[0] - acc*dt/2.;    //methode de Runge !
   
  for (int i=1; i<STEPS; i++) {
    acc = - angle[i-1]/abs(angle[i-1])*k * pow(abs(angle[i-1]),alpha);
    omega[i] = omega[i-1] + acc*dt;
    angle[i] = angle[i-1] + omega[i]*dt;
    time[i] = i*dt;
  }

  //dessin de la trajectoire
  color("RED");
  curve(time, angle, STEPS);	
  disfin();     //fin de DISLIN
  
  return 0;
}