#include<iostream>
#include<fstream>
#include<math.h>

using namespace std;

double func(double x0, double alpha , double x){
  double val1 = (x-x0)/alpha  * (x-x0)/alpha;
  double x1 = x/1000.;
  double val = x1*x1*x1*exp(-val1);
  return(val);
}

double func_energy( double x0, double alpha , double E, double mass){
   // x0 is still a velocity in m/s
   // alpha is still velocity spread m/s
   
   //engies in eV
   double Es = 0.5*mass * x0 * x0 * 1.0364269e-08; // energy spread 
   double dEs  = Es * (2.0*alpha)/x0     ;
  
   double arg = (sqrt(E) - sqrt(Es)) / dEs;

   double val = E * exp(- 4.0 * Es * arg * arg);

   return(val);

}


int main(void){


   double Emin = 0;
   double Emax = 4.; //eV -> Hartree
   int N = 1e06;
   double dE = (Emax-Emin)/(N);
 
   
   long double sum = 0.;
   long double sum2 = 0.;
   double E = 0.;
   double x = 0.;

   double alpha = 193.6;  // spread in m/s
   double x0    = 1932.3; //v0 m/s

   char hmass;
   
   cout << "GIVE ME VELOCITY SPREAD ALPHA   [m/s]: \n";
   cin >> alpha;
   cout << "GIVE ME STREAM VELOCITY V0      [m/s]: \n"; 
   cin >> x0;
   cout << "H2 (press h) or D2 (press d) \n";
   cin >> hmass;
   

   double mass = 0.;
   
   if(hmass == 'd'){ 
     mass = 2.01410178 * 2.0; // D2 in amu 
   }
   else{
     mass = 1.0078 * 2.0;     // H2
   }
   cout << "THE BEAM PARAMETERS ARE FOR " << hmass << "2 \n";
   cout << "Mass in [amu] is " << mass << "\n";


   for(int i = 0; i < N; i++){
     E = i*dE + Emin;
     x = sqrt(2.0*E/mass);
     //convert from atomicmass units to m/s
     x *= 9822.6952; // convert sqrt(eV / amu) to m/s

     double val = func(x0, alpha , x);
     sum  += val;
     sum2 += E*val;
//     cout << x << " " << val << "\n";
   }

//   cout << "INTEGRATION OF THE VELOCITY DISTRIBUTION : ";
//   cout << sum << " " << sum2 << " " << (sum2/sum) << " "  << (sum2/sum)*96.485341  << "\n";
//   cout << (sum2/sum) << " "  << (sum2/sum)*96.485341  << "\n";

   sum = sum2 = 0.;

   //intergrate ENERGY DISTRIBUTION   
   for(int i = 0; i < N; i++){
     E = i*dE + Emin;
     double val = func_energy( x0, alpha , E, mass);
     sum  += val;
     sum2 += E*val;
   }
   cout << "INTEGRATION OF THE ENERGY DISTRIBUTION   : ";
//   cout << sum << " " << sum2 << " " << (sum2/sum) << " "  << (sum2/sum)*96.485341  << "\n";
   cout << (sum2/sum) << " meV "  << (sum2/sum)*96.485341  << " KJ/mol\n";



}