1
0
Fork 0
mirror of https://github.com/gwm17/catima.git synced 2024-11-22 18:28:51 -05:00

refactoring

This commit is contained in:
hrocho 2019-05-13 20:51:09 +02:00
parent 955334a0b6
commit 7871489ba6
14 changed files with 176 additions and 214 deletions

View File

@ -113,31 +113,31 @@ if(PYTHON_MODULE)
$<INSTALL_INTERFACE:include>)
target_link_libraries(pycatima PRIVATE catima)
find_program(CYTHON_EXECUTABLE
NAMES cython cython2 cython3 cython.bat
DOC "path to the cython executable"
)
if(NOT CYTHON_EXECUTABLE)
MESSAGE(SEND_ERROR "Cython not found, it is required to build nurex python modules")
endif(NOT CYTHON_EXECUTABLE)
MESSAGE(STATUS "Cython found: " ${CYTHON_EXECUTABLE})
# find_program(CYTHON_EXECUTABLE
# NAMES cython cython2 cython3 cython.bat
# DOC "path to the cython executable"
# )
# if(NOT CYTHON_EXECUTABLE)
# MESSAGE(SEND_ERROR "Cython not found, it is required to build nurex python modules")
# endif(NOT CYTHON_EXECUTABLE)
# MESSAGE(STATUS "Cython found: " ${CYTHON_EXECUTABLE})
### build libraries string
foreach(entry ${EXTRA_LIBS} ${GSL_LIBRARIES} catima)
LIST (APPEND EXTRA_PYTHON_LIBS \"${entry}\")
endforeach(entry ${EXTRA_LIBS} ${GSL_LIBRARIES} catima)
string (REPLACE ";" "," EXTRA_PYTHON_LIBS "${EXTRA_PYTHON_LIBS}")
# foreach(entry ${EXTRA_LIBS} ${GSL_LIBRARIES} catima)
# LIST (APPEND EXTRA_PYTHON_LIBS \"${entry}\")
# endforeach(entry ${EXTRA_LIBS} ${GSL_LIBRARIES} catima)
# string (REPLACE ";" "," EXTRA_PYTHON_LIBS "${EXTRA_PYTHON_LIBS}")
# if(THREADS)
# set (CYTHON_DEFINES "-DUSE_THREADS=1")
# endif(THREADS)
### insert libraries string and create setup.py
FILE(COPY catimac.pxd catima.pyx DESTINATION ${PROJECT_BINARY_DIR})
set(CATIMA_LIB ${CMAKE_SHARED_LIBRARY_PREFIX}catima${CMAKE_SHARED_LIBRARY_SUFFIX})
configure_file(${CMAKE_CURRENT_SOURCE_DIR}/setup.py.in ${PROJECT_BINARY_DIR}/setup.py)
add_custom_target(target_python ALL DEPENDS catima)
add_custom_command(TARGET target_python
COMMAND ${PYTHON_EXECUTABLE} ${PROJECT_BINARY_DIR}/setup.py build_ext ${CYTHON_DEFINES} -i
)
# FILE(COPY catimac.pxd catima.pyx DESTINATION ${PROJECT_BINARY_DIR})
# set(CATIMA_LIB ${CMAKE_SHARED_LIBRARY_PREFIX}catima${CMAKE_SHARED_LIBRARY_SUFFIX})
# configure_file(${CMAKE_CURRENT_SOURCE_DIR}/setup.py.in ${PROJECT_BINARY_DIR}/setup.py)
# add_custom_target(target_python ALL DEPENDS catima)
# add_custom_command(TARGET target_python
# COMMAND ${PYTHON_EXECUTABLE} ${PROJECT_BINARY_DIR}/setup.py build_ext ${CYTHON_DEFINES} -i
# )
endif(PYTHON_MODULE )
########## Sub Directories ###########

View File

@ -26,10 +26,9 @@ bool operator==(const Config &a, const Config&b){
}
double dedx(Projectile &p, double T, const Material &mat, const Config &c){
double dedx(Projectile &p, const Material &mat, const Config &c){
double sum = 0;
if(T<=0)return 0.0;
p.T = T;
if(p.T<=0)return 0.0;
sum += dedx_n(p,mat);
double se=0;
if(p.T<=10){
@ -74,18 +73,18 @@ double da2dx(Projectile &p, double T, const Material &mat, const Config &c){
}
double range(Projectile &p, double T, const Material &t, const Config &c){
double range(Projectile &p, const Material &t, const Config &c){
auto& data = _storage.Get(p,t,c);
//Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num);
spline_type range_spline = get_range_spline(data);
return range_spline(T);
return range_spline(p.T);
}
double dedx_from_range(Projectile &p, double T, const Material &t, const Config &c){
double dedx_from_range(Projectile &p, const Material &t, const Config &c){
auto& data = _storage.Get(p,t,c);
//Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num);
spline_type range_spline = get_range_spline(data);
return p.A/range_spline.derivative(T);
return p.A/range_spline.derivative(p.T);
}
std::vector<double> dedx_from_range(Projectile &p, const std::vector<double> &T, const Material &t, const Config &c){
@ -177,11 +176,11 @@ double energy_out(double T, double thickness, const Interpolator &range_spline){
return -1;
}
double energy_out(Projectile &p, double T, const Material &t, const Config &c){
double energy_out(Projectile &p, const Material &t, const Config &c){
auto& data = _storage.Get(p,t,c);
//Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num);
spline_type range_spline = get_range_spline(data);
return energy_out(T,t.thickness(),range_spline);
return energy_out(p.T,t.thickness(),range_spline);
}
std::vector<double> energy_out(Projectile &p, const std::vector<double> &T, const Material &t, const Config &c){
@ -306,74 +305,12 @@ Result calculate(double pa, int pz, double T, double ta, double tz, double thick
return calculate(p(T),m);
}
std::vector<double> calculate_range(Projectile p, const Material &t, const Config &c){
double res;
std::vector<double>values;
values.reserve(max_datapoints);
auto fdedx = [&](double x)->double{return 1.0/dedx(p,x,t);};
//calculate 1st point to have i-1 element ready for loop
res = integrator.integrate(fdedx,Ezero,energy_table(0));
res = p.A*res;
values.push_back(res);
for(int i=1;i<max_datapoints;i++){
res = integrator.integrate(fdedx,energy_table(i-1),energy_table(i));
res = p.A*res;
res += values[i-1];
values.push_back(res);
}
return values;
}
std::vector<double> calculate_range_straggling(Projectile p, const Material &t, const Config &c){
double res;
std::vector<double>values;
values.reserve(max_datapoints);
auto function = [&](double x)->double{return 1.0*domega2dx(p,x,t)/pow(dedx(p,x,t),3);};
//auto function = [&](double x)->double{
//double de = dedx(p,x,t);
//return 1.0*domega2dx(p,x,t)/(de*de*de);
//};
//calculate 1st point to have i-1 element ready for loop
res = integrator.integrate(function,Ezero,energy_table(0));
res = p.A*res;
values.push_back(res);
for(int i=1;i<max_datapoints;i++){
res = integrator.integrate(function,energy_table(i-1),energy_table(i));
res = p.A*res;
res += values[i-1];
values.push_back(res);
}
return values;
}
std::vector<double> calculate_da2dx(Projectile p, const Material &t, const Config &c){
double res;
std::vector<double>values;
values.reserve(max_datapoints);
//auto function = [&](double x)->double{return p.A*da2dx(p,x,t)/dedx(p,x,t);};
auto function = [&](double x)->double{return 1.0/dedx(p,x,t);};
res = integrator.integrate(function,Ezero,energy_table(0));
res = p.A*da2dx(p,energy_table(0),t)*res;
values.push_back(res);
for(int i=1;i<max_datapoints;i++){
res = integrator.integrate(function,energy_table(i-1),energy_table(i));
res = p.A*da2dx(p,energy_table(i),t)*res;
res += values[i-1];
values.push_back(res);
}
return values;
}
std::vector<double> calculate_tof(Projectile p, const Material &t, const Config &c){
double res;
std::vector<double> values;
values.reserve(max_datapoints);
auto function = [&](double x)->double{return 1.0/(dedx(p,x,t)*beta_from_T(x));};
auto function = [&](double x)->double{return 1.0/(dedx(p(x),t,c)*beta_from_T(x));};
res = integrator.integrate(function,Ezero,energy_table(0));
res = res*10.0*p.A/(c_light*t.density());
values.push_back(res);
@ -392,11 +329,11 @@ DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){
dp.range_straggling.resize(max_datapoints);
dp.angular_variance.resize(max_datapoints);
auto fdedx = [&](double x)->double{
return 1.0/dedx(p,x,t,c);
return 1.0/dedx(p(x),t,c);
};
auto fomega = [&](double x)->double{
//return 1.0*domega2dx(p,x,t)/pow(dedx(p,x,t),3);
return domega2dx(p,x,t,c)/catima::power(dedx(p,x,t,c),3);
//return 1.0*domega2dx(p,x,t)/pow(dedx(p(x),t),3);
return domega2dx(p,x,t,c)/catima::power(dedx(p(x),t,c),3);
};
double res;
@ -434,11 +371,11 @@ DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){
dp.range_straggling.resize(max_datapoints);
dp.angular_variance.resize(max_datapoints);
auto fdedx = [&](double x)->double{
return 1.0/dedx(p,x,t,c);
return 1.0/dedx(p(x),t,c);
};
auto fomega = [&](double x)->double{
//return 1.0*domega2dx(p,x,t)/pow(dedx(p,x,t),3);
return domega2dx(p,x,t,c)/catima::power(dedx(p,x,t,c),3);
//return 1.0*domega2dx(p,x,t)/pow(dedx(p(x),t,c),3);
return domega2dx(p,x,t,c)/catima::power(dedx(p(x),t,c),3);
};
double res=0.0;
@ -471,7 +408,7 @@ double calculate_tof_from_E(Projectile p, double Eout, const Material &t, const
double res;
//double beta_in = beta_from_T(p.T);
//double beta_out = beta_from_T(Eout);
auto function = [&](double x)->double{return 1.0/(dedx(p,x,t)*beta_from_T(x));};
auto function = [&](double x)->double{return 1.0/(dedx(p(x),t,c)*beta_from_T(x));};
res = integrator.integrate(function,Eout,p.T);
res = res*10.0*p.A/(c_light*t.density());
return res;

View File

@ -37,12 +37,12 @@ namespace catima{
* @param mat - Material
* @return dEdx
*/
double dedx(Projectile &p, double T, const Material &mat, const Config &c=default_config);
double dedx(Projectile &p, const Material &mat, const Config &c=default_config);
/**
* calculate energy loss straggling variance for projectile-Material combination
* @param p - Projectile
* @param mat - Material
* @param t - Material
* @return dOmega^2/dx
*/
double domega2dx(Projectile &p, double T, const Material &t, const Config &c=default_config);
@ -55,20 +55,18 @@ namespace catima{
/**
* returns the range of the Projectile in Material calculated from range spline
* @param p - Projectile
* @param T - energy in MeV/u
* @param mat - Material
* @param t - Material
* @return range
*/
double range(Projectile &p, double T, const Material &t, const Config &c=default_config);
double range(Projectile &p, const Material &t, const Config &c=default_config);
/**
* returns the dEdx calculated from range spline as derivative
* @param p - Projectile
* @param T - energy in MeV/u
* @param mat - Material
* @param t - Material
* @return range
*/
double dedx_from_range(Projectile &p, double T, const Material &t, const Config &c=default_config);
double dedx_from_range(Projectile &p, const Material &t, const Config &c=default_config);
/**
* returns the dEdx calculated from range spline as derivative
@ -151,7 +149,7 @@ namespace catima{
* @param T - incoming energy
* @return outcoming energy after the material in Mev/u
*/
double energy_out(Projectile &p, double T, const Material &t, const Config &c=default_config);
double energy_out(Projectile &p, const Material &t, const Config &c=default_config);
/**
* calculates outcoming energy
@ -194,11 +192,8 @@ namespace catima{
return calculate(p, layers, c);
}
/// the following functions are used to calculates array of data points for whole range of energies
/// usually used to construct splines
std::vector<double> calculate_range(Projectile p, const Material &t, const Config &c=default_config);
std::vector<double> calculate_range_straggling(Projectile p, const Material &t, const Config &c=default_config);
std::vector<double> calculate_angular_variance(Projectile p, const Material &t, const Config &c=default_config);
/// this calculate tof spline, at the moment it is not used
std::vector<double> calculate_tof(Projectile p, const Material &t, const Config &c=default_config);
/**

View File

@ -399,11 +399,12 @@ def projectile_range(Projectile projectile, Material material, energy = None, Co
if(isinstance(energy,numpy.ndarray)):
res = numpy.empty(energy.size)
for i,e in enumerate(energy):
res[i] = catimac.range(projectile.cbase, e, material.cbase, config.cbase)
projectile.T(e)
res[i] = catimac.range(projectile.cbase, material.cbase, config.cbase)
return res
if(energy is None):
energy = projectile.T()
return catimac.range(projectile.cbase, energy, material.cbase, config.cbase);
elif(energy is not None):
projectile.T(energy)
return catimac.range(projectile.cbase, material.cbase, config.cbase);
def dedx_from_range(Projectile projectile, Material material, energy = None, Config config = default_config):
if(isinstance(energy,numpy.ndarray)):
@ -411,9 +412,9 @@ def dedx_from_range(Projectile projectile, Material material, energy = None, Con
return numpy.asarray(res);
if(isinstance(energy,list)):
return catimac.dedx_from_range(projectile.cbase, <vector[double]>energy, material.cbase, config.cbase)
if(energy is None):
energy = projectile.T()
return catimac.dedx_from_range(projectile.cbase, <double>energy, material.cbase, config.cbase);
if(energy is not None):
projectile.T(energy)
return catimac.dedx_from_range(projectile.cbase, material.cbase, config.cbase);
def domega2de(Projectile projectile, Material material, energy = None, Config config = default_config):
if(isinstance(energy,numpy.ndarray)):
@ -439,11 +440,12 @@ def dedx(Projectile projectile, Material material, energy = None, Config config
if(isinstance(energy,numpy.ndarray)):
res = numpy.empty(energy.size)
for i,e in enumerate(energy):
res[i] = catimac.dedx(projectile.cbase, e, material.cbase, config.cbase)
projectile.T(e)
res[i] = catimac.dedx(projectile.cbase, material.cbase, config.cbase)
return res
if(energy is None):
energy = projectile.T()
return catimac.dedx(projectile.cbase, energy, material.cbase, config.cbase)
if(energy is not None):
projectile.T(energy)
return catimac.dedx(projectile.cbase, material.cbase, config.cbase)
def domega2dx(Projectile projectile, Material material, energy = None, Config config = default_config):
if(isinstance(energy,numpy.ndarray)):

View File

@ -72,10 +72,10 @@ cdef extern from "catima/config.h" namespace "catima":
char calculation;
cdef extern from "catima/catima.h" namespace "catima":
cdef double dedx(Projectile &p, double T, const Material &t,const Config &c)
cdef double dedx(Projectile &p, const Material &t,const Config &c)
cdef double domega2dx(Projectile &p, double T, const Material &mat, const Config &c)
cdef double range(Projectile &p, double T, const Material &t, const Config &c);
cdef double dedx_from_range(Projectile &p, double T, const Material &t, const Config &c);
cdef double range(Projectile &p, const Material &t, const Config &c);
cdef double dedx_from_range(Projectile &p, const Material &t, const Config &c);
cdef vector[double] dedx_from_range(Projectile &p, vector[double] &T, const Material &t, const Config &c);
cdef double energy_out(Projectile &p, double T, const Material &t, const Config &c);
cdef vector[double] energy_out(Projectile &p, vector[double] &T, const Material &t, const Config &c);

View File

@ -82,19 +82,19 @@ namespace catima{
};
inline void set_config_lowenergy(Config c, low_energy_types lt){
c.corrections = c.corrections & (lt<<2);
c.calculation = c.calculation & (lt<<2);
}
inline unsigned char config_lowenergy(const Config c){
return (c.corrections>>2) & 0x7;
return (c.calculation>>2) & 0x7;
}
inline void set_config_omega(Config c, omega_types ot){
c.corrections = c.corrections & ot;
c.calculation = c.calculation & ot;
}
inline unsigned char config_omega(const Config c){
return c.corrections & 0x3;
return c.calculation & 0x3;
}

View File

@ -1,8 +1,8 @@
#include <exception>
#include <stdexcept>
#include <pybind11/pybind11.h>
#include <pybind11/operators.h>
#include <pybind11/numpy.h>
#include <pybind11/stl.h>
#include "catima/catima.h"
#include "catima/srim.h"
#include "catima/nucdata.h"
@ -56,7 +56,7 @@ py::list get_data(Projectile& p, const Material &m, const Config& c=default_conf
py::list ran;
py::list rans;
py::list av;
for(double e:data.range){ran.append(py::float_(e));}
for(double e:data.range){ran.append(e);}
for(double e:data.range_straggling)rans.append(e);
for(double e:data.angular_variance)av.append(e);
r.append(ran);
@ -65,11 +65,12 @@ py::list get_data(Projectile& p, const Material &m, const Config& c=default_conf
return r;
}
Material py_make_material(py::list d, double density=0.0, double ipot=0.0, double mass=0.0){
Material py_make_material(py::list d, double density=0.0, double thickness=0.0, double ipot=0.0, double mass=0.0){
Material m;
if(density>0.0)m.density(density);
if(ipot>0.0)m.I(ipot);
if(mass>0.0)m.M(mass);
if(thickness>0.0)m.thickness(thickness);
for(int i=0;i<d.size();i++){
py::list e(d[i]);
if(e.size() != 3)throw std::invalid_argument("invalid Material constructor argument");
@ -99,7 +100,9 @@ PYBIND11_MODULE(pycatima,m){
py::class_<Material>(m,"Material")
.def(py::init<>(),"constructor")
.def(py::init(&py_make_material),"constructor", py::arg("elements"),py::arg("density")=0.0,py::arg("i_potential")=0.0, py::arg("mass")=0.0)
.def(py::init<const Material&>(),"constructor")
.def(py::init<double, int, double, double, double>(),"constructor", py::arg("A"),py::arg("Z"),py::arg("density")=0.0,py::arg("thickness")=0.0,py::arg("i_potential")=0.0)
.def(py::init(&py_make_material),"constructor", py::arg("elements"),py::arg("density")=0.0,py::arg("thickness")=0.0,py::arg("i_potential")=0.0, py::arg("mass")=0.0)
.def("add_element",&Material::add_element)
.def("ncomponents",&Material::ncomponents)
.def("density",py::overload_cast<>(&Material::density, py::const_), "get density")
@ -115,9 +118,16 @@ PYBIND11_MODULE(pycatima,m){
.def(py::init<>(),"constructor")
.def("add",&Layers::add)
.def("num",&Layers::num)
.def("__getitem__",&Layers::operator[], py::is_operator())
// .def("__getitem__",&Layers::operator[], py::is_operator())
.def("__getitem__",[](Layers &r, int i)->Material*
{
if(i>=r.num()){
throw std::invalid_argument("index out of range");}
return &r[i];
}, py::is_operator(),py::return_value_policy::automatic_reference)
.def("get",&Layers::operator[])
.def(py::self + py::self);
.def(py::self + py::self)
.def("__add__",[](const Layers s, const Material& m){return s+m;});
py::class_<Result>(m,"Result")
.def(py::init<>(),"constructor")
@ -152,6 +162,27 @@ PYBIND11_MODULE(pycatima,m){
.def(py::init<>(),"constructor")
.def_readwrite("total_result", &MultiResult::total_result)
.def_readwrite("results", &MultiResult::results)
// .def_readwrite("Eout",&MultiResult::total_result.Eout)
.def("__getitem__",[](MultiResult &r, int i){
return py::cast(r.results[i]);
})
.def("__getattr__",[](MultiResult &r, std::string& k){
if(k.compare("Eout")==0){
return py::cast(r.total_result.Eout);
}
else if(k.compare("sigma_a")==0){
return py::cast(r.total_result.sigma_a);
}
else if(k.compare("tof")==0){
return py::cast(r.total_result.tof);
}
else if(k.compare("Eloss")==0){
return py::cast(r.total_result.Eloss);
}
else{
return py::cast(NULL);
}
})
.def("getJSON",[](const MultiResult &r){
py::dict d;
py::list p;
@ -235,16 +266,21 @@ PYBIND11_MODULE(pycatima,m){
});
m.def("srim_dedx_e",&srim_dedx_e);
m.def("sezi_dedx_e",&sezi_dedx_e, "sezi_dedx_e", py::arg("Projectile"), py::arg("Material"), py::arg("Config")=default_config);
m.def("calculate",py::overload_cast<Projectile&, const Material&, const Config&>(&calculate),"calculate",py::arg("Projectile"), py::arg("Material"), py::arg("Config")=default_config);
m.def("dedx_from_range",py::overload_cast<Projectile&, double, const Material&, const Config&>(&dedx_from_range),"calculate",py::arg("Projectile"), py::arg("energy") ,py::arg("Material"), py::arg("config")=default_config);
m.def("dedx_from_range",py::overload_cast<Projectile&, const std::vector<double>&, const Material&, const Config&>(&dedx_from_range),"calculate",py::arg("Projectile"), py::arg("energy") ,py::arg("Material"), py::arg("config")=default_config);
//m.def("calculate",py::overload_cast<Projectile&, const Material&, double, const Config&>(&calculate), "calculate",py::arg("Projectile"), py::arg("Material"), py::arg("energy"), py::arg("config")=default_config);
m.def("sezi_dedx_e",&sezi_dedx_e, "sezi_dedx_e", py::arg("projectile"), py::arg("material"), py::arg("config")=default_config);
m.def("calculate",py::overload_cast<Projectile&, const Material&, const Config&>(&calculate),"calculate",py::arg("projectile"), py::arg("material"), py::arg("config")=default_config);
m.def("calculate_layers",py::overload_cast<Projectile&, const Layers&, const Config&>(&calculate),"calculate_layers",py::arg("projectile"), py::arg("material"), py::arg("config")=default_config);
m.def("dedx_from_range",py::overload_cast<Projectile&, const Material&, const Config&>(&dedx_from_range),"calculate",py::arg("projectile") ,py::arg("material"), py::arg("config")=default_config);
m.def("dedx_from_range",py::overload_cast<Projectile&, const std::vector<double>&, const Material&, const Config&>(&dedx_from_range),"calculate",py::arg("projectile"), py::arg("energy") ,py::arg("material"), py::arg("config")=default_config);
m.def("dedx",py::overload_cast<Projectile&, const Material&, const Config&>(&dedx), "dedx",py::arg("projectile"), py::arg("material"), py::arg("config")=default_config);
m.def("range",py::overload_cast<Projectile&, const Material&, const Config&>(&range), "range",py::arg("projectile"), py::arg("material"), py::arg("config")=default_config);
m.def("energy_out",py::overload_cast<Projectile&, const std::vector<double>&, const Material&, const Config&>(&energy_out),"energy_out",py::arg("projectile"), py::arg("energy") ,py::arg("material"), py::arg("config")=default_config);
m.def("energy_out",py::overload_cast<Projectile&, const Material&, const Config&>(&energy_out),"energy_out",py::arg("projectile"), py::arg("material"), py::arg("config")=default_config);
m.def("get_material",py::overload_cast<int>(&get_material));
m.def("get_data",py::overload_cast<Projectile&, const Material&, const Config&>(get_data),"list of data",py::arg("projectile"),py::arg("material"),py::arg("Config")=default_config);
m.def("get_data",py::overload_cast<Projectile&, const Material&, const Config&>(get_data),"list of data",py::arg("projectile"),py::arg("material"),py::arg("config")=default_config);
m.def("catima_info",&catima_info);
m.def("storage_info",&storage_info);
m.def("get_energy_table",&get_energy_table);
m.def("energy_table",[](int i){return energy_table(i);});
m.attr("max_datapoints") = max_datapoints;
m.attr("max_storage_data") = max_storage_data;
m.attr("logEmin")=logEmin;

View File

@ -111,12 +111,6 @@ void Data::Add(const Projectile &p, const Material &t, const Config &c){
if(e==dp)return;
}
if(index==storage.end())index=storage.begin();
#if(0)
*index = dp;
index->range = calculate_range(p,t,c);
index->range_straggling = calculate_range_straggling(p,t,c);
index->angular_variance = calculate_angular_variance(p,t,c);
#else
*index = calculate_DataPoint(p,t,c);
#ifdef STORE_SPLINES
//index->range_spline = Interpolator(energy_table.values,index->range);
@ -126,8 +120,6 @@ void Data::Add(const Projectile &p, const Material &t, const Config &c){
index->range_straggling_spline = Interpolator(energy_table, index->range_straggling);
index->angular_variance_spline = Interpolator(energy_table, index->angular_variance);
#endif
#endif
index++;
}

View File

@ -42,7 +42,7 @@ Material::Material(std::initializer_list<std::array<double,3>>list,double _densi
}
}
Material::Material(double _a, int _z, double _rho, double _th):rho(_rho),th(_th){
Material::Material(double _a, int _z, double _rho, double _th, double _ipot):rho(_rho),th(_th),i_potential(_ipot){
add_element(_a,_z,1.0);
}

View File

@ -78,7 +78,7 @@ namespace catima{
* @param _rho - density of the material in g/cm3, default 0.0
* @param _th - thickness of the material in g/cm2, default 0.0
*/
Material(double _a, int _z, double _rho=0.0, double _th=0.0);
Material(double _a, int _z, double _rho=0.0, double _th=0.0, double _ipot = 0.0);
/**

View File

@ -1,7 +1,7 @@
import sys
sys.path.insert(0,"../build")
import unittest
import catima
import pycatima as catima
import math
class TestStructures(unittest.TestCase):
@ -18,10 +18,8 @@ class TestStructures(unittest.TestCase):
self.assertEqual(p.Q(),90)
p.T(1000)
self.assertEqual(p.T(),1000)
self.assertEqual(p(),1000)
p(500)
self.assertEqual(p.T(),500)
self.assertEqual(p(),500)
p = catima.Projectile(238,92,90, T=100)
self.assertEqual(p.T(),100)
@ -33,27 +31,27 @@ class TestStructures(unittest.TestCase):
mat.add_element(1,1,2)
self.assertEqual(mat.ncomponents(),2)
mat2 = catima.Material([12.01,6,1])
mat2 = catima.Material(12.01,6)
self.assertEqual(mat2.ncomponents(),1)
self.assertEqual(mat2.molar_mass(),12.01)
mat3 = catima.Material([12,6,1])
mat3 = catima.Material([[12,6,1]])
self.assertEqual(mat3.ncomponents(),1)
self.assertEqual(mat3.molar_mass(),12)
water = catima.Material([[1,1,2],[16,8,1]])
self.assertEqual(water.molar_mass(),18)
mat2 = catima.Material([0,6,1])
mat2 = catima.Material(0,6)
self.assertEqual(mat2.ncomponents(),1)
self.assertAlmostEqual(mat2.molar_mass(),12,1)
mat5 = catima.Material([0,6,1],density=1.9, thickness=0.5)
mat5 = catima.Material(0,6,density=1.9, thickness=0.5)
self.assertEqual(mat5.ncomponents(),1)
self.assertEqual(mat5.thickness(),0.5)
self.assertEqual(mat5.density(),1.9)
mat6 = catima.Material([0,6,1],density=1.9, thickness=0.5,i_potential=80.0)
mat6 = catima.Material(0,6,density=1.9, thickness=0.5,i_potential=80.0)
self.assertEqual(mat6.ncomponents(),1)
self.assertEqual(mat6.thickness(),0.5)
self.assertEqual(mat6.density(),1.9)
@ -61,7 +59,7 @@ class TestStructures(unittest.TestCase):
# copy
mat3.density(1.8)
matc = mat3.copy()
matc = catima.Material(mat3)
self.assertEqual(matc.ncomponents(),1)
self.assertEqual(matc.molar_mass(),12)
self.assertEqual(matc.density(),1.8)
@ -76,7 +74,7 @@ class TestStructures(unittest.TestCase):
self.assertEqual(m1.ncomponents(),1)
self.assertAlmostEqual(m1.density(),2.0,1)
m2 = catima.get_material(catima.material.WATER)
m2 = catima.get_material(catima.material.Water)
self.assertEqual(m2.ncomponents(),2)
self.assertAlmostEqual(m2.molar_mass(),18,1)
self.assertAlmostEqual(m2.density(),1.0,1)
@ -114,14 +112,15 @@ class TestStructures(unittest.TestCase):
self.assertAlmostEqual(mat[0].density(),2.0,1)
self.assertAlmostEqual(mat[2].thickness(),1.0,1)
self.assertAlmostEqual(mat[2].density(),1.8,1)
mat[2].thickness(1.2)
mat[2].density(1.9)
self.assertAlmostEqual(mat.materials[2].thickness(),1.2,1)
self.assertAlmostEqual(mat.materials[2].density(),1.9,1)
#self.assertAlmostEqual(mat.materials[0].thickness(),0.5,1)
#self.assertAlmostEqual(mat.materials[0].density(),2.0,1)
self.assertEqual(mat[3],None)
self.assertEqual(mat["a"],None)
a = mat[2]
a.thickness(1.2)
a.density(1.9)
self.assertAlmostEqual(mat[2].thickness(),1.2,1)
self.assertAlmostEqual(mat[2].density(),1.9,1)
self.assertEqual(mat.num(),3)
with self.assertRaises(ValueError):
mat[3]
mat2 = catima.Layers()
mat2.add(n2)
@ -143,14 +142,13 @@ class TestStructures(unittest.TestCase):
self.assertEqual(mats[4].thickness(),0.5)
def test_material_calculation(self):
water = catima.get_material(catima.material.WATER)
water = catima.get_material(catima.material.Water)
p = catima.Projectile(1,1)
p(1000)
res = catima.calculate(p,water)
res2 = catima.dedx(p,water)
self.assertAlmostEqual(res.dEdxi,2.23,1)
self.assertAlmostEqual(res["dEdxi"],2.23,1)
self.assertAlmostEqual(res.dEdxi,res2,3)
res = catima.calculate(p(500),water)
res2 = catima.dedx(p,water)
@ -171,14 +169,16 @@ class TestStructures(unittest.TestCase):
self.assertAlmostEqual(res.dEdxi,res2,3)
def test_config(self):
water = catima.get_material(catima.material.WATER)
water = catima.get_material(catima.material.Water)
water.density(1.0)
water.thickness(1.0)
p = catima.Projectile(1,1)
conf = catima.Config()
conf.dedx_straggling(catima.omega_type.bohr)
conf.calculation = catima.omega_types.bohr
conf2 = catima.Config()
conf2.dedx_straggling(catima.omega_type.atima)
conf2.calculation = catima.omega_types.atima
self.assertEqual(conf.calculation, 1)
self.assertEqual(conf2.calculation, 0)
p(1000)
res = catima.calculate(p,water,config=conf)
res2 = catima.calculate(p,water,config=conf2)
@ -194,7 +194,6 @@ class TestStructures(unittest.TestCase):
res = catima.calculate(p(1000),graphite)
res2 = catima.energy_out(p(1000),graphite)
self.assertAlmostEqual(res.Eout,997.077,1)
self.assertAlmostEqual(res["Eout"],997.077,1)
self.assertAlmostEqual(res.Eout,res2,3)
def test_eout_list(self):
@ -204,7 +203,7 @@ class TestStructures(unittest.TestCase):
energies = [100,500,1000]
res = catima.calculate(p(1000),graphite)
self.assertAlmostEqual(res.Eout,997.077,1)
res2 = catima.energy_out(p,graphite,energy=energies)
res2 = catima.energy_out(p,energies, graphite)
self.assertEqual(len(res2),len(energies))
self.assertAlmostEqual(res2[2], 997.077,1)
self.assertAlmostEqual(res2[0], catima.calculate(p(energies[0]),graphite).Eout ,1)
@ -215,16 +214,16 @@ class TestStructures(unittest.TestCase):
graphite.thickness(0.5)
p = catima.Projectile(12,6)
energies = [100,500,1000]
res2 = catima.dedx_from_range(p,graphite,energy=energies)
res2 = catima.dedx_from_range(p,energies, graphite)
self.assertEqual(len(res2),len(energies))
self.assertEqual(len(res2),3)
for i,e in enumerate(energies):
r = catima.dedx_from_range(p, graphite, energy=e)
r = catima.dedx_from_range(p(e), graphite)
self.assertAlmostEqual(res2[i], r, 0.1)
def test_layer_calculation(self):
p = catima.Projectile(12,6)
water = catima.get_material(catima.material.WATER)
water = catima.get_material(catima.material.Water)
water.thickness(10.0)
graphite = catima.get_material(6)
graphite.thickness(1.0)
@ -237,17 +236,17 @@ class TestStructures(unittest.TestCase):
self.assertEqual(len(res.results),2)
self.assertAlmostEqual(res.total_result.Eout,926.3,1)
self.assertAlmostEqual(res.total_result.sigma_a,0.00269,1)
self.assertAlmostEqual(res["Eout"],926.3,1)
self.assertAlmostEqual(res["sigma_a"],0.00269,4)
self.assertAlmostEqual(res["tof"],0.402,2)
self.assertAlmostEqual(res["Eloss"],884,0)
self.assertAlmostEqual(res.Eout,926.3,1)
self.assertAlmostEqual(res.sigma_a,0.00269,4)
self.assertAlmostEqual(res.tof,0.402,2)
self.assertAlmostEqual(res.Eloss,884,0)
self.assertAlmostEqual(res[0]["Eout"],932.24,0)
self.assertAlmostEqual(res[1]["Eout"],926.3,0)
self.assertAlmostEqual(res[0]["sigma_a"],0.00258,4)
self.assertAlmostEqual(res[1]["sigma_a"],0.000774,4)
self.assertAlmostEqual(res[0]["range"],107.1,0)
self.assertAlmostEqual(res[1]["range"],111.3,0)
self.assertAlmostEqual(res[0].Eout,932.24,0)
self.assertAlmostEqual(res[1].Eout,926.3,0)
self.assertAlmostEqual(res[0].sigma_a,0.0025,3)
self.assertAlmostEqual(res[1].sigma_a,0.000774,4)
self.assertAlmostEqual(res[0].range,107.1,0)
self.assertAlmostEqual(res[1].range,111.3,0)
def test_energy_table(self):
table = catima.get_energy_table()
@ -257,7 +256,7 @@ class TestStructures(unittest.TestCase):
def test_storage(self):
p = catima.Projectile(12,6)
water = catima.get_material(catima.material.WATER)
water = catima.get_material(catima.material.Water)
water.thickness(10.0)
graphite = catima.get_material(6)
graphite.thickness(1.0)
@ -270,33 +269,33 @@ class TestStructures(unittest.TestCase):
res = catima.calculate(p(et[10]),water)
self.assertAlmostEqual(res.range,data[0][10],6)
self.assertAlmostEqual(catima.projectile_range(p,water),data[0][10],6)
self.assertAlmostEqual(catima.range(p,water),data[0][10],6)
#self.assertAlmostEqual(catima.domega2de(p,water),data[1][10],6)
res = catima.calculate(p(et[100]),water)
self.assertAlmostEqual(res.range,data[0][100],6)
self.assertAlmostEqual(catima.projectile_range(p,water),data[0][100],6)
self.assertAlmostEqual(catima.range(p,water),data[0][100],6)
#self.assertAlmostEqual(catima.domega2de(p,water),data[1][100],6)
res = catima.calculate(p(et[200]),water)
self.assertAlmostEqual(res.range,data[0][200],6)
self.assertAlmostEqual(catima.projectile_range(p,water),data[0][200],6)
self.assertAlmostEqual(catima.range(p,water),data[0][200],6)
#self.assertAlmostEqual(catima.domega2de(p,water),data[1][200],6)
res = catima.calculate(p(et[401]),water)
self.assertAlmostEqual(res.range,data[0][401],6)
self.assertAlmostEqual(catima.projectile_range(p,water),data[0][401],6)
self.assertAlmostEqual(catima.range(p,water),data[0][401],6)
#self.assertAlmostEqual(catima.domega2de(p,water),data[1][401],6)
def test_python_storage_access(self):
p = catima.Projectile(12,6)
water = catima.get_material(catima.material.WATER)
water = catima.get_material(catima.material.Water)
water.thickness(10.0)
graphite = catima.get_material(6)
graphite.thickness(1.0)
data = catima.get_data(p, water)
self.assertEqual(catima.max_storage_data,100) # assuming 50, this has to be changed manually
self.assertEqual(catima.max_storage_data,60) # assuming 60, this has to be changed manually
r = catima.storage_info()
#self.assertAlmostEqual(catima.da2de(p,water,et[100]),data[2][100],6)

View File

@ -131,9 +131,9 @@ const lest::test specification[] =
{15.9994,8,1}
});
EXPECT( catima::dedx(p,1000, water) == approx(2.23).R(5e-3));
EXPECT( catima::dedx(p,500, water) == approx(2.76).R(5e-3));
EXPECT( catima::dedx(p,9, water) == approx(51.17).R(5e-3));
EXPECT( catima::dedx(p(1000), water) == approx(2.23).R(5e-3));
EXPECT( catima::dedx(p(500), water) == approx(2.76).R(5e-3));
EXPECT( catima::dedx(p(9), water) == approx(51.17).R(5e-3));
},
CASE("dEdx from spline vs dEdx"){
catima::Projectile p{238,92,92,1000};
@ -142,13 +142,13 @@ const lest::test specification[] =
});
auto res = catima::calculate(p(1000),graphite);
EXPECT(catima::dedx(p,1000, graphite) == approx(res.dEdxi).R(0.001) );
EXPECT(catima::dedx(p(1000), graphite) == approx(res.dEdxi).R(0.001) );
res = catima::calculate(p,graphite,500);
EXPECT(catima::dedx(p,500, graphite) == approx(res.dEdxi).R(0.001) );
EXPECT(catima::dedx(p(500), graphite) == approx(res.dEdxi).R(0.001) );
res = catima::calculate(p,graphite,9);
EXPECT(catima::dedx(p,9, graphite) == approx(res.dEdxi).R(0.001) );
EXPECT(catima::dedx(p(9), graphite) == approx(res.dEdxi).R(0.001) );
},
// CASE("dOmega2dx Ferisov test"){
@ -256,10 +256,10 @@ const lest::test specification[] =
auto re = catima::energy_straggling_from_E(p,res1.Ein,res1.Eout,graphite);
EXPECT(res1.sigma_E == re);
auto eo1 = catima::energy_out(p,1000,graphite);
auto eo1 = catima::energy_out(p(1000),graphite);
EXPECT(res1.Eout == eo1);
auto de1 = catima::dedx_from_range(p,1000,graphite);
auto de1 = catima::dedx_from_range(p(1000),graphite);
EXPECT(res1.dEdxi == de1);
},
@ -382,15 +382,15 @@ const lest::test specification[] =
auto res2 = catima::energy_out(p,energies, graphite);
EXPECT(res2.size()==energies.size());
EXPECT(res2[2] == approx(997.07,01));
EXPECT(res2[0] == approx(catima::energy_out(p,energies[0],graphite),0.1));
EXPECT(res2[1] == approx(catima::energy_out(p,energies[1],graphite),0.1));
EXPECT(res2[2] == approx(catima::energy_out(p,energies[2],graphite),0.1));
EXPECT(res2[0] == approx(catima::energy_out(p(energies[0]),graphite),0.1));
EXPECT(res2[1] == approx(catima::energy_out(p(energies[1]),graphite),0.1));
EXPECT(res2[2] == approx(catima::energy_out(p(energies[2]),graphite),0.1));
auto res3 = catima::dedx_from_range(p,energies,graphite);
EXPECT(res3.size()==energies.size());
EXPECT(res3[0] == approx(catima::dedx_from_range(p,energies[0],graphite),0.1));
EXPECT(res3[1] == approx(catima::dedx_from_range(p,energies[1],graphite),0.1));
EXPECT(res3[2] == approx(catima::dedx_from_range(p,energies[2],graphite),0.1));
EXPECT(res3[0] == approx(catima::dedx_from_range(p(energies[0]),graphite),0.1));
EXPECT(res3[1] == approx(catima::dedx_from_range(p(energies[1]),graphite),0.1));
EXPECT(res3[2] == approx(catima::dedx_from_range(p(energies[2]),graphite),0.1));
},
CASE("constants"){
using namespace catima;

View File

@ -58,7 +58,7 @@ void comp_dedx(catima::Projectile p, catima::Material t, double epsilon = 0.001,
double v1, v2;
for(auto _e : etable.values){
p.T = _e;
v1 = catima::dedx(p,_e,t);
v1 = catima::dedx(p(_e),t);
auto ares = catima::calculate(p(_e),t);
v2 = ares.dEdxi;
dif = rdif(v1,v2);

View File

@ -200,6 +200,7 @@ const lest::test specification[] =
catima::Material mat6;
mat6 = mat5;
EXPECT(mat5==mat6);
EXPECT(mat5.density() == mat6.density());
EXPECT(mat5.ncomponents()==mat6.ncomponents());
EXPECT(mat5.get_element(0).A==mat6.get_element(0).A);
EXPECT(mat5.get_element(1).A==mat6.get_element(1).A);