From b4516ecd05592f6c8c5f46c042c3651e944c46b7 Mon Sep 17 00:00:00 2001 From: hrocho Date: Mon, 15 Jan 2018 15:11:51 +0100 Subject: [PATCH] mean charge calculation --- CMakeLists.txt | 5 ++ build_config.in | 1 + calculations.cpp | 92 +++++++++++++++++++++++++++++++++---- calculations.h | 28 +++++++++++ catima.pyx | 23 ++++++++++ catimac.pxd | 7 +++ config.h | 5 +- constants.h | 2 + tests/test_calculations.cpp | 39 ++++++++++++++++ 9 files changed, 192 insertions(+), 10 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 251cbf6..6ffe505 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -10,6 +10,7 @@ option(GSL_INTEGRATION "use GSL integration" ON) option(GENERATE_DATA "make data tables generator" OFF) option(THIN_TARGET_APPROXIMATION "thin target approximation" ON) option(DOCS "build documentation (requires doxygen)" OFF) +option(GLOBAL "build with global, sources are required" OFF) ######## build type ############ set(CMAKE_BUILD_TYPE Release) @@ -58,6 +59,10 @@ configure_file("${PROJECT_SOURCE_DIR}/init.sh.in" ############### main build ########################### file(GLOB SOURCES *.cpp) +if(GLOBAL) +file(GLOB GLOBAL_SOURCES global/*.c) +LIST (APPEND SOURCES ${GLOBAL_SOURCES}) +endif(GLOBAL) file(GLOB HEADERS *.h) add_library(catima SHARED ${SOURCES}) add_library(catima_static STATIC ${SOURCES}) diff --git a/build_config.in b/build_config.in index 1186513..87d7280 100644 --- a/build_config.in +++ b/build_config.in @@ -3,6 +3,7 @@ #cmakedefine THIN_TARGET_APPROXIMATION #cmakedefine GSL_INTEGRATION +#cmakedefine GLOBAL #endif diff --git a/calculations.cpp b/calculations.cpp index f92c86e..9c30527 100644 --- a/calculations.cpp +++ b/calculations.cpp @@ -8,6 +8,13 @@ #include "catima/generated_LS_coeff.h" #include "catima/nucdata.h" #include "catima/storage.h" +#ifdef GLOBAL +extern "C" +{ + #include "global/globallib.h" +} +#endif + namespace catima{ @@ -659,32 +666,43 @@ double z_effective(const Projectile &p,const Target &t, const Config &c){ if(c.z_effective == z_eff_type::pierce_blann){ return z_eff_Pierce_Blann(p.Z, beta); } - - if(c.z_effective == z_eff_type::anthony_landorf){ + else if(c.z_effective == z_eff_type::anthony_landorf){ return z_eff_Anthony_Landford(p.Z, beta, t.Z); } - - if(c.z_effective == z_eff_type::hubert){ + + else if(c.z_effective == z_eff_type::hubert){ return z_eff_Hubert(p.Z, p.T, t.Z); } - - return 0.0; + else if(c.z_effective == z_eff_type::winger){ + return z_eff_Winger(p.Z, beta, t.Z); + } + else if(c.z_effective == z_eff_type::global){ + return z_eff_global(p.Z, p.T, t.Z); + } + else if(c.z_effective == z_eff_type::schiwietz){ + return z_eff_Schiwietz(p.Z, beta, t.Z); + } + else{ + return 0.0; + } } double z_eff_Pierce_Blann(double z, double beta){ - return z*(1.0-exp(-130.1842*beta/pow(z,2.0/3.0))); + return z*(1.0-exp(-0.95*fine_structure_inverted*beta/pow(z,2.0/3.0))); } double z_eff_Anthony_Landford(double pz, double beta, double tz){ double B = 1.18-tz*(7.5e-03 - 4.53e-05*tz); double A = 1.16-tz*(1.91e-03 - 1.26e-05*tz); - return pz*(1.0-exp(-137.035999139*B*beta/pow(pz,2.0/3.0))*A); + return pz*(1.0-(A*exp(-fine_structure_inverted*B*beta/pow(pz,2.0/3.0)))); } double z_eff_Hubert(double pz, double E, double tz){ double lntz = log(tz); double x1,x2,x3,x4; + if(E<2.5) + return 0.0; if(tz == 4){ x1 = 2.045 + 2.0*exp(-0.04369*pz); x2 = 7.0; @@ -703,10 +721,66 @@ double z_eff_Hubert(double pz, double E, double tz){ x3 = 0.314 + 0.01072*lntz; x4 = 0.5218 + 0.02521*lntz; } - + return pz*(1-x1*exp(-x2*catima::power(E,x3)*catima::power(pz,-x4))); } +double z_eff_Winger(double pz, double beta, double tz){ + double c0,c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13; + double x, lnz, lnzt, a0,a1,a2,a3,a4; + + c0 = 0.4662; + c1 = 0.5491; + c2 = 0.7028; + c3 = 0.1089; + c4 = 0.001644; + c5 = 0.5155; + c6 = 0.05633; + c7 = 0.005447; + c8 = 0.8795; + c9 = 1.091; + c10= 0.0008261; + c11= 2.848; + c12= 0.2442; + c13= 0.00009293; + + x = beta /0.012 /pow(pz,0.45); + lnz =log(pz); + lnzt=log(tz); + + a0 = -c0; + a1 = -c1 * exp( c2 *lnz - c3 *lnz*lnz +c4*lnz*lnz*lnz -c5*lnzt + c6 *lnzt*lnzt); + a2 = c7 * exp( c8 *lnz - c9 *lnzt); + a3 = -c10 * exp( c11*lnz - c12*lnz*lnz*lnz); + a4 = -c13; + + return pz * (1. - exp(a0 +a1*x +a2*x*x +a3*x*x*x +a4*x*x*x*x)); + } + +double z_eff_global(double pz, double E, double tz){ + if(E>2000) + return pz; + else + #ifdef GLOBAL + return global_qmean(pz, tz, E); + #else + return -1; + #endif +} + +double z_eff_Schiwietz(double pz, double beta, double tz){ + double scaled_velocity; + double c1, c2; + double x; + + scaled_velocity = catima::power(pz,-0.543)*beta/bohr_velocity; + c1 = 1-0.26*exp(-tz/11.0)*exp(-(tz-pz)*(tz-pz)/9); + c2 = 1+0.030*scaled_velocity*log(tz); + x = c1*catima::power(scaled_velocity/c2/1.54,1+(1.83/pz)); + return pz*((8.29*x) + (x*x*x*x))/((0.06/x) + 4 + (7.4*x) + (x*x*x*x) ); + +} + std::complex hyperg(const std::complex &a, const std::complex &b, const std::complex &z){ diff --git a/calculations.h b/calculations.h index 2373611..f612d59 100644 --- a/calculations.h +++ b/calculations.h @@ -134,7 +134,35 @@ namespace catima{ * @return effective charge */ double z_eff_Hubert(double pz, double E, double tz); + + /** + * calculates effective charge + * @param pz - proton number of projectile + * @param beta - velocity of projectile + * @param tz - proton number of target material + * @return effective charge + */ + double z_eff_Winger(double pz, double beta, double tz); + + /** + * calculates effective charge + * @param pz - proton number of projectile + * @param beta - velocity of projectile + * @param tz - proton number of target material + * @return effective charge + */ + double z_eff_global(double pz, double E, double tz); + /** + * calculates effective charge + * @param pz - proton number of projectile + * @param beta - velocity of projectile + * @param tz - proton number of target material + * @return effective charge + */ + double z_eff_Schiwietz(double pz, double beta, double tz); + + //helper double gamma_from_T(double T); diff --git a/catima.pyx b/catima.pyx index 62368c9..9dc9360 100644 --- a/catima.pyx +++ b/catima.pyx @@ -263,6 +263,9 @@ class z_eff_type(IntEnum): pierce_blann = 1 anthony_landorf = 2 hubert = 3 + winger = 4 + schiwietz = 5 + global_code = 6 class skip_calculation(IntEnum): skip_none = 0 @@ -416,7 +419,27 @@ def z_effective(Projectile p, Target t, Config c = default_config): def z_eff_Pierce_Blann(double z, double beta): return catimac.z_eff_Pierce_Blann(z,beta) + +def z_eff_Anthony_Landford(double pz, double beta, double tz): + return catimac.z_eff_Anthony_Landford(pz, beta, tz); +def z_eff_Hubert(double pz, double E, double tz): + return catimac.z_eff_Hubert(pz, E, tz); + +def z_eff_Winger(double pz, double beta, double tz): + return catimac.z_eff_Winger(pz, beta, tz); + +def z_eff_global(double pz, double E, double tz): + return catimac.z_eff_global(pz, E, tz); + +def z_eff_Schiwietz(double pz, double beta, double tz): + return catimac.z_eff_Schiwietz(pz, beta, tz); + +def gamma_from_T(double T): + return catimac.gamma_from_T(T); + +def beta_from_T(double T): + return catimac.beta_from_T(T); def get_data(Projectile projectile, Material material, Config config = default_config): data = catimac.get_data(projectile.cbase, material.cbase, config.cbase) diff --git a/catimac.pxd b/catimac.pxd index 20b9569..173ba8f 100644 --- a/catimac.pxd +++ b/catimac.pxd @@ -85,6 +85,13 @@ cdef extern from "catima/catima.h" namespace "catima": cdef extern from "catima/calculations.h" namespace "catima": cdef double z_effective(const Projectile &p, const Target &t, const Config &c); cdef double z_eff_Pierce_Blann(double z, double beta); + cdef double z_eff_Anthony_Landford(double pz, double beta, double tz); + cdef double z_eff_Hubert(double pz, double E, double tz); + cdef double z_eff_Winger(double pz, double beta, double tz); + cdef double z_eff_global(double pz, double E, double tz); + cdef double z_eff_Schiwietz(double pz, double beta, double tz); + cdef double gamma_from_T(double T); + cdef double beta_from_T(double T); cdef extern from "catima/constants.h" namespace "catima": int max_datapoints "catima::max_datapoints" diff --git a/config.h b/config.h index 39cd76e..59456db 100644 --- a/config.h +++ b/config.h @@ -14,7 +14,10 @@ namespace catima{ atima = 1, // the same as Pierce Blann pierce_blann = 1, anthony_landorf = 2, - hubert = 3 + hubert = 3, + winger = 4, + schiwietz = 5, + global = 6 }; /** diff --git a/constants.h b/constants.h index af37aa5..2aa1487 100644 --- a/constants.h +++ b/constants.h @@ -24,7 +24,9 @@ constexpr double electron_mass = 0.510998928; // MeV/c^2 constexpr double atomic_mass_unit = 931.4940954; // MeV/c^2 constexpr double classical_electron_radius = 2.8179403227; //fm constexpr double fine_structure = 1/137.035999139; +constexpr double fine_structure_inverted = 1/fine_structure; constexpr double c_light = 299.792458; //Mm/s +constexpr double bohr_velocity = 2.19 / c_light; // in c unit constexpr double dedx_constant = 0.3070749187; //4*pi*Na*me*c^2*r_e^2 //MeV cm^2 constexpr double domega2dx_constant = dedx_constant*electron_mass; //4*pi*Na*me*c^2*r_e^2 //MeV^2 cm^2 diff --git a/tests/test_calculations.cpp b/tests/test_calculations.cpp index 51ca74a..36891d4 100644 --- a/tests/test_calculations.cpp +++ b/tests/test_calculations.cpp @@ -4,6 +4,8 @@ using namespace std; using catima::approx; #include "catima/catima.h" +#include "catima/calculations.h" + bool rcompare(double a, double b,double eps){ if(fabs((a-b)/fabs(b))