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

using namespace std;

// Déclaration des paramètres physiques
#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

// Définition de la zone de travail en unité caractéristique
#define XMIN       -5
#define XMAX       5
#define YMIN       0
#define YMAX       10                  

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

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


//************************************************************* 
// Variables globales
//*************************************************************
double K;

//******************************************************************************
// 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 - Diffusion de He2+");
    name("Distance du noyau","X");
    name("Paramètre d'impact","Y");
    height(30);           // hauteur des caractères 
    hname(30);            // hauteur des labels d'axe
    htitle(40);           // hauteur du titre
    titlin("Experience de Rutherford",1);
}

//******************************************************************************
// Calcul du critère d'arrêt de la boucle de calcul
//******************************************************************************
bool Arret(double x, double y)
{
    bool Fin;
    
    Fin = false;
    if ((x < (XMIN+0.0005)) ||
        (x >= (XMAX - 0.001)) ||
        (y < YMIN) ||
        (y >= (YMAX - 0.001))) Fin = true;
    return Fin;    
}

//******************************************************************************
// Description du système différentiel pour RK4
//******************************************************************************
void Derivee(double y[], double dydt[])
{
    double d;
    
    d = (y[0]*y[0]+y[1]*y[1])*Distance(y[0],y[1]);
    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;
}

  
//******************************************************************************
// Programme principal
//******************************************************************************
int main(int argc, char *argv[])
{
    double lc, tc, k, kprime;
    double deltat;
    double x, y;
    double vx, vy;
    double b, DMin, DMinC,d;
    double x0, y0, v0,thetaT, thetaM, E0;
    double Yin[N], Yout[N];
    int i;
    char buffer[80];

    // En-tête
    cout << "Diffusion coulombienne d'une particule alpha par un noyau Au\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
    x0 = XMIN + 0.001;         // unité caractéristique       
    deltat = 0.001;            // unité caractéristique

    // Saisie de l'énergie cinétique de la particule alpha
    cout << "Energie cinetique de la particule alpha (en MeV) : " << flush;
    cin >> E0;
    v0 = sqrt(2*E0*1.6e-13/Malpha)/c;   // vitesse initiale en unité caractéristique
    
    // Saisie du paramètre d'impact
    cout << "Valeur du parametre d'impact (en unite caracteristique): " << flush;
    cin >> b;

    // 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,0.5,YMIN,YMAX,YMIN,0.5);
    title();
    
    // tracé du noyau à l'origine du repère rayon du noyau = 1 UC
    color("red");
    rlcirc(0.0,0.0,1);

    // initialisation des données
    x = x0; y = b;
    vx = v0; vy = 0.0;
    DMin = Distance(x,y); 
    rlstrt(x0,b); 

    // boucle de calcul et tracé
    color("blue");
    while (!Arret(x,y)) 
     {
       // Calcul de la trajectoire par RK4
       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];
       // calcul de la distance minimum particule-noyau
       d = Distance(x,y);
       if (d < DMin) DMin = d;
       // tracé du point
       rlconn(x,y);      
     }

    // Calcule de l'angle de diffusion et de la distance minimum
    thetaM = atan2(vy,vx);             // angle mesuré sur la simulation
    thetaT = 2*atan(K/(b*pow(v0,2)));  // angle calculé
    DMinC = 2*kprime/pow(v0*c,2);      // distance min en cas de choc frontal     
        
    // Affichage des paramètres dans le header de la page graphique
    height(40);
    sprintf(buffer," - Angle de diffusion: %.2lf degres  Distance mini du noyau: %.2lf fermis",
            thetaM*180/PI,DMin*lc*1e15);
    paghdr("Dominique Lefebvre - ",buffer,2,0); 
  
    // affichage de l'angle de diffusion
    printf("Angle de diffusion calcule: %.2lf degres\n",thetaT*180/PI);

    // Fin du programme
    disfin();
    return EXIT_SUCCESS;
}