From 3c0ce8c131cc76795a6a81623035909b03ff5b00 Mon Sep 17 00:00:00 2001
From: hrocho <hrocho@vodacionline.sk>
Date: Fri, 12 Jan 2018 01:25:23 +0100
Subject: [PATCH 1/2] energy straggling is like ATIMA now, befor from LS

---
 calculations.cpp | 20 ++++++++++++++------
 catima.cpp       | 17 +++++++++--------
 constants.h      |  2 +-
 3 files changed, 24 insertions(+), 15 deletions(-)

diff --git a/calculations.cpp b/calculations.cpp
index 7eb0720..125b761 100644
--- a/calculations.cpp
+++ b/calculations.cpp
@@ -532,8 +532,6 @@ double beta_from_T(double T){
     return sqrt(1.0-1.0/(gamma*gamma));
 }
 
- 
-
 double energy_straggling_firsov(double z1,double energy, double z2, double m2){
     double gamma = gamma_from_T(energy);
     double beta2=1.0-1.0/(gamma*gamma);
@@ -634,10 +632,20 @@ double precalculated_lindhard_X(const Projectile &p){
 double dedx_variance(Projectile &p, Target &t, const Config &c){
     double zp_eff = z_effective(p,t,c);
     double gamma = gamma_from_T(p.T);
-    double f = domega2dx_constant*pow(zp_eff,2)*t.Z*gamma*gamma;
-    //double X = bethek_lindhard_X(p); 
-    double X = precalculated_lindhard_X(p); 
-    return f*X/t.A;
+    double f = domega2dx_constant*pow(zp_eff,2)*t.Z/t.A;
+    double cor;
+    double beta2 = pow(beta_from_T(p.T),2);
+    cor = 24.89 * pow(t.Z,1.2324)/(electron_mass*1e6 * beta2)*
+			log( 2.0*electron_mass*1e6*beta2/(33.05*pow(t.Z,1.6364)));
+	cor = std::max(cor, 0.0 );
+	
+	//double X = bethek_lindhard_X(p);
+    double X = precalculated_lindhard_X(p);
+    X *= gamma*gamma;
+    if(p.T<1.0)
+		return std::min(f*(X+cor), energy_straggling_firsov(p.Z, p.T, t.Z,t.A));
+	else
+		return f*(X+cor);
 }
 
 double z_effective(const Projectile &p,const Target &t, const Config &c){
diff --git a/catima.cpp b/catima.cpp
index 95ce76d..e23bdc9 100644
--- a/catima.cpp
+++ b/catima.cpp
@@ -336,15 +336,16 @@ DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){
 
     double res;
     //calculate 1st point to have i-1 element ready for loop
-    res = integrator.integrate(fdedx,Ezero,energy_table(0));
-    res = p.A*res;
-    dp.range[0] = res;
-    res = da2dx(p,energy_table(0),t)*res;
-    dp.angular_variance[0] = res;
+    //res = integrator.integrate(fdedx,Ezero,energy_table(0));
+    //res = p.A*res;
+    dp.range[0] = 0.0;
     
-    res = integrator.integrate(fomega,Ezero,energy_table(0));
-    res = p.A*res;
-    dp.range_straggling[0]=res;
+    //res = da2dx(p,energy_table(0),t)*res;
+    dp.angular_variance[0] = 0.0;
+    
+    //res = integrator.integrate(fomega,Ezero,energy_table(0));
+    //res = p.A*res;
+    dp.range_straggling[0]=0.0;
 
     for(int i=1;i<max_datapoints;i++){
         res = p.A*integrator.integrate(fdedx,energy_table(i-1),energy_table(i));
diff --git a/constants.h b/constants.h
index 4353c67..af37aa5 100644
--- a/constants.h
+++ b/constants.h
@@ -3,7 +3,7 @@
 
 namespace catima {
 
-constexpr double Ezero = 9E-4; // lowest E to calculate, below taken as 0
+constexpr double Ezero = 1E-3; // lowest E to calculate, below taken as 0
 constexpr double logEmin = -3; // log of minimum energy
 constexpr double logEmax = 5.0;  // log of max energy
 constexpr int max_datapoints = 500; // how many datapoints between logEmin and logEmax

From 397473360fdbe60d3fc3f24012324fad3f907e37 Mon Sep 17 00:00:00 2001
From: hrocho <hrocho@vodacionline.sk>
Date: Fri, 12 Jan 2018 01:34:50 +0100
Subject: [PATCH 2/2] X precalculation fix

---
 calculations.cpp | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/calculations.cpp b/calculations.cpp
index 125b761..f92c86e 100644
--- a/calculations.cpp
+++ b/calculations.cpp
@@ -613,8 +613,9 @@ double precalculated_lindhard_X(const Projectile &p){
     double T = p.T;
     int z = (int)p.Z ;
     if(z>LS_MAX_Z)z=LS_MAX_Z;
-    //if(p.T<ls_coefficients::ls_energy_points[0])T=ls_coefficients::ls_energy_points[0];
-    if(p.T<ls_coefficients::ls_energy_table(0))T=ls_coefficients::ls_energy_table(0);
+    //if(p.T<ls_coefficients::ls_energy_table(0))T=ls_coefficients::ls_energy_table(0);
+    if(p.T<ls_coefficients::ls_energy_table(0))
+		return 1.0;
     double da = (p.A - element_atomic_weight(z))/element_atomic_weight(z);
     z = z-1;