diff --git a/CMakeLists.txt b/CMakeLists.txt
index 86f506c..5395dab 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,27 +1,30 @@
-cmake_minimum_required(VERSION 3.2.0)
+ cmake_minimum_required(VERSION 3.2.0)
 project(catima)
 
 ############ options #############
-#option(THREADS "Use multi-threading" ON)
+option(BUILD_SHARED_LIBS "build as shared library" ON)
 option(PYTHON_MODULE "compile the Catima python module(requires numpy and cython installed)" OFF)
 option(TESTS "build tests" OFF)
 option(EXAMPLES "build examples" ON)
-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(APPS "build catima applications" ON)
 option(GLOBAL "build with global, sources are required" OFF)
 option(REACTIONS "enable/disable nuclear reaction rate" ON)
-option(APPS "build catima applications" ON)
+option(STORE_SPLINES "store splines, if disables splines are always recreated" ON)
+option(GSL_INTEGRATION "use GSL integration" OFF)
+option(GSL_INTERPOLATION "use GSL inteRPOLATION" OFF)
+option(THIN_TARGET_APPROXIMATION "thin target approximation" ON)
+option(GENERATE_DATA "make data tables generator" OFF)
+option(DOCS "build documentation (requires doxygen)" OFF)
 
 ######## build type ############
 if(NOT CMAKE_BUILD_TYPE  STREQUAL "Debug")
     set(CMAKE_BUILD_TYPE "Release")
+    set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE}")
     MESSAGE(STATUS "Build type Release")
 else()
     set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -g")
     if (${CMAKE_CXX_COMPILER_ID} STREQUAL "Clang" OR ${CMAKE_CXX_COMPILER_ID} STREQUAL "GNU")
-        set(CMAKE_CXX_FLAGS_DEBUG "-Wall -Wextra -Wfatal-errors -g")
+        set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -Wall -Wextra -Wfatal-errors -Wno-unused-parameter -Wno-sign-compare")
     endif()
     MESSAGE(STATUS "Build type Debug")
 endif()
@@ -41,18 +44,10 @@ else()
 endif()
 
 ############# Requirements ##################
-find_package(GSL REQUIRED)
-MESSAGE(STATUS "GSL include dirs: " ${GSL_INCLUDE_DIRS})
-#if(THREADS)
-#    find_package(Threads REQUIRED)
-#    set (EXTRA_LIBS ${EXTRA_LIBS} ${CMAKE_THREAD_LIBS_INIT})
-#    set (USE_THREADS ON)
-#    MESSAGE(STATUS "Nurex will use threads")
-#endif(THREADS)
-
-find_package(PythonInterp)
-if(PYTHONINTERP_FOUND)
-    message("-- Python found: ${PYTHON_EXECUTABLE}")
+if(GSL_INTEGRATION OR GSL_INTERPOLATION)
+    find_package(GSL REQUIRED)
+    MESSAGE(STATUS "GSL include dirs: " ${GSL_INCLUDE_DIRS})
+    list(APPEND EXTRA_LIBS ${GSL_LIBRARIES} )
 endif()
 
 find_package(nurex QUIET)
@@ -81,29 +76,19 @@ endif(GLOBAL)
 file(GLOB HEADERS *.h)
 
 add_library(catima SHARED ${SOURCES})
-add_library(catima_static STATIC ${SOURCES})
 set_target_properties(catima PROPERTIES
 		      POSITION_INDEPENDENT_CODE ON
                       LIBRARY_OUTPUT_DIRECTORY ${PROJECT_BINARY_DIR}/lib
                       )
-set_target_properties(catima_static 
-                        PROPERTIES OUTPUT_NAME catima 
-                        POSITION_INDEPENDENT_CODE ON
-                        ARCHIVE_OUTPUT_DIRECTORY ${PROJECT_BINARY_DIR}/lib
-                    )
-target_link_libraries(catima ${EXTRA_LIBS} ${GSL_LIBRARIES})
-target_link_libraries(catima_static ${EXTRA_LIBS} ${GSL_LIBRARIES})
+
+target_link_libraries(catima ${EXTRA_LIBS})
 
 target_include_directories(catima 
                            PUBLIC  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
                                    $<BUILD_INTERFACE:${GSL_INCLUDE_DIRS}>
                           )
-target_include_directories(catima_static
-                           PUBLIC  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
-                                   $<BUILD_INTERFACE:${GSL_INCLUDE_DIRS}>
-                          )
+
 add_library(catima::catima ALIAS catima)
-add_library(catima::catima_static ALIAS catima_static)
 
 FILE(COPY ${HEADERS} DESTINATION ${PROJECT_BINARY_DIR}/include/catima)
 
@@ -111,10 +96,23 @@ FILE(COPY ${HEADERS} DESTINATION ${PROJECT_BINARY_DIR}/include/catima)
 MESSAGE( STATUS "CMAKE_CXX_COMPILER: " ${CMAKE_CXX_COMPILER} )
 
 ######## for python module
+find_package(PythonInterp)
+if(PYTHONINTERP_FOUND)
+    message("-- Python found: ${PYTHON_EXECUTABLE}")
+endif()
 if(PYTHON_MODULE)
     if(NOT PYTHONINTERP_FOUND)
         MESSAGE(SEND_ERROR "Python is required to build nurex python modules") 
     endif(NOT PYTHONINTERP_FOUND)
+    #find_package(pybind11 REQUIRED)
+    #set(PYBIND11_CPP_STANDARD -std=c++14)
+    #pybind11_add_module(pycatima pymodule/pycatima)
+    #target_include_directories(pycatima PUBLIC
+    #                                $<BUILD_INTERFACE:${CMAKE_CURRENT_BINARY_DIR}/include>
+    #                                $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/libs>
+    #                                $<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"
@@ -182,7 +180,7 @@ endif(APPS)
 ####### install part #######
 FILE(GLOB headers "*.h")
 include(GNUInstallDirs)
-install (TARGETS catima catima_static 
+install (TARGETS catima
         EXPORT catimaConfig
         LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}
         ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR})
@@ -194,7 +192,7 @@ install(EXPORT catimaConfig
         DESTINATION lib/cmake/catima
          )
          
-export(TARGETS catima catima_static NAMESPACE catima:: FILE catimaConfig.cmake)
+export(TARGETS catima NAMESPACE catima:: FILE catimaConfig.cmake)
 export(PACKAGE catima)
 
 ###### packaging #######
diff --git a/README.md b/README.md
index 1e08fd4..604bda5 100644
--- a/README.md
+++ b/README.md
@@ -26,17 +26,18 @@ compile options, enable or disable with cmake:
 > cmake ../ -D[OPTION]
 
 available options:
+  * BUILD_SHARED_LIBS - if ON shared library is build, otherwise static
   * PYTHON_MODULE - enable/disable building of the python bindigs, cython and numpy are required to build the catima python module, default OFF
   * TESTS - build tests
   * EXAMPLES - build examples
   * DOCS - prepare doxygen documentation (after cmake, __make docs__ needs to be executed)
   * GENERATE_DATA - makes program to re-generate precalculated tables (ie precalculated LS coefficients), default:OFF
   * THIN_TARGET_APPROXIMATION - compile the library with thin target approximation, default: ON
-  * GSL_INTEGRATION - use GSL integration functions, otherwise use built-in integrator, default: ON
+  * GSL_INTEGRATION - use GSL integration functions, otherwise use built-in integrator, default: OFF
   * GLOBAL - compile with GLOBAL code (source not included at the moment, needs to be manually added to __global__ directory, default:OFF)
 
 ie:
-> cmake -DCATIMA_PYTHON=ON -DEXAMPLES=ON ../
+> cmake -DPYTHON_MODULE=ON -DEXAMPLES=ON ../
 
 
 after the compilation the libraries and headers must be either installed system-wide by make install or PATH and LD_LIBRARY_PATH must be adjusted to point to headers and library files.
diff --git a/bin/catima_calculator.cpp b/bin/catima_calculator.cpp
index ec27df6..fff7a69 100644
--- a/bin/catima_calculator.cpp
+++ b/bin/catima_calculator.cpp
@@ -18,6 +18,19 @@ void help(){
         std::cout<<"usage: catima_calculator config_file.json\n";
 }
 
+inline std::vector<double> linspace_vector(double a, double b, unsigned int num){
+    std::vector<double> res;
+    if(num>=2 && a<b){
+    res.resize(num);
+    double step = (b-a)/(num-1);
+    for(unsigned int i=0;i<(num-1);i++){
+        res[i]=a+(i*step);
+        }
+    res[num-1] = b;
+    }
+    return res;
+    }
+
 json load_json(const char *fname);
 char* getCmdOption(char ** begin, char ** end, const std::string & option);
 Material json_material(json &j);
diff --git a/build_config.in b/build_config.in
index b770cb4..5d03161 100644
--- a/build_config.in
+++ b/build_config.in
@@ -3,6 +3,8 @@
 
 #cmakedefine THIN_TARGET_APPROXIMATION
 #cmakedefine GSL_INTEGRATION
+#cmakedefine GSL_INTERPOLATION
+#cmakedefine STORE_SPLINES
 #cmakedefine GLOBAL
 #cmakedefine REACTIONS
 #cmakedefine NUREX
diff --git a/calculations.cpp b/calculations.cpp
index 001d5e6..d8d5b73 100644
--- a/calculations.cpp
+++ b/calculations.cpp
@@ -43,11 +43,17 @@ double dedx_n(const Projectile &p, const Target &t){
     double asum = p.A + t.A;
     double epsilon = 32.53*t.A*1000*p.T*p.A/(p.Z*t.Z*asum*zpowers); //projectile energy is converted from MeV/u to keV
     double sn=0;
-    
-    if(epsilon<=30){
+    if(epsilon<=0){
+        return 0.0;
+    }
+    else if(epsilon<=30){
+        assert(p.A>0);
+        assert(epsilon>0);
         sn = log(1+(1.1383*epsilon))/ (2*(epsilon + 0.01321*pow(epsilon,0.21226) + 0.19593*pow(epsilon,0.5)));
     }
     else{
+        assert(p.A>0);
+        assert(epsilon>0);
         sn = log(epsilon)/(2*epsilon);
     }
     sn = 100*8.4621*p.Z*t.Z*p.A*sn*Avogadro/(asum*zpowers*t.A);
@@ -798,6 +804,7 @@ double z_eff_global(double pz, double E, double tz){
         #ifdef GLOBAL
         return global_qmean(pz, tz, E);
         #else
+        assert(false);
         return -1;
         #endif
 }
@@ -871,6 +878,8 @@ double z_eff_atima14(double pz, double T, double tz){
                 }
         }
     }
+    #else
+    assert(false);
     #endif
     return qmean;
     }
diff --git a/calculations.h b/calculations.h
index 119ea2b..2f29511 100644
--- a/calculations.h
+++ b/calculations.h
@@ -39,9 +39,13 @@ namespace catima{
     double reduced_energy_loss_unit(const Projectile &p, const Target &t);
     
 
-
+    /**
+     * @brief bethek_dedx_e - electronics stopping power
+     * @return stopping power
+     */
     double bethek_dedx_e(Projectile &p,const Target &t, const Config &c=default_config, double I=0.0);
     double bethek_dedx_e(Projectile &p,const Material &mat, const Config &c=default_config);
+
     double bethek_barkas(double zp_eff,double eta, double zt);
     double bethek_density_effect(double beta, int zt);
     
diff --git a/catima.cpp b/catima.cpp
index 5032a03..4eefc7a 100644
--- a/catima.cpp
+++ b/catima.cpp
@@ -17,7 +17,6 @@ namespace catima{
 
 Config default_config;
 
-
 bool operator==(const Config &a, const Config&b){
     if(std::memcmp(&a,&b,sizeof(Config)) == 0){
         return true;
@@ -30,11 +29,9 @@ bool operator==(const Config &a, const Config&b){
 double dedx(Projectile &p, double T, const Material &mat, const  Config &c){
     double sum = 0;
     if(T<=0)return 0.0;
-
-    sum += dedx_n(p,mat);
-
-    double se=0;
     p.T = T;
+    sum += dedx_n(p,mat);
+    double se=0;
     if(p.T<=10){
         se = sezi_dedx_e(p,mat);
     }
@@ -78,21 +75,23 @@ double da2dx(Projectile &p, double T, const Material &mat, const Config &c){
 
 
 double range(Projectile &p, double T, 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);
+    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);
 }
 
 double dedx_from_range(Projectile &p, double T, 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);
+    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);
 }
 
 std::vector<double> dedx_from_range(Projectile &p, const std::vector<double> &T, 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);
-
+    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);
     std::vector<double> dedx;
     dedx.reserve(T.size());
     for(auto e:T){
@@ -107,45 +106,52 @@ std::vector<double> dedx_from_range(Projectile &p, const std::vector<double> &T,
 }
 
 double range_straggling(Projectile &p, double T, const Material &t, const Config &c){
-    auto data = _storage.Get(p,t,c);
-    Interpolator range_straggling_spline(energy_table.values,data.range_straggling.data(),energy_table.num);
+    auto& data = _storage.Get(p,t,c);
+    //Interpolator range_straggling_spline(energy_table.values,data.range_straggling.data(),energy_table.num);
+    spline_type range_straggling_spline = get_range_straggling_spline(data);
     return sqrt(range_straggling_spline(T));
 }
 
 double range_variance(Projectile &p, double T, const Material &t, const Config &c){
-    auto data = _storage.Get(p,t,c);
-    Interpolator range_straggling_spline(energy_table.values,data.range_straggling.data(),energy_table.num);
+    auto& data = _storage.Get(p,t,c);
+    //Interpolator range_straggling_spline(energy_table.values,data.range_straggling.data(),energy_table.num);
+    spline_type range_straggling_spline = get_range_straggling_spline(data);
     return range_straggling_spline(T);
 }
 
 double domega2de(Projectile &p, double T, const Material &t, const Config &c){
-    auto data = _storage.Get(p,t,c);
-    Interpolator range_straggling_spline(energy_table.values,data.range_straggling.data(),energy_table.num);
+    auto& data = _storage.Get(p,t,c);
+    //Interpolator range_straggling_spline(energy_table.values,data.range_straggling.data(),energy_table.num);
+    spline_type range_straggling_spline = get_range_straggling_spline(data);
     return range_straggling_spline.derivative(T);
 }
 
 double da2de(Projectile &p, double T, const Material &t, const Config &c){
-    auto data = _storage.Get(p,t,c);
-    Interpolator angular_variance_spline(energy_table.values,data.angular_variance.data(),energy_table.num);
+    auto& data = _storage.Get(p,t,c);
+    //Interpolator angular_variance_spline(energy_table.values,data.angular_variance.data(),energy_table.num);
+    spline_type angular_variance_spline = get_angular_variance_spline(data);
     return angular_variance_spline.derivative(T);
 }
 
 double angular_straggling_from_E(Projectile &p, double T, double Tout, const Material &t, const Config &c){
-    auto data = _storage.Get(p,t,c);
-    Interpolator angular_straggling_spline(energy_table.values,data.angular_variance.data(),energy_table.num);
-    return sqrt(angular_straggling_spline(T) - angular_straggling_spline(Tout));
+    auto& data = _storage.Get(p,t,c);
+    //Interpolator angular_straggling_spline(energy_table.values,data.angular_variance.data(),energy_table.num);
+    spline_type angular_variance_spline = get_angular_variance_spline(data);
+    return sqrt(angular_variance_spline(T) - angular_variance_spline(Tout));
 }
 
 double energy_straggling_from_E(Projectile &p, double T, double Tout,const Material &t, const Config &c){
-    auto data = _storage.Get(p,t,c);
+    auto& data = _storage.Get(p,t,c);
 
-    Interpolator range_straggling_spline(energy_table.values,data.range_straggling.data(),energy_table.num);
-    Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num);
+    //Interpolator range_straggling_spline(energy_table.values,data.range_straggling.data(),energy_table.num);
+    //Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num);
+    spline_type range_spline = get_range_spline(data);
+    spline_type range_straggling_spline = get_range_straggling_spline(data);
     double dEdxo = p.A/range_spline.derivative(Tout);
     return dEdxo*sqrt(range_straggling_spline(T) - range_straggling_spline(Tout))/p.A;
 }
 
-double energy_out(double T, double thickness, Interpolator &range_spline){
+double energy_out(double T, double thickness, const Interpolator &range_spline){
     constexpr double epsilon = 1E-5;
     int counter = 0;
     double range;
@@ -172,14 +178,16 @@ double energy_out(double T, double thickness, Interpolator &range_spline){
 }
 
 double energy_out(Projectile &p, double T, 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);
+    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);
     }
 
 std::vector<double> energy_out(Projectile &p, const std::vector<double> &T, 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);
+    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);
 
     std::vector<double> eout;
     eout.reserve(T.size());
@@ -199,15 +207,18 @@ Result calculate(Projectile &p, const Material &t, const Config &c){
     Result res;
     double T = p.T;
     if(T<catima::Ezero && T<catima::Ezero-catima::numeric_epsilon){return res;}
-    auto data = _storage.Get(p,t,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);
 
-    Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num);
     res.Ein = T;
     res.range = range_spline(T);
     res.dEdxi = p.A/range_spline.derivative(T);
     res.Eout = energy_out(T,t.thickness(),range_spline);
 
-    Interpolator range_straggling_spline(energy_table.values,data.range_straggling.data(),energy_table.num);
+    //Interpolator range_straggling_spline(energy_table.values,data.range_straggling.data(),energy_table.num);
+    spline_type range_straggling_spline = get_range_straggling_spline(data);
 
     if(res.Eout<Ezero){
         res.dEdxo = 0.0;
@@ -225,20 +236,23 @@ Result calculate(Projectile &p, const Material &t, const Config &c){
             double s2 = range_straggling_spline.derivative(res.Eout);
             res.sigma_E = res.dEdxo*sqrt(edif*0.5*(s1+s2))/p.A;
         
-            Interpolator angular_variance_spline(energy_table.values,data.angular_variance.data(),energy_table.num);
+            //Interpolator angular_variance_spline(energy_table.values,data.angular_variance.data(),energy_table.num);
+            spline_type angular_variance_spline = get_angular_variance_spline(data);
             s1 = angular_variance_spline.derivative(T);
             s2 = angular_variance_spline.derivative(res.Eout);
             res.sigma_a = sqrt(0.5*(s1+s2)*edif);
         }
         else{
             res.sigma_E = res.dEdxo*sqrt(range_straggling_spline(T) - range_straggling_spline(res.Eout))/p.A;
-            Interpolator angular_variance_spline(energy_table.values,data.angular_variance.data(),energy_table.num);
+            //Interpolator angular_variance_spline(energy_table.values,data.angular_variance.data(),energy_table.num);
+            spline_type angular_variance_spline = get_angular_variance_spline(data);
             res.sigma_a = sqrt(angular_variance_spline(T) - angular_variance_spline(res.Eout));    
         }
         
         #else
         res.sigma_E = res.dEdxo*sqrt(range_straggling_spline(T) - range_straggling_spline(res.Eout))/p.A;
-        Interpolator angular_variance_spline(energy_table.values,data.angular_variance.data(),energy_table.num);
+        //Interpolator angular_variance_spline(energy_table.values,data.angular_variance.data(),energy_table.num);
+        spline_type angular_variance_spline = get_angular_variance_spline(data);
         res.sigma_a = sqrt(angular_variance_spline(T) - angular_variance_spline(res.Eout));
         #endif
         if( !(c.skip&skip_tof) && t.thickness()>0){
@@ -371,7 +385,7 @@ std::vector<double> calculate_tof(Projectile p, const Material &t, const Config
     }
     return values;
 }
-
+/*
 DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){
     DataPoint dp(p,t,c);
     dp.range.resize(max_datapoints);
@@ -411,6 +425,47 @@ DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){
     }
     return dp;
 }
+*/
+
+
+DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){
+    DataPoint dp(p,t,c);
+    dp.range.resize(max_datapoints);
+    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);
+            };
+    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);
+            };
+
+    double res=0.0;
+    //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;
+    dp.range[0] = 0.0;
+
+    dp.angular_variance[0] = 0.0;
+
+    //res = integrator.integrate(fomega,Ezero,energy_table(0));
+    //res = p.A*res;
+    dp.range_straggling[0]=0.0;
+    //p.T = energy_table(0);
+    for(int i=1;i<max_datapoints;i++){
+        res = p.A*integrator.integrate(fdedx,energy_table(i-1),energy_table(i));
+        dp.range[i] = res + dp.range[i-1];
+        res = da2dx(p,energy_table(i),t)*res;
+        dp.angular_variance[i] = res + dp.angular_variance[i-1];
+
+        res = integrator.integrate(fomega,energy_table(i-1),energy_table(i));
+        res = p.A*res;
+        dp.range_straggling[i] = res + dp.range_straggling[i-1];
+    }
+    return dp;
+}
 
 double calculate_tof_from_E(Projectile p, double Eout, const Material &t, const Config &c){
     double res;
diff --git a/catima.h b/catima.h
index 12256ca..f841f3b 100644
--- a/catima.h
+++ b/catima.h
@@ -21,6 +21,7 @@
 #include <vector>
 
 // #define NDEBUG
+#include "catima/build_config.h"
 #include "catima/config.h"
 #include "catima/constants.h"
 #include "catima/structures.h"
@@ -141,7 +142,7 @@ namespace catima{
       * @range_spline - precaclulated range spline for material 
       * @return outcoming energy after the thickness in Mev/u
       */
-    double energy_out(double T, double thickness, Interpolator &range_spline);
+    double energy_out(double T, double thickness, const Interpolator &range_spline);
 
     /**
       * calculates outcoming energy 
diff --git a/catimac.pxd b/catimac.pxd
index 56e8011..8b643e6 100644
--- a/catimac.pxd
+++ b/catimac.pxd
@@ -114,7 +114,7 @@ cdef extern from "catima/constants.h" namespace "catima":
     int logEmin "catima::logEmin"
     int logEmax "catima::logEmax"
     bool reactions "catima::reactions"
-    
+
 cdef extern from "catima/storage.h" namespace "catima":
     cdef cppclass Interpolator:
         Interpolator(const double *x, const double *y, int num) except +
@@ -141,4 +141,4 @@ cdef extern from "catima/storage.h" namespace "catima":
         
     cdef EnergyTableType energy_table;
     cdef Data _storage;
-    cdef DataPoint& get_data(const Projectile &p, const Material &t, Config c);
+    cdef DataPoint& get_data(const Projectile &p, const Material &t, const Config c);
diff --git a/integrator.cpp b/integrator.cpp
index b0b7dcb..2c8a176 100644
--- a/integrator.cpp
+++ b/integrator.cpp
@@ -1,12 +1,13 @@
 #include "integrator.h"
+
+//#ifdef GSL_INTEGRATION
 #include "gsl/gsl_integration.h"
 #include "gsl/gsl_errno.h"
+//#endif
 
 namespace catima{
-    
     integrator_type integrator;
-    IntegratorGSL integratorGSL(true);
-    
+#ifdef GSL_INTEGRATION
     double funcwrapper3(double x, void *_c){
     std::function<double(double)> *f = (std::function<double(double)> *)_c;
     return (*f)(x);
@@ -44,6 +45,5 @@ namespace catima{
             }
         return result;
     };
-
-
+#endif
 }
diff --git a/integrator.h b/integrator.h
index 3bad088..56bd398 100644
--- a/integrator.h
+++ b/integrator.h
@@ -17,15 +17,20 @@
 #ifndef INTEGRATOR_H
 #define INTEGRATOR_H
 #include "catima/build_config.h"
-#include "gsl/gsl_integration.h"
 #include <functional>
 #include <array>
-#ifdef USE_THREADS
-#include <mutex>
+
+//#ifdef USE_THREADS
+//#include <mutex>
+//#endif
+
+#ifdef GSL_INTEGRATION
+#include "gsl/gsl_integration.h"
 #endif
 
 namespace catima{
     
+#ifdef GSL_INTEGRATION
 /// helper class to integrate functions using the GSL library
 class IntegratorGSL{
 	public:
@@ -45,6 +50,7 @@ class IntegratorGSL{
 	std::mutex integration_mutex;
 	#endif
     };
+#endif
 
 template<int order>
 struct GL_data{
@@ -176,7 +182,6 @@ using integrator_type = GaussLegendreIntegration<8>;
 #endif
 
 extern integrator_type integrator;
-extern IntegratorGSL integratorGSL;
 }
 
 #endif
diff --git a/reactions.cpp b/reactions.cpp
index 127670f..5812b2f 100644
--- a/reactions.cpp
+++ b/reactions.cpp
@@ -23,8 +23,8 @@ double nonreaction_rate(Projectile &projectile, const Material &target, const Co
     int ap = lround(projectile.A);
     int zp = lround(projectile.Z);
 
-    auto data = _storage.Get(projectile,target,c);
-    Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num);
+    auto& data = _storage.Get(projectile,target,c);
+    spline_type range_spline = get_range_spline(data);
     if(energy_out(projectile.T, target.thickness(), range_spline) < emin_reaction)return -1.0;
     
     auto sigma_r = [&](double th){
diff --git a/setup.py.in b/setup.py.in
index f43777a..6715d8d 100644
--- a/setup.py.in
+++ b/setup.py.in
@@ -10,8 +10,8 @@ setup(
                 library_dirs=["${CMAKE_CURRENT_BINARY_DIR}/lib"],
 		include_dirs=["${CMAKE_CURRENT_BINARY_DIR}/include"],
 #		extra_objects=["${CATIMA_LIB}"],
-		extra_compile_args=["-std=c++14"],
-		extra_link_args=["-std=c++14"]
+                extra_compile_args=["-std=c++14"],
+                extra_link_args=["-std=c++14"]
 		),
 	    ])
 )
diff --git a/spline.cpp b/spline.cpp
index 9ec79ba..76f4cf7 100644
--- a/spline.cpp
+++ b/spline.cpp
@@ -10,247 +10,6 @@
 
 namespace catima{
 
-band_matrix::band_matrix(int dim)
-{
-    resize(dim);
-}
-void band_matrix::resize(int dim)
-{
-    assert(dim>0);
-    a.resize(dim);
-    d.resize(dim);
-    c.resize(dim);
-}
-int band_matrix::dim() const
-{
-    return d.size();
-}
-
-
-// defines the new operator (), so that we can access the elements
-// by A(i,j), index going from i=0,...,dim()-1
-double & band_matrix::operator () (int i, int j)
-{
-    int k=j-i;       // what band is the entry
-    assert( (i>=0) && (i<dim()) && (j>=0) && (j<dim()) );
-    assert(k<2 && k>-2);
-    if(k>0)return c[i];
-    else if(k==0) return d[i];
-    else return a[i];
-}
-double band_matrix::operator () (int i, int j) const
-{
-    int k=j-i;       // what band is the entry
-    assert( (i>=0) && (i<dim()) && (j>=0) && (j<dim()) );
-    if(k>0)return c[i];
-    else if(k==0) return d[i];
-    else return a[i];
-}
-
-
-
-std::vector<double> band_matrix::trig_solve(const std::vector<double>& b) const
-{
-    assert( this->dim()==(int)b.size() );
-    std::vector<double> x(this->dim());
-    std::vector<double> g(this->dim());
-    
-    assert(d[0]!=0.0);
-    x[0] = b[0]/d[0];
-    double bet = d[0];
-    for(int j=1;j<this->dim();j++){
-        g[j] = c[j-1]/bet;
-        bet = d[j] - (a[j]*g[j]);
-        assert(bet != 0.0);
-        x[j] = (b[j]-a[j]*x[j-1])/bet;
-    }
-    for(int j=this->dim()-2;j>=0;j--){
-        x[j] -= g[j+1]*x[j+1];
-    }
-    
-    return x;
-}
-
-
-// spline implementation
-// -----------------------
-
-void spline::set_boundary(spline::bd_type left, double left_value,
-                          spline::bd_type right, double right_value,
-                          bool force_linear_extrapolation)
-{
-    assert(n==0);          // set_points() must not have happened yet
-    m_left=left;
-    m_right=right;
-    m_left_value=left_value;
-    m_right_value=right_value;
-    m_force_linear_extrapolation=force_linear_extrapolation;
-}
-
-
-void spline::set_points(const double *x,
-                        const double *y,
-                        const size_t num
-                        )
-{
-    assert(num>2);
-    m_x=x;
-    m_y=y;
-    n=num;
-    // TODO: maybe sort x and y, rather than returning an error
-    for(int i=0; i<n-1; i++) {
-        assert(m_x[i]<m_x[i+1]);
-    }
-
-    
-        // setting up the matrix and right hand side of the equation system
-        // for the parameters b[]
-        band_matrix A(n);
-        std::vector<double>  rhs(n);
-        for(int i=1; i<n-1; i++) {
-            A(i,i-1)=1.0/3.0*(x[i]-x[i-1]);
-            A(i,i)=2.0/3.0*(x[i+1]-x[i-1]);
-            A(i,i+1)=1.0/3.0*(x[i+1]-x[i]);
-            rhs[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
-        }
-        // boundary conditions
-        if(m_left == spline::bd_type::second_deriv) {
-            // 2*b[0] = f''
-            A(0,0)=2.0;
-            A(0,1)=0.0;
-            rhs[0]=m_left_value;
-        } else{
-            // c[0] = f', needs to be re-expressed in terms of b:
-            // (2b[0]+b[1])(x[1]-x[0]) = 3 ((y[1]-y[0])/(x[1]-x[0]) - f')
-            A(0,0)=2.0*(x[1]-x[0]);
-            A(0,1)=1.0*(x[1]-x[0]);
-            rhs[0]=3.0*((y[1]-y[0])/(x[1]-x[0])-m_left_value);
-        } 
-
-        if(m_right == spline::bd_type::second_deriv) {
-            // 2*b[n-1] = f''
-            A(n-1,n-1)=2.0;
-            A(n-1,n-2)=0.0;
-            rhs[n-1]=m_right_value;
-        } else{
-            // c[n-1] = f', needs to be re-expressed in terms of b:
-            // (b[n-2]+2b[n-1])(x[n-1]-x[n-2])
-            // = 3 (f' - (y[n-1]-y[n-2])/(x[n-1]-x[n-2]))
-            A(n-1,n-1)=2.0*(x[n-1]-x[n-2]);
-            A(n-1,n-2)=1.0*(x[n-1]-x[n-2]);
-            rhs[n-1]=3.0*(m_right_value-(y[n-1]-y[n-2])/(x[n-1]-x[n-2]));
-        }
-
-        // solve the equation system to obtain the parameters b[]
-        //m_b=A.lu_solve(rhs);
-        m_b=A.trig_solve(rhs);
-
-        // calculate parameters a[] and c[] based on b[]
-        m_a.resize(n);
-        m_c.resize(n);
-        for(int i=0; i<n-1; i++) {
-            m_a[i]=1.0/3.0*(m_b[i+1]-m_b[i])/(x[i+1]-x[i]);
-            m_c[i]=(y[i+1]-y[i])/(x[i+1]-x[i])
-                   - 1.0/3.0*(2.0*m_b[i]+m_b[i+1])*(x[i+1]-x[i]);
-        }
-    
-    // for left extrapolation coefficients
-    m_b0 = (m_force_linear_extrapolation==false) ? m_b[0] : 0.0;
-    m_c0 = m_c[0];
-
-    // for the right extrapolation coefficients
-    // f_{n-1}(x) = b*(x-x_{n-1})^2 + c*(x-x_{n-1}) + y_{n-1}
-    double h=x[n-1]-x[n-2];
-    // m_b[n-1] is determined by the boundary condition
-    m_a[n-1]=0.0;
-    m_c[n-1]=3.0*m_a[n-2]*h*h+2.0*m_b[n-2]*h+m_c[n-2];   // = f'_{n-2}(x_{n-1})
-    if(m_force_linear_extrapolation==true)
-        m_b[n-1]=0.0;
-}
-
-double spline::operator() (double x) const
-{
-    assert(n>2);
-    // find the closest point m_x[idx] < x, idx=0 even if x<m_x[0]
-    auto it=std::lower_bound(m_x,m_x+n,x);
-    //int idx=std::max( int(it-m_x)-1, 0);
-    if(it!=m_x)it--;
-    int idx = std::distance(m_x,it);
-    double mx = *it;
-    double h=x-mx;
-    double interpol;
-    if(x<m_x[0]) {
-        // extrapolation to the left
-        interpol=(m_b0*h + m_c0)*h + m_y[0];
-    } else if(x>m_x[n-1]) {
-        // extrapolation to the right
-        interpol=(m_b[n-1]*h + m_c[n-1])*h + m_y[n-1];
-    } else {
-        // interpolation
-        interpol=((m_a[idx]*h + m_b[idx])*h + m_c[idx])*h + m_y[idx];
-    }
-    return interpol;
-}
-
-double spline::deriv(int order, double x) const
-{
-    assert(order>0);
-    // find the closest point m_x[idx] < x, idx=0 even if x<m_x[0]
-    auto it=std::lower_bound(m_x,m_x+n,x);
-    //int idx=std::max( int(it-m_x)-1, 0);
-    if(it!=m_x)it--;
-    int idx = std::distance(m_x,it);
-    double mx = *it;
-    double h=x-mx;
-    double interpol;
-    if(x<m_x[0]) {
-        // extrapolation to the left
-        switch(order) {
-        case 1:
-            interpol=2.0*m_b0*h + m_c0;
-            break;
-        case 2:
-            interpol=2.0*m_b0*h;
-            break;
-        default:
-            interpol=0.0;
-            break;
-        }
-    } else if(x>m_x[n-1]) {
-        // extrapolation to the right
-        switch(order) {
-        case 1:
-            interpol=2.0*m_b[n-1]*h + m_c[n-1];
-            break;
-        case 2:
-            interpol=2.0*m_b[n-1];
-            break;
-        default:
-            interpol=0.0;
-            break;
-        }
-    } else {
-        // interpolation
-        switch(order) {
-        case 1:
-            interpol=(3.0*m_a[idx]*h + 2.0*m_b[idx])*h + m_c[idx];
-            break;
-        case 2:
-            interpol=6.0*m_a[idx]*h + 2.0*m_b[idx];
-            break;
-        case 3:
-            interpol=6.0*m_a[idx];
-            break;
-        default:
-            interpol=0.0;
-            break;
-        }
-    }
-    return interpol;
-}
-
-
-
 } // namespace tk
 
 
diff --git a/spline.h b/spline.h
index 1623575..2abea56 100644
--- a/spline.h
+++ b/spline.h
@@ -17,77 +17,208 @@
  *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-#ifndef TK_SPLINE_H
-#define TK_SPLINE_H
+#ifndef CATIMA_SPLINE_H
+#define CATIMA_SPLINE_H
 
 #include <cstdio>
 #include <cassert>
 #include <vector>
 #include <algorithm>
+#include <array>
+#include "catima/constants.h"
+
+#ifdef GSL_INTERPOLATION
+#include <gsl/gsl_spline.h>
+#endif
 
 namespace catima
 {
 
-// band matrix solver
-class band_matrix
+enum interpolation_t {cspline, linear};
+
+/**
+ * Tridiagonal matrix solver
+ */
+
+template<int N>
+class tridiagonal_matrix
 {
 private:
-    std::vector<double> a;
-    std::vector<double> d;
-    std::vector<double> c;
-    std::vector<double> save;
+    std::array<double,N> a;
+    std::array<double,N> d;
+    std::array<double,N> c;
 public:
-    band_matrix() {};                             // constructor
-    band_matrix(int dim);       // constructor
-    ~band_matrix() {};                            // destructor
-    void resize(int dim);      // init with dim,n_u,n_l
-    int dim() const;                             // matrix dimension
+    tridiagonal_matrix() {}
 
     // access operator
-    double & operator () (int i, int j);            // write
-    double   operator () (int i, int j) const;      // read
-    // we can store an additional diogonal (in m_lower)
-    std::vector<double> trig_solve(const std::vector<double>& b) const;
+    double & operator () (unsigned int i, unsigned int j);            // write
+    double   operator () (unsigned int i, unsigned int j) const;      // read
+    std::array<double, N> trig_solve(const std::array<double, N>& b) const;
 };
 
-
-// spline interpolation
-class spline
+template<int N>
+double & tridiagonal_matrix<N>::operator () (unsigned int i, unsigned int j)
 {
-public:
-    enum class bd_type {
-        first_deriv = 1,
-        second_deriv = 2
-    };
+    int k=j-i;
+    if(k == -1)return c[i];
+    else if(k==0) return d[i];
+    else return a[i];
+}
 
-private:
-    const double *m_x, *m_y;            // x,y coordinates of points
-    size_t n=0;
-    // interpolation parameters
-    // f(x) = a*(x-x_i)^3 + b*(x-x_i)^2 + c*(x-x_i) + y_i
-    std::vector<double> m_a,m_b,m_c;        // spline coefficients
-    double  m_b0, m_c0;                     // for left extrapol
-    bd_type m_left = bd_type::second_deriv;
-    bd_type m_right = bd_type::second_deriv;
-    double  m_left_value = 0.0;
-    double  m_right_value = 0.0;
-    bool    m_force_linear_extrapolation = false;
+template<int N>
+double tridiagonal_matrix<N>::operator () (unsigned int i, unsigned int j) const
+{
+    int k=j-i;
+    if(k==-1)return c[i];
+    else if(k==0) return d[i];
+    else if(k==1)return a[i];
+    else return 0.0;
+}
 
-public:
-    // set default boundary condition to be zero curvature at both ends
-    spline(){}
+template<int N>
+std::array<double, N> tridiagonal_matrix<N>::trig_solve(const std::array<double, N>& b) const
+{
+    std::array<double, N> x;
+    if(d[0] == 0.0){return x;}
+    std::array<double, N> g;
+    x[0] = b[0]/d[0];
+    double bet = d[0];
+    for(std::size_t j=1, max=N;j<max;j++){
+        g[j] = c[j-1]/bet;
+        bet = d[j] - (a[j]*g[j]);
+        if(bet == 0.0){
+            x.fill(0.0);
+            return x;
+        }
+        x[j] = (b[j]-a[j]*x[j-1])/bet;
+    }
+    for(int j=N-2;j>=0;j--){
+        x[j] -= g[j+1]*x[j+1];
+    }
+    return x;
+}
 
-    // optional, but if called it has to come be before set_points()
-    void set_boundary(bd_type left, double left_value,
-                      bd_type right, double right_value,
-                      bool force_linear_extrapolation=false);
-    void set_points(const double *x,
-                    const double *y,
-                    const size_t num);
-    double operator() (double x) const;
-    double deriv(int order, double x) const;
+
+/**
+ * Cubic Spline class, accepting EnergyTable type as x-variable
+ */
+template<typename T>
+struct cspline_special{
+    constexpr static int N = T::size();
+    cspline_special(const T& x,
+                    const std::vector<double>& y,
+                    bool boundary_second_deriv = true);
+    cspline_special() = default;
+
+    const T *table;
+    const double *m_x;
+    const double *m_y;
+    std::array<double,N> m_a,m_b,m_c;
+    double  m_b0, m_c0;
+
+
+    double operator()(double x)const{return evaluate(x);}
+    double evaluate(double x) const
+    {
+        int idx=std::max( table->index(x), 0);
+        double h=x-m_x[idx];
+        double interpol;
+        if(x<m_x[0]) {
+            // extrapolation to the left
+            interpol=(m_b0*h + m_c0)*h + m_y[0];
+        } else if(x>m_x[N-1]) {
+            // extrapolation to the right
+            interpol=(m_b[N-1]*h + m_c[N-1])*h + m_y[N-1];
+        } else {
+            // interpolation
+            interpol=((m_a[idx]*h + m_b[idx])*h + m_c[idx])*h + m_y[idx];
+        }
+        return interpol;
+    }
+
+
+    double deriv(double x) const
+    {
+        int idx=std::max( table->index(x), 0);
+
+        double h=x-m_x[idx];
+        double interpol;
+        if(x<m_x[0]) {
+            // extrapolation to the left
+            interpol=2.0*m_b0*h + m_c0;
+        } else if(x>m_x[N-1]) {
+            // extrapolation to the right
+            interpol=2.0*m_b[N-1]*h + m_c[N-1];
+        } else {
+            // interpolation
+            interpol=(3.0*m_a[idx]*h + 2.0*m_b[idx])*h + m_c[idx];
+        }
+        return interpol;
+    }
+    static_assert (T::size()>2, "N must be > 2");
 };
 
-} // namespace tk
+template<typename T>
+cspline_special<T>::cspline_special(const T &x,
+                      const std::vector<double>& y,
+                      bool boundary_second_deriv
+                      ):table(&x),m_y(y.data()),m_x(x.values)
+{
+    static_assert (N>2, "N must be > 2");
+    tridiagonal_matrix<N> A{};
+    std::array<double, N> rhs;
+    for(std::size_t i=1; i<N-1; i++) {
+        A(i,i-1)=1.0/3.0*(x[i]-x[i-1]);
+        A(i,i)=2.0/3.0*(x[i+1]-x[i-1]);
+        A(i,i+1)=1.0/3.0*(x[i+1]-x[i]);
+        rhs[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
+    }
+        // boundary conditions
+    if(boundary_second_deriv) {
+        // 2*b[0] = f''
+        A(0,0)=2.0;
+        A(0,1)=0.0;
+        rhs[0]=0.0; // 0.0 is value of derivative
+        A(N-1,N-1)=2.0;
+        A(N-1,N-2)=0.0;
+        rhs[N-1]=0.0; // 0.0 is value of derivative
+       } else {
+        // c[0] = f', needs to be re-expressed in terms of b:
+        // (2b[0]+b[1])(x[1]-x[0]) = 3 ((y[1]-y[0])/(x[1]-x[0]) - f')
+        A(0,0)=2.0*(x[1]-x[0]);
+        A(0,1)=1.0*(x[1]-x[0]);
+        rhs[0]=3.0*((y[1]-y[0])/(x[1]-x[0])-0.0); // 0.0 is deriv value
 
-#endif /* TK_SPLINE_H */
+        // c[n-1] = f', needs to be re-expressed in terms of b:
+        // (b[n-2]+2b[n-1])(x[n-1]-x[n-2])
+        // = 3 (f' - (y[n-1]-y[n-2])/(x[n-1]-x[n-2]))
+        A(N-1,N-1)=2.0*(x[N-1]-x[N-2]);
+        A(N-1,N-2)=1.0*(x[N-1]-x[N-2]);
+        rhs[N-1]=3.0*(0.0-(y[N-1]-y[N-2])/(x[N-1]-x[N-2]));
+    }
+
+
+        m_b=A.trig_solve(rhs);
+
+        // calculate parameters a[] and c[] based on b[]
+        for(int i=0; i<N-1; i++) {
+            m_a[i]=1.0/3.0*(m_b[i+1]-m_b[i])/(x[i+1]-x[i]);
+            m_c[i]=(y[i+1]-y[i])/(x[i+1]-x[i])
+                   - 1.0/3.0*(2.0*m_b[i]+m_b[i+1])*(x[i+1]-x[i]);
+        }
+
+
+    // for left extrapolation coefficients
+    //s.m_b0 = (m_force_linear_extrapolation==false) ? s.m_b[0] : 0.0;
+    m_b0 =  0.0;
+    m_c0 = m_c[0];
+
+    double h=x[N-1]-x[N-2];
+    m_a[N-1]=0.0;
+    m_c[N-1]=3.0*m_a[N-2]*h*h+2.0*m_b[N-2]*h+m_c[N-2];   // = f'_{n-2}(x_{n-1})
+    m_b[N-1]=0.0;
+}
+
+} // namespace end
+
+#endif
diff --git a/storage.cpp b/storage.cpp
index f3bdafb..a84d033 100644
--- a/storage.cpp
+++ b/storage.cpp
@@ -29,11 +29,72 @@ namespace catima {
 	    return false;
         }
     }
-    
-    DataPoint::~DataPoint(){
-	
+
+#ifdef GSL_INTERPOLATION
+//////////// Interpolator ////////////////////////////////
+InterpolatorGSL::InterpolatorGSL(const EnergyTable<max_datapoints>& x, const std::vector<double>& y, interpolation_t type){
+    acc = gsl_interp_accel_alloc ();
+    const int num = y.size();
+    if(type==cspline)
+    spline = gsl_spline_alloc (gsl_interp_cspline, num);
+    else
+    spline = gsl_spline_alloc (gsl_interp_linear, num);
+
+    gsl_spline_init (spline, x.values, y.data(), num);
+    min= x[0];
+    max= x[num-1];
+
+}
+
+InterpolatorGSL::~InterpolatorGSL(){
+    gsl_interp_accel_free (acc);
+    gsl_spline_free (spline);
+}
+
+double InterpolatorGSL::eval(double x) const{
+    if(x<min)x=min;
+    if(x>max)x=max;
+    return gsl_spline_eval(spline, x, acc);
+}
+
+double InterpolatorGSL::derivative(double x)const{
+    if(x<min)x=min;
+    if(x>max)x=max;
+    return gsl_spline_eval_deriv (spline, x, acc);
+}
+#endif
+
+#ifdef STORE_SPLINES
+    const Interpolator& get_range_spline(const DataPoint &data){
+        return data.range_spline;
     }
-    
+
+    const Interpolator& get_range_straggling_spline(const DataPoint &data){
+        return data.range_straggling_spline;
+    }
+
+    const Interpolator& get_angular_variance_spline(const DataPoint &data){
+        return data.angular_variance_spline;
+    }
+#else
+    Interpolator get_range_spline(const DataPoint &data){
+        //return Interpolator(energy_table.values,data.range);
+        //return data.range_spline;
+        return Interpolator(energy_table,data.range);
+    }
+
+    Interpolator get_range_straggling_spline(const DataPoint &data){
+        //return Interpolator(energy_table.values,data.range_straggling);
+        //return data.range_straggling_spline;
+        return Interpolator(energy_table,data.range_straggling);
+    }
+
+    Interpolator get_angular_variance_spline(const DataPoint &data){
+        //return Interpolator(energy_table.values,data.angular_variance);
+        //return data.angular_variance_spline;
+        return Interpolator(energy_table,data.angular_variance);
+    }
+#endif
     Data::Data(){
         //storage.reserve(max_storage_data); // disabled because of "circular" storage
         storage.resize(max_storage_data);
@@ -57,15 +118,23 @@ void Data::Add(const Projectile &p, const Material &t, const Config &c){
     index->angular_variance = calculate_angular_variance(p,t,c);
     #else
     *index = calculate_DataPoint(p,t,c);
-    #endif
+#ifdef STORE_SPLINES
+    //index->range_spline = Interpolator(energy_table.values,index->range);
+    //index->range_straggling_spline = Interpolator(energy_table.values,index->range_straggling);
+    //index->angular_variance_spline = Interpolator(energy_table.values,index->angular_variance);
+    index->range_spline = Interpolator(energy_table, index->range);
+    index->range_straggling_spline = Interpolator(energy_table, index->range_straggling);
+    index->angular_variance_spline = Interpolator(energy_table, index->angular_variance);
+#endif
+#endif
 
     index++;
     }
     
-DataPoint& Data::Get(const Projectile &p, const Material &t, const Config &c){
+ DataPoint& Data::Get(const Projectile &p, const Material &t, const Config &c){
 	for(auto &e:storage){
 	    if( (e.p==p) && (e.m==t) && (e.config==c)){
-		return e;
+        return e;
 	    }
 	}
     Add(p,t,c);
@@ -73,71 +142,4 @@ DataPoint& Data::Get(const Projectile &p, const Material &t, const Config &c){
     return *std::prev(index);
     }
 
-//////////// Interpolator ////////////////////////////////
-InterpolatorGSL::InterpolatorGSL(const double *x, const double *y, int num,interpolation_t type){
-    acc = gsl_interp_accel_alloc ();
-
-    if(type==cspline)
-	spline = gsl_spline_alloc (gsl_interp_cspline, num);
-    else
-	spline = gsl_spline_alloc (gsl_interp_linear, num);
-
-    gsl_spline_init (spline, x, y, num);
-    min= x[0];
-    max= x[num-1];
-
-}
-InterpolatorGSL::InterpolatorGSL(const std::vector<double>& x, const std::vector<double>& y,interpolation_t type){
-    //Interpolator(x.data(),y.data(),x.size());
-    acc = gsl_interp_accel_alloc ();
-    if(type==cspline)
-	spline = gsl_spline_alloc (gsl_interp_cspline, x.size());
-    else
-	spline = gsl_spline_alloc (gsl_interp_linear, x.size());
-
-    gsl_spline_init (spline, x.data(), y.data(), x.size());
-    min= x[0];
-    max= x[x.size()-1];
-}
-
-InterpolatorGSL::~InterpolatorGSL(){
-    gsl_interp_accel_free (acc);
-    gsl_spline_free (spline);
-}
-
-double InterpolatorGSL::eval(double x){
-    if(x<min)x=min;
-    if(x>max)x=max;
-    return gsl_spline_eval(spline, x, acc);
-}
-
-double InterpolatorGSL::derivative(double x){
-    if(x<min)x=min;
-    if(x>max)x=max;
-    return gsl_spline_eval_deriv (spline, x, acc);
-}
-
-
-//////////// Interpolator2 ////////////////////////////////
-#ifdef BUILTIN_SPLINE
-Interpolator2::Interpolator2(const double *x, const double *y, int num){
-    ss.set_points(x,y,num);
-    min= x[0];
-    max= x[num-1];
-
-}
-
-double Interpolator2::eval(double x){
-    if(x<min)x=min;
-    if(x>max)x=max;
-    return ss(x);
-}
-
-double Interpolator2::derivative(double x){
-    if(x<min)x=min;
-    if(x>max)x=max;
-    return ss.deriv(1,x);
-}
-#endif
-
 }
diff --git a/storage.h b/storage.h
index 6baeba0..3abb641 100644
--- a/storage.h
+++ b/storage.h
@@ -18,24 +18,23 @@
 #define STORAGE
 
 #include <vector>
+#include <array>
 #include <iterator>
 #include <cmath>
 //#include <unordered_set>
-#include <gsl/gsl_spline.h>
 #include "catima/build_config.h"
 #include "catima/constants.h"
 #include "catima/structures.h"
 #include "catima/config.h"
 
-#ifdef BUILTIN_SPLINE
 #include "catima/spline.h"
-#endif
 
 
 namespace catima{
 
-	enum interpolation_t {cspline, linear};
-
+    /**
+     * Class to store energy points, log spaced from logmin to logmax.
+     */
 	template<int N>
     struct EnergyTable{
 		EnergyTable(double logmin, double logmax):values(),step(0.0),num(N){
@@ -45,10 +44,18 @@ namespace catima{
 		}
 	    }
 	double operator()(int i)const{return values[i];}
+    double operator[](int i)const{return values[i];}
+    static constexpr int size() {return N;};
 	double values[N];
 	double step;
 	double* begin(){return values;}
-	double* end(){return &values[num-1];}
+    double* end(){return &values[num];}
+    int index(double v)const noexcept{
+        double lxval = (log(v/values[0])/M_LN10);
+        if(v<values[0] || step==0.0)return -1;
+        if(v>=values[N-1])return N-1;
+        return static_cast<int> (std::floor(lxval/step));
+    };
 	std::size_t num;
     };
 
@@ -75,109 +82,146 @@ namespace catima{
 	    return r;
 	}
 
-    /*
-    template<int N>
-    struct EnergyTableLinear{
-	constexpr EnergyTableLinear():values(),num(N){
-	    for(auto i=0;i<N;i++){
-		values[i]=exp(M_LN10*(logEmin + ((double)i)*(logEmax-logEmin)/(N - 1.0)));
-		}
-	    }
-	double operator()(int i)const{return values[i];}
-	double values[N];
-	std::size_t num;
+    //////////////////////////////////////////////////////////////////////////////////////
+    #ifdef GSL_INTERPOLATION
+    /// Interpolation class, to store interpolated values
+        class InterpolatorGSL{
+            public:
+            InterpolatorGSL(){};
+            InterpolatorGSL(const EnergyTable<max_datapoints>& x, const std::vector<double>& y, interpolation_t type=cspline);
+            ~InterpolatorGSL();
+            double operator()(double x)const{return eval(x);};
+            double eval(double x) const;
+            double derivative(double x) const;
+            double get_min()const{return min;};
+            double get_max()const{return max;};
+
+            private:
+            double min=0;
+            double max=0;
+            gsl_interp_accel *acc;
+            gsl_spline *spline;
+        };
+    #endif
+
+    class InterpolatorCSpline{
+    public:
+        using xtype = EnergyTable<max_datapoints>;
+        InterpolatorCSpline()=default;
+        InterpolatorCSpline(const xtype &table, const std::vector<double> &y):
+            min(table.values[0]), max(table.values[max_datapoints-1]), ss(table,y){}
+        double operator()(double x)const{return eval(x);}
+        double eval(double x)const{return ss.evaluate(x);}
+        double derivative(double x)const{return ss.deriv(x);}
+        double get_min()const{return min;}
+        double get_max()const{return max;}
+
+    private:
+        double min=0;
+        double max=0;
+        cspline_special<xtype> ss;
     };
-    */
 
-    // return vector with lineary spaced elements from a to b, num is number of elements
-inline std::vector<double> linspace_vector(double a, double b, unsigned int num){
-    std::vector<double> res;
-    if(num>=2 && a<b){
-	res.resize(num);
-	double step = (b-a)/(num-1);
-	for(unsigned int i=0;i<(num-1);i++){
-	    res[i]=a+(i*step);
-	    }
-	res[num-1] = b;
-	}
-    return res;
-    }
+#ifdef GSL_INTERPOLATION
+using Interpolator = InterpolatorGSL;
+#else
+using Interpolator = InterpolatorCSpline;
+#endif
 
-    class DataPoint{
+#ifdef STORE_SPLINES
+    using spline_type = const Interpolator&;
+#else
+    using spline_type = Interpolator;
+#endif
+
+// return vector with lineary spaced elements from a to b, num is number of elements
+
+/**
+ * @brief structure to store calculated data points and optionally also splines
+ */
+class DataPoint{
 	public:
 	Projectile p;
 	Material m;
 	Config config;
-	std::vector<double> range;
-	std::vector<double> range_straggling;
-	std::vector<double> angular_variance;
 
-	DataPoint(){};
-	DataPoint(const Projectile _p, const Material _m,const Config &_c=default_config):p(_p),m(_m),config(_c){};
-	~DataPoint();
+    std::vector<double> range;
+    std::vector<double> range_straggling;
+    std::vector<double> angular_variance;
+#ifdef STORE_SPLINES
+    Interpolator range_spline;
+    Interpolator range_straggling_spline;
+    Interpolator angular_variance_spline;
+#endif
+    DataPoint()=default;
+    DataPoint(const Projectile _p, const Material _m,const Config &_c=default_config):p(_p),m(_m),config(_c){}
+    DataPoint(const DataPoint&)=delete;
+    DataPoint(DataPoint&&)=default;
+    DataPoint& operator=(const DataPoint&)=default;
+    DataPoint& operator=(DataPoint&&)=default;
 	friend bool operator==(const DataPoint &a, const DataPoint &b);
     };
 
-    class Data{
-	public:
-	Data();
-	~Data();
-	void Add(const Projectile &p, const Material &t, const Config &c=default_config);
-	int GetN() const {return storage.size();};
-	void Reset(){storage.clear();storage.resize(max_storage_data);index=storage.begin();};
-	DataPoint& Get(const Projectile &p, const Material &t, const Config &c=default_config);
-	DataPoint& Get(unsigned int i){return storage[i];};
-	int get_index() {return std::distance(storage.begin(),index);}
-	private:
-	std::vector<DataPoint> storage;
-	std::vector<DataPoint>::iterator index;
-    };
-
-/// Interpolation class, to store interpolated values
-	class InterpolatorGSL{
-        public:
-        InterpolatorGSL(const double *x, const double *y, int num,interpolation_t type=cspline);
-        InterpolatorGSL(const std::vector<double>& x, const std::vector<double>& y,interpolation_t type=cspline);
-        ~InterpolatorGSL();
-        double operator()(double x){return eval(x);};
-        double eval(double x);
-        double derivative(double x);
-        double get_min(){return min;};
-        double get_max(){return max;};
-
-        private:
-        double min=0;
-        double max=0;
-        gsl_interp_accel *acc;
-        gsl_spline *spline;
-    };
-
-#ifdef BUILTIN_SPLINE
-	class Interpolator2{
-        public:
-        Interpolator2(const double *x, const double *y, int num);
-        double operator()(double x){return eval(x);};
-        double eval(double x);
-        double derivative(double x);
-        double get_min(){return min;};
-        double get_max(){return max;};
-
-        private:
-        double min=0;
-        double max=0;
-        spline ss;
-    };
+#ifdef STORE_SPLINES
+    const Interpolator& get_range_spline(const DataPoint &data);
+    const Interpolator& get_range_straggling_spline(const DataPoint &data);
+    const Interpolator& get_angular_variance_spline(const DataPoint &data);
+#else
+    Interpolator get_range_spline(const DataPoint &data);
+    Interpolator get_range_straggling_spline(const DataPoint &data);
+    Interpolator get_angular_variance_spline(const DataPoint &data);
 #endif
+
+/**
+ * @brief The Data class to store DataPoints
+ */
+    class Data{
+    public:
+        Data();
+        ~Data();
+
+        /**
+         * @brief Add new DataPoint
+         * @param p - Projectile
+         * @param t - Material
+         * @param c - Config
+         */
+        void Add(const Projectile &p, const Material &t, const Config &c=default_config);
+
+        int GetN() const {return storage.size();};
+        void Reset(){storage.clear();storage.resize(max_storage_data);index=storage.begin();};
+
+        /**
+         * @brief Get DataPoint reference for projectile-target-config combination
+         * @param p - Projectile
+         * @param t - Material
+         * @param c - Config
+         * @return reference to DataPoint
+         */
+        DataPoint& Get(const Projectile &p, const Material &t, const Config &c=default_config);
+        DataPoint& Get(unsigned int i){return storage[i];};
+        int get_index() {return std::distance(storage.begin(),index);}
+
+    private:
+        std::vector<DataPoint> storage;
+        std::vector<DataPoint>::iterator index;
+    };
+
     extern Data _storage;
 
-	inline DataPoint& get_data(const Projectile &p, const Material &t, const Config &c=default_config){
+    /**
+     * @brief get_data - Get DataPoint from the global storage class
+     * @param p - Projectile
+     * @param t - Material
+     * @param c - Config
+     * @return const reference to DataPoint
+     */
+    inline const DataPoint& get_data(const Projectile &p, const Material &t, const Config &c=default_config){
 		return _storage.Get(p, t, c);
 	}
 
     bool operator==(const DataPoint &a, const DataPoint &b);
 
-    using InterpolatorLinear = InterpolatorGSL;
-    using Interpolator = InterpolatorGSL;
 }
 
 #endif
diff --git a/tests/test.py b/tests/test.py
index da604c6..89cd010 100644
--- a/tests/test.py
+++ b/tests/test.py
@@ -217,6 +217,7 @@ class TestStructures(unittest.TestCase):
         energies = [100,500,1000]
         res2 = catima.dedx_from_range(p,graphite,energy=energies)
         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)
             self.assertAlmostEqual(res2[i], r, 0.1)
diff --git a/tests/test2.py b/tests/test2.py
new file mode 100644
index 0000000..af6c087
--- /dev/null
+++ b/tests/test2.py
@@ -0,0 +1,309 @@
+import sys
+sys.path.insert(0,"../build")
+import unittest
+import pycatima as catima
+import math
+
+class TestStructures(unittest.TestCase):
+
+    def test_Projectile(self):
+        print(catima.storage_info())
+        p = catima.Projectile(238,92)
+        self.assertEqual(p.A(),238)
+        self.assertEqual(p.Z(),92)
+        self.assertEqual(p.Q(),92)
+
+        p = catima.Projectile(238,92,90)
+        self.assertEqual(p.A(),238)
+        self.assertEqual(p.Z(),92)
+        self.assertEqual(p.Q(),90)
+        p.T(1000)
+        self.assertEqual(p.T(),1000)
+        p(500)
+        self.assertEqual(p.T(),500)
+
+        p = catima.Projectile(238,92,90, T=100)
+        self.assertEqual(p.T(),100)
+    
+    def test_Material(self):
+        mat = catima.Material()
+        mat.add_element(12,6,1)
+        self.assertEqual(mat.ncomponents(),1)
+        mat.add_element(1,1,2)
+        self.assertEqual(mat.ncomponents(),2)
+
+        mat2 = catima.Material([[12.01,6,1]])
+        self.assertEqual(mat2.ncomponents(),1)
+        self.assertEqual(mat2.molar_mass(),12.01)
+
+        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]])
+        self.assertEqual(mat2.ncomponents(),1)
+        self.assertAlmostEqual(mat2.molar_mass(),12,1)
+
+        mat5 = catima.Material([[0,6,1]],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)
+        self.assertEqual(mat6.ncomponents(),1)
+        self.assertEqual(mat6.thickness(),0.5)
+        self.assertEqual(mat6.density(),1.9)
+        self.assertEqual(mat6.I(),80.0)
+
+        # copy
+        mat3.density(1.8)
+        matc = mat3.copy()
+        self.assertEqual(matc.ncomponents(),1)
+        self.assertEqual(matc.molar_mass(),12)
+        self.assertEqual(matc.density(),1.8)
+        mat3.density(2.0)
+        self.assertEqual(matc.density(),1.8)
+        self.assertEqual(mat3.density(),2.0)
+
+    
+    def test_default_material(self):
+        m1 = catima.get_material(6);
+        self.assertAlmostEqual(m1.molar_mass(),12,1)
+        self.assertEqual(m1.ncomponents(),1)
+        self.assertAlmostEqual(m1.density(),2.0,1)
+
+        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)
+
+        m3 = catima.get_material(3001)
+        self.assertEqual(m3.ncomponents(),0)
+        self.assertAlmostEqual(m3.molar_mass(),0,1)
+        self.assertAlmostEqual(m3.density(),0.0,1)
+
+    def test_layers(self):
+        graphite = catima.get_material(6)
+        graphite.thickness(0.5)
+        p10 = catima.get_material(catima.material.P10)
+        p10.thickness(0.01)
+        n2 = catima.get_material(7)
+        n2.thickness(0.02)
+
+        mat= catima.Layers()
+        self.assertEqual(mat.num(),0)
+        mat.add(graphite)
+        self.assertEqual(mat.num(),1)
+        self.assertAlmostEqual(mat[0].molar_mass(),12,1)
+        self.assertAlmostEqual(mat[0].thickness(),0.5,1)
+        self.assertAlmostEqual(mat[0].density(),2.0,1)
+        
+        mat.add(p10)
+        self.assertEqual(mat.num(),2)
+
+        graphite.thickness(1.0)
+        graphite.density(1.8)
+        mat.add(graphite)
+        self.assertEqual(mat.num(),3)
+        self.assertAlmostEqual(mat[2].molar_mass(),12,1)
+        self.assertAlmostEqual(mat[0].thickness(),0.5,1)
+        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)
+
+        mat2 = catima.Layers()
+        mat2.add(n2)
+        self.assertEqual(mat2.num(),1)
+
+        mats = mat2 + mat
+        self.assertEqual(mats.num(),4)
+        self.assertAlmostEqual(mats[0].molar_mass(),14,1)
+        self.assertEqual(mats[0].thickness(),0.02)
+        self.assertAlmostEqual(mats[1].molar_mass(),12,1)
+        self.assertAlmostEqual(mats[3].molar_mass(),12,1)
+
+        n2.thickness(0.5)
+        mats = mats + n2
+        self.assertEqual(mats.num(),5)
+        self.assertAlmostEqual(mats[0].molar_mass(),14,1)
+        self.assertEqual(mats[0].thickness(),0.02)
+        self.assertAlmostEqual(mats[4].molar_mass(),14,1)
+        self.assertEqual(mats[4].thickness(),0.5)
+
+    def test_material_calculation(self):
+        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)
+        self.assertAlmostEqual(res.dEdxi,2.76,1)
+        self.assertAlmostEqual(res.dEdxi,res2,3)
+
+        res = catima.calculate(p(9),Water)
+        res2 = catima.dedx(p,Water)
+        self.assertAlmostEqual(res.dEdxi,51.17,1)
+        self.assertAlmostEqual(res.dEdxi,res2,3)
+        res = catima.calculate(p(9),Water)
+        res = catima.calculate(p(9),Water)
+        self.assertAlmostEqual(res.dEdxi,51.17,1)
+        
+        p(900000)
+        res = catima.calculate(p,Water)
+        res2 = catima.dedx_from_range(p,Water)
+        self.assertAlmostEqual(res.dEdxi,res2,3)
+        
+    def test_config(self):
+        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
+        conf2 = catima.Config()
+        conf2.dedx_straggling = catima.omega_type.atima
+        p(1000)
+        res = catima.calculate(p,Water,config=conf)
+        res2 = catima.calculate(p,Water,config=conf2)
+        self.assertAlmostEqual(res.dEdxi,res2.dEdxi,delta=1e-6)
+        self.assertNotAlmostEqual(res.sigma_E,res2.sigma_E,delta=1e-4)
+        self.assertNotAlmostEqual(res.sigma_r,res2.sigma_r,delta=1e-4)
+        
+    
+    def test_eout(self):
+        graphite = catima.get_material(6)
+        graphite.thickness(0.5)
+        p = catima.Projectile(12,6)
+        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):
+        graphite = catima.get_material(6)
+        graphite.thickness(0.5)
+        p = catima.Projectile(12,6)
+        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)
+        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)
+        self.assertAlmostEqual(res2[1], catima.calculate(p(energies[1]),graphite).Eout ,1)
+    
+    def test_dedx_from_range_list(self):
+        graphite = catima.get_material(6)
+        graphite.thickness(0.5)
+        p = catima.Projectile(12,6)
+        energies = [100,500,1000]
+        res2 = catima.dedx_from_range(p,graphite,energy=energies)
+        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)
+            print(r)
+            print(res2)
+            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.thickness(10.0)
+        graphite = catima.get_material(6)
+        graphite.thickness(1.0)
+        graphite.density(2.0)
+
+        mat = catima.Layers()
+        mat.add(Water)
+        mat.add(graphite)
+        res = catima.calculate_layers(p(1000),mat)
+        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[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)
+    
+    def test_energy_table(self):
+        table = catima.get_energy_table()
+        self.assertEqual(table[0],math.log(catima.logEmin))
+        #self.assertEqual(table[10],catima.energy_table(10))
+        self.assertEqual(len(table),catima.max_datapoints)
+        
+    def test_storage(self):
+        p = catima.Projectile(12,6)
+        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)
+        print(data[0])
+        et = catima.get_energy_table()
+        
+        self.assertEqual(len(data),3)
+        self.assertEqual(len(data[0]),len(et))
+        self.assertEqual(len(data[0]),catima.max_datapoints)
+        
+        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.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.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.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.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.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
+        r = catima.storage_info()
+        
+        #self.assertAlmostEqual(catima.da2de(p,Water,et[100]),data[2][100],6)
+        #self.assertAlmostEqual(catima.da2de(p,Water,et[400]),data[2][400],6)
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/tests/test_storage.cpp b/tests/test_storage.cpp
index cf95738..2908f7c 100644
--- a/tests/test_storage.cpp
+++ b/tests/test_storage.cpp
@@ -57,18 +57,18 @@ const lest::test specification[] =
       EXPECT(catima::_storage.get_index()==0);
       
       catima::_storage.Add(p,water);
-      auto dp = catima::_storage.Get(0);
+      auto& dp = catima::_storage.Get(0);
       EXPECT(catima::_storage.get_index()==1);
       EXPECT(dp.p.A==12);
       EXPECT(dp.m.ncomponents()==2);
       catima::_storage.Add(p,water);
-      auto dp2 = catima::_storage.Get(1);
+      auto& dp2 = catima::_storage.Get(1);
       EXPECT(catima::_storage.get_index()==1);
       EXPECT(dp2.p.A==0);
       EXPECT(dp2.m.ncomponents()==0);
       
       catima::_storage.Add(p,graphite);
-      auto dp3 = catima::_storage.Get(1);
+      auto&  dp3 = catima::_storage.Get(1);
       EXPECT(catima::_storage.get_index()==2);
       EXPECT(dp3.p.A==12);
       EXPECT(dp3.m.ncomponents()==1);