#include <iostream>
#include <cmath>
#include <iomanip>

using namespace std;

int main() {   
  const double nAvogadro = 6.022141e23;
  const double secParAn = 365.*24.*3600.;
  double tau = 5.271/log(2.);     //temps de demi-vie du Co-60 (5.271 ans)
  
  cout << "***********************************************************\n";
  cout << "Programme pour le calcul de l'activite d'une source de Co60\n";
  const double A0=59.9338222e-3;
  cout << "La masse du Cobalt est A0 = " << A0 << " kg per mole.\n";
  cout << "***********************************************************\n";

  //saisi des parametres
  double M0;
  cout << "Quelle est la masse initiale de la source (en kg) ? ";
  cin >> M0;
  double t0;
  cout << "Apres combien d'ans voulez vous connaitre l'activite de la source ? ";
  cin >> t0;
  int n;
  cout << "Combien d'iterations voulez-vous faire pour ce calcul ? ";
  cin >> n;

  double S = M0/A0;   //nombre initial d'atomes / nAvogadro
  double dt = t0/n;
  cout << "valeur exacte " << exp(-t0/tau)/tau*M0/A0*nAvogadro << endl;
  cout << "valeur approx "
       << (exp((-t0+dt)/tau)-exp((-t0-dt)/tau))/(2.*dt)*M0/A0*nAvogadro
       << " (calcul incorrect !)" << endl;
  //faux parce que si on connait exp(-t0+dt)
  //on connait aussi la valeur exacte !
  for(int i=0; i< n; i++)
    S = S*(1.0-dt/tau);

  double activite = S/tau*nAvogadro; 
  cout << "Fraction remanente "<< S/M0*A0<<endl;
  cout << "L'activite de la source apres " << t0 << " ans est de " 
       << activite << " desintegrations par an" << endl; 
  cout << "equivalent a " << setprecision(10) << activite/secParAn 
       << " Bq (desintegrations par seconde)." << endl;

  return 0;
}