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

using namespace std;

int main () {
  const float dt = 0.1;
  const int nPasMax= int(1000/dt) + 1;
	
  double A[nPasMax];
  double B[nPasMax];
  double C[nPasMax];
  double temps[nPasMax];

  //conditiones initiales	
  A[0]=10000;
  B[0]=10000;
  C[0]=0;
  const double tauA = 100.;
  const double tauB = 200.;
  temps[0]=0.;

  double nucleimA;
  double nucleimB;
  for (int i=0; i<nPasMax-1; i++) {
  	//pour calculer N(t_i + Dt/2) on utilse la methode d'Euler
    nucleimA = A[i] - 1./tauA*A[i]*dt/2.;
    A[i+1] = A[i] - 1./tauA *nucleimA*dt;
    nucleimB = B[i] - 1./tauB*B[i]*dt/2.;
    B[i+1] = B[i] - 1./tauB *nucleimB*dt;
    //on calcule C a partir de nayaux B
	//qui se sont desintegres entre t et t + Delta t
    C[i+1] = C[i] +  1./tauB*B[i]*dt;
    //on ajout a B le noyaoux A
	//qui se sont desintegres entre t et t + Delta t
    B[i+1] = B[i+1] + 1./tauA*A[i]*dt;

    temps[i+1] = temps[i]+dt;
  }

  //initialisation DISLIN	
  metafl("XWIN"); 
  disini();
  name("temps (s)", "X");
  name("N(t)", "Y");
  graf(0., 1000., 0., 200., 0., 20000., 0., 2000.);
  title();
  //dessin noyaux A
  color("RED");
  marker(1); 
  curve(temps, A, nPasMax-1);
  //dessin noyaux B  	
  color("BLUE");
  marker(1); 
  curve(temps, B, nPasMax-1);
  //dessin noyaux C  	
  color("GREEN");
  marker(1); 
  curve(temps, C, nPasMax-1);
  //fin de DISLIN	
  disfin();
	
  return 0;
}