From 556a2a2feacf753104cdb4569a8d729a1ca938f1 Mon Sep 17 00:00:00 2001
From: hrocho <hrocho@vodacionline.sk>
Date: Fri, 16 Feb 2018 16:13:09 +0100
Subject: [PATCH] wp

---
 calculations.cpp |  6 +++++-
 calculations.h   |  1 +
 catima.cpp       | 12 +++++++-----
 catima.h         |  2 +-
 catimac.pxd      |  3 ++-
 5 files changed, 16 insertions(+), 8 deletions(-)

diff --git a/calculations.cpp b/calculations.cpp
index de3f17e..eba95d4 100644
--- a/calculations.cpp
+++ b/calculations.cpp
@@ -556,6 +556,10 @@ double beta_from_T(double T){
     return sqrt(1.0-1.0/(gamma*gamma));
 }
 
+double p_from_T(double T, double M){
+    return sqrt(T*(T+2*atomic_mass_unit))*M;
+}
+
 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);
@@ -566,7 +570,7 @@ double energy_straggling_firsov(double z1,double energy, double z2, double m2){
 double angular_scattering_variance(Projectile &p, Target &t){
     if(p.T<=0)return 0.0;
     double e=p.T;
-    double _p = sqrt(e*(e+2*atomic_mass_unit))*p.A;
+    double _p = p_from_T(e,p.A);
     double beta = _p/((e+atomic_mass_unit)*p.A);
     double lr = radiation_length(t.Z,t.A);
     return 198.81 * pow(p.Z,2)/(lr*pow(_p*beta,2));
diff --git a/calculations.h b/calculations.h
index a9b0449..119ea2b 100644
--- a/calculations.h
+++ b/calculations.h
@@ -174,6 +174,7 @@ namespace catima{
     //helper
     double gamma_from_T(double T);
     double beta_from_T(double T);
+    double p_from_T(double T, double M);
     std::complex<double> lngamma( const std::complex<double> &z );
     std::complex<double> hyperg(const std::complex<double> &a,
                                                 const std::complex<double> &b,
diff --git a/catima.cpp b/catima.cpp
index 3592c20..dcabd09 100644
--- a/catima.cpp
+++ b/catima.cpp
@@ -418,18 +418,20 @@ double calculate_tof_from_E(Projectile p, double Eout, const Material &t, const
     return res;
 }
 
-double w_magnification(Projectile p, double Ein, const Material &t, const Config &c){
-    double res = 1.0;
+std::pair<double,double> w_magnification(Projectile p, double Ein, const Material &t, const Config &c){
+    std::pair<double, double> res{1.0,1.0};
     if(t.density()<= 0.0 || t.thickness()<=0){
         return res;
     }
-    std::vector<double> energies{0.99*Ein, Ein, 1.1*Ein};
+    std::vector<double> energies{0.99*Ein, Ein, 1.01*Ein};
     auto eres = energy_out(p,energies,t,c);
     if(eres[0]>0.0 && eres[1]>0.0 && eres[2]>0.0){
-        res = energies[1]*(eres[2]-eres[0])/(eres[1]*(energies[2]-energies[0]));
+        res.first = energies[1]*(eres[2]-eres[0])/(eres[1]*(energies[2]-energies[0]));
+        res.second = p_from_T(energies[1],p.A)*(p_from_T(eres[2],p.A)-p_from_T(eres[0],p.A))/( p_from_T(eres[1],p.A)*( p_from_T(energies[2],p.A)-p_from_T(energies[0],p.A) ) );
         }
     else {
-        res = 0.0;
+        res.first = 0.0;
+        res.second = 0.0;
     }
     return res;
 }
diff --git a/catima.h b/catima.h
index 8608d27..bc5a249 100644
--- a/catima.h
+++ b/catima.h
@@ -212,7 +212,7 @@ namespace catima{
     /**
      * returns energy magnification after passing material t
      */
-    double w_magnification(Projectile p, double Ein, const Material &t, const Config &c=default_config);
+    std::pair<double,double> w_magnification(Projectile p, double Ein, const Material &t, const Config &c=default_config);
 
     class DataPoint;
     /**
diff --git a/catimac.pxd b/catimac.pxd
index ba77aba..f7f8e00 100644
--- a/catimac.pxd
+++ b/catimac.pxd
@@ -6,6 +6,7 @@
 """
 
 from libcpp.vector cimport vector
+from libcpp.pair cimport pair
 from libcpp cimport bool
 
 cdef extern from "catima/structures.h" namespace "catima":
@@ -87,7 +88,7 @@ cdef extern from "catima/catima.h" namespace "catima":
     cdef Result calculate(Projectile &p, const Material &t, const Config &c);
     cdef MultiResult calculate(Projectile &p, const Layers &layers, const Config &c);
 
-    cdef double w_magnification(Projectile p, double E, const Material &t, const Config &c);
+    cdef pair[double,double] w_magnification(Projectile p, double E, const Material &t, const Config &c);
 
 cdef extern from "catima/calculations.h" namespace "catima":
     cdef double bethek_lindhard(const Projectile &p);