//******************************************************************************
//* Programme Rutherford1 de simulation de l'expérience de Rutherford
//* Dominique Lefebvre  novembre 2010
//* www.tangenteX.com
//******************************************************************************
#include <cstdlib>
#include <iostream>
#include <fstream>
#include <stdio.h>
#include <math.h>
#include <dislin_d.h>

using namespace std;

// Paramètres physiques de l'expérience
#define PI         3.141592
#define epsilon0   8.85419e-12
#define e          1.60218e-19
#define Zalpha     2.                  // nb de charges du noyau He++
#define Zau        79.                 // nb de charges du noyau Au
#define A          197                 // nombre de masse de l'or
#define Malpha     6.696e-27           // masse du noyau He++
#define r0         1.27e-15            // rayon du noyau d'hydrogène
#define c          3e+8                // vitesse de la lumière

// Zone d'affichage
#define XMIN       -600
#define XMAX       500
#define YMIN       -500
#define YMAX       500   

// Paramètres du système de détecteurs
#define RD          400      // rayon du détecteur en u.c.
#define OUV         10       // ouverture angulaire d'un détecteur en °
#define NBD         36       // Nombre de détecteurs

// Paramètres de la fonction ran3
#define MBIG 1000000000
#define MSEED 161803398
#define MZ 0
#define FAC (1.0/MBIG)

#define N          4          // nombre d'équations du système différentiel

#define Distance(x,y)  (sqrt((x)*(x) + (y)*(y)))

//************************************************************* 
// Variables globales
//*************************************************************
int Comptage[NBD];
double K;
bool XDetecteur;

//******************************************************************************
// Définition des subroutines graphiques
//******************************************************************************
void InitGraphDislin(void)
{
    metafl("XWIN");
    disini();
    pagera();
    pagfll(255);
    color("black");
    winfnt("arial");  
}

void LibelleGraphique(void)
{
    wintit("Experience de Rutherford");
    height(30);           // hauteur des caractères 
    hname(30);            // hauteur des labels d'axe
    htitle(40);           // hauteur du titre
    titlin("Experience de Rutherford",1);
}

//******************************************************************************
// Routine de calcul du germe pour rand(), à partir du nombre de ticks de la
// CPU depuis son démarrage - ATTENTION: valable seulement sur archi x86
//******************************************************************************
int rdtsc()
{
    __asm__ __volatile__("rdtsc");
}


//******************************************************************************
// Routine de production d'un nombre aléatoire distribué entre 0.0 et 1.0
// D'après Numerical Recipes in C §7.1
//******************************************************************************
float ran3(int *idum)
{
	static int inext,inextp;
	static long ma[56];
	static int iff=0;
	long mj,mk;
	int i,ii,k;

	if (*idum < 0 || iff == 0) {
		iff=1;
		mj=MSEED-(*idum < 0 ? -*idum : *idum);
		mj %= MBIG;
		ma[55]=mj;
		mk=1;
		for (i=1;i<=54;i++) {
			ii=(21*i) % 55;
			ma[ii]=mk;
			mk=mj-mk;
			if (mk < MZ) mk += MBIG;
			mj=ma[ii];
		}
		for (k=1;k<=4;k++)
			for (i=1;i<=55;i++) {
				ma[i] -= ma[1+(i+30) % 55];
				if (ma[i] < MZ) ma[i] += MBIG;
			}
		inext=0;
		inextp=31;
		*idum=1;
	}
	if (++inext == 56) inext=1;
	if (++inextp == 56) inextp=1;
	mj=ma[inext]-ma[inextp];
	if (mj < MZ) mj += MBIG;
	ma[inext]=mj;

	return mj*FAC;
}

//******************************************************************************
// Calcul du critère d'arrêt de la boucle de calcul
//******************************************************************************
bool Arret(double x, double y)
{
    bool Fin=false;
    
    // la particule entre dans la chambre de comptage
    if (x >= -5) XDetecteur = true;
    // la particule heurte le détecteur - Arrêt du calcul
    if (XDetecteur)
      if (Distance(x,y) >= RD) Fin = true;
    
    return Fin;    
}

//******************************************************************************
// Description du système différentiel pour RK4
//******************************************************************************
void Derivee(double y[], double dydt[])
{
    double d;
    
    d = pow((y[0]*y[0]+y[1]*y[1]),1.5);
    dydt[0] = y[2];            // dx/dt
    dydt[1] = y[3];            // dy/dt
    dydt[2] = K*y[0]/d;        // d²x/dt²
    dydt[3] = K*y[1]/d;        // d²y/dt²    
}

//******************************************************************************
// Routine RK4
//******************************************************************************
void RK4(double yin[], double yout[], double h)
{
   int i;
   double dydx[N], dydxt[N], yt[N], k1[N], k2[N], k3[N], k4[N];
  
   Derivee(yin, dydx); 
   for (i = 0; i < N; i++) 
    {
      k1[i] = h*dydx[i];
      yt[i] = yin[i] + 0.5*k1[i];
    }
   Derivee(yt, dydxt); 
   for (i = 0; i < N; i++)
    {
      k2[i] = h*dydxt[i];
      yt[i] = yin[i] + 0.5*k2[i];
    }   
   Derivee(yt, dydxt); 
   for (i = 0; i < N; i++)
    {
      k3[i] = h*dydxt[i];
      yt[i] = yin[i] + k3[i];
    }
   Derivee(yt, dydxt); 
   for (i = 0; i < N; i++)
    {
      k4[i] = h*dydxt[i];
      yout[i] = yin[i] + k1[i]/6. + k2[i]/3. + k3[i]/3. + k4[i]/6.;
    }
  return;
}

//******************************************************************************
// Initialisation du système de détecteurs
//******************************************************************************
void InitialisationDetecteur(void)
{
  double theta;
  int i;
  
  // Tracé du système de détecteurs
  color("green");    // couleur du système de détecteurs
  theta = 0.0;
  rlstrt(RD*cos(theta),RD*sin(theta));
  for(theta=0.0;theta<=2*PI;theta += OUV*PI/180) 
     rlconn(RD*cos(theta),RD*sin(theta));  
  // Initialisation du système de comptage
  for(i=0;i<NBD;i++) Comptage[i] = 0;
  
  return;   
}

//******************************************************************************
// Sélectionne le compteur à partir de la valeur de l'angle en degrés
//******************************************************************************
int SelectionneCompteur(int angle)
{
  int i,j,compteur;
  
  compteur = 0;
  if ((angle == 180) || (angle == -180))
    compteur = 35;
  else if (angle >= 0)
    {
      for(i=0,j=0;i<18;i++,j+=OUV)
          if ((angle >= j)&&(angle <= (j + 9))) 
            {
              compteur = i;
              break;
            }
    }
  else   // angle négatif
    {
      for(i=0,j=0;i<18;i++,j+=OUV)
          if ((angle < -j)&&(angle >= -(j + 9))) 
            {
              compteur = 18 + i;
              break;
            }
    }    

  return compteur;    
}
//******************************************************************************
// Routine de comptage des particules selon l'ange de diffusion
//******************************************************************************
void ComptageParticules(double vx, double vy)
{
  int theta;
  
  // Calcul de l'angle de diffusion en degré
  theta = (int)floor(atan2(vy,vx)*180/PI);
  // Cumul dans le détecteur correspondant
  Comptage[SelectionneCompteur(theta)]++;
 
  return;
}
  
//******************************************************************************
// Programme principal
//******************************************************************************
int main(int argc, char *argv[])
{
    double lc, tc, k, kprime;
    double deltat,RFaisceau;
    double x, y;
    double vx, vy;
    double x0, y0, v0,E0;
    double Yin[N], Yout[N];
    double b;
    int i,idum,NbParticules,NbPartDetectee;

    // En-tête
    cout << "Simulation de l'experience de Rutherford\n";
    
    // paramètres physiques
    k = Zalpha*Zau*e*e/(4.*PI*epsilon0);
    kprime = k/Malpha;
    lc = r0*pow(A,0.33333);    // longueur caractéristique = rayon du noyau
    tc = lc/c;                 // temps caractéristique
    K = kprime/(pow(c,2)*lc);  // coefficient adimensionné  

    // définition des paramètres de la simulation
    deltat = 0.01;    // unité caractéristique

    // Saisie des données de l'expérience
    cout << "Nombre de particules dans le faisceau: " << flush;
    cin >> NbParticules;
    cout << "Rayon du faisceau de particules (en unite caracteristique): " << flush;
    cin >> RFaisceau;
    cout << "Energie cinetique des particules alpha (en MeV): " << flush;
    cin >> E0;
    v0 = sqrt(2*E0*1.6e-13/Malpha)/c;   // vitesse initiale en unité caractéristique    
    
    // Initialisation de la page graphique
    InitGraphDislin();
    LibelleGraphique();

    // Définition de la zone de tracé graduée en unités caractéristiques
    graf(XMIN,XMAX,XMIN,200,YMIN,YMAX,YMIN,100);
    title();
    
    // tracé et initialisation du système de détecteurs
    color("red");rlcirc(0.0,0.0,1);           // tracé du noyau d'or cible
    InitialisationDetecteur();                             

    // Tir des particules alpha du faisceau
    printf("Experience en cours ...\n");
    color("blue");     // couleur du faisceau
    for (i=0; i<NbParticules;i++)    
      {
        // Détermination du paramètre d'impact
        idum = -rdtsc();   // initialisation du générateur aléatoire
        b = -RFaisceau + 2*RFaisceau*ran3(&idum);
        // initialisation des données
        x0 = XMIN;
        x = x0; y = b;
        vx = v0; vy = 0.0;
        rlstrt(x0,b); 
        // boucle de calcul et tracé
        XDetecteur = false;
        while (!Arret(x,y))
          {
            Yin[0] = x;Yin[1] = y;Yin[2] = vx;Yin[3] = vy;
            RK4(Yin, Yout, deltat);
            x = Yout[0];y = Yout[1];vx = Yout[2];vy = Yout[3];
            rlconn(x,y);  // tracé du point     
          }
        // Comptage dans le détecteur correspondant
        ComptageParticules(vx,vy);
      }
     
    // Comptage du nombre de particules détectées
    NbPartDetectee = 0;
    for(i=0;i<NBD;i++) NbPartDetectee += Comptage[i];
    
    //Enregistrement des compteurs
    ofstream SauvCompteurs("RutherfordCompteurs.txt");
    SauvCompteurs << NbPartDetectee << endl;
    for (i=0;i<NBD;i++) SauvCompteurs << Comptage[i] << endl;   
    SauvCompteurs.close();

    // Fin du programme
    printf("Fin de l'experience\n");
    disfin();
    return EXIT_SUCCESS;
}