//etude de la desintegration d'une surce radioactive
//avec la methode d'Euler
#include <iostream>
#include <cmath>
#include <iomanip>
#include "dislin.h"     //bibliotheque DISLIN

using namespace std;

int main() {  
  const int STEPMAX = 100;    //nombre de pas
  double dt;                //pas de discretisation
     
  double nuclei[STEPMAX];   //nombre d' atomes
  double temps[STEPMAX];    //intervalles temporels
  double tau;               //temps de disintegration
  
//conditions intiales   
  nuclei[0] = 100.;
  temps[0] = 0.;
  tau = 1.;

//initialisation de DISLIN
  metafl("XWIN");     //XWIN ou PDF
  disini();
  smxalf("GREEK", "{", "}", 1);
  name("temps (s)", "X");
  name("N(t)", "Y");
  titlin("Desintegration radioactive", 1);
  titlin("dN/dt = N(0)*exp(-t/{t})", 2);
  titlin("N(0) = 100, {t} = 1 s", 3);
  graf(0., 5., 0., 1., 0., 100., 0., 20.);
  title();

//solution exacte
  dt = 0.05;
  for (int step=0; step<STEPMAX-1; step++) {
    temps[step+1] = temps[step] + dt;
    nuclei[step+1] = nuclei[0] * exp(-temps[step]/tau);
  }
//dessin de la solution exacte
  thkcrv(10);
  color("BLUE");
  curve(temps, nuclei, STEPMAX);

//solution methode d'Euler (premiere iteration)
  dt = 0.5;
  for (int step=0; step<STEPMAX-1; step++) {
    nuclei[step+1] = nuclei[step] - 1./tau * nuclei[step] * dt;
    temps[step+1] = temps[step] + dt;
  }
//dessin
  color("GREEN");
  incmrk(-1);   //symboles au lieu de la curbe
  marker(21);   //symbole pour le plot
  curve(temps, nuclei, 10);


//solution methode d'Euler (deuxieme iteration)
  dt = 0.05;
  for (int step=0; step<STEPMAX-1; step++) {
    nuclei[step+1] = nuclei[step] - (1./tau * nuclei[step] * dt;
    temps[step+1] = temps[step] + dt;
  }  
//dessin
  color("RED");
  marker(0);   //symbole pour le plot
  curve(temps, nuclei, 100);

  disfin();     //fin de DISLIN

  return 0;
}