diff --git a/python/thor_scsi/utils/accelerator.py b/python/thor_scsi/utils/accelerator.py
index a3f75bfa..2d7784e6 100644
--- a/python/thor_scsi/utils/accelerator.py
+++ b/python/thor_scsi/utils/accelerator.py
@@ -101,10 +101,11 @@ def extract_orbit_from_standard_observers(
     tps_tmp = [ob.get_truncated_power_series_a() for ob in ob_sub]
     logger.debug(tps_tmp)
     tps = xr.DataArray(
-        data=[t.jacobian() for t in tps_tmp],
+        data=np.array([t.jacobian() for t in tps_tmp], float),
         name="tps_jacobian",
         dims=["index", "phase_coordinate_row", "phase_coordinate_col"],
-        coords=[indices, phase_space_coords_names, phase_space_coords_names],
+# todo: don't understand where this comes from
+        coords=[indices, phase_space_coords_names, phase_space_coords_names + ["unknown"]],
     )
     ps = xr.DataArray(
         data=[np.array(t.cst().iloc) for t in tps_tmp],
diff --git a/src/thor_scsi/CMakeLists.txt b/src/thor_scsi/CMakeLists.txt
index 89e5fff1..c3922467 100644
--- a/src/thor_scsi/CMakeLists.txt
+++ b/src/thor_scsi/CMakeLists.txt
@@ -179,55 +179,12 @@ target_link_libraries(test_transform_phase_space
 
 add_test(transform_phase_space test_transform_phase_space)
 
-## add_executable(a_kick  elements/a_kick.cc)
-## target_include_directories(a_kick
-##     PUBLIC
-##     "$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/../>"
-##     "$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>"
-## )
-## target_link_libraries(a_kick
-##   thor_scsi
-##   thor_scsi_core
-##   tpsa_lin
-##   # PRIVATE Boost::boost
-##   #  quadmath
-##     ${Boost_PRG_EXEC_MONITOR_LIBRARY}
-##     ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}
-## )
-
-#
-# Non linear part needs still to be ported to new interface
-## set(tpsa_nlin_FILES
-##     tpsa_for_pm.cc
-##     # tpsa_for.cc
-## )
-##
-## add_library(tpsa_nlin SHARED
-##   ${tpsa_common_FILES}
-##   ${tpsa_nlin_FILES}
-##   ${tpsa_HEADERS}
-## )
-
-
-
-## set(thor_scsi_CORE_HEADERS
-##   thor/core/util.h
-##   flame/core/base.h
-##   flame/core/config.h
-##   )
-##
-##
-##
-## add_library(thor_scsi SHARED
-##   ${thor_scsi_CORE_files}
-##   ${thor_scsi_CORE_HEADERS}
-## )
-
 
 add_subdirectory(elements)
 add_subdirectory(std_machine)
 add_subdirectory(examples)
 add_subdirectory(custom)
+add_subdirectory(core)
 
 # ---- install helpers --
 
diff --git a/src/thor_scsi/core/CMakeLists.txt b/src/thor_scsi/core/CMakeLists.txt
index 72f92f09..dc2aa76a 100644
--- a/src/thor_scsi/core/CMakeLists.txt
+++ b/src/thor_scsi/core/CMakeLists.txt
@@ -1,23 +1,39 @@
 # could be added to file above
 # but not yet required
 # lets leave it here
+message(STATUS "thor_scsi::core gtpsa cpp include dir ${gtpsa_cpp_INCLUDE_DIR}")
+
+add_executable(test_math_comb
+        test_math_comb.cc
+)
+target_include_directories(test_math_comb
+        PUBLIC
+        "$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/../../>"
+)
+target_link_libraries(test_math_comb
+        thor_scsi_core
+        ${Boost_PRG_EXEC_MONITOR_LIBRARY}
+        ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}
+)
+
+add_test(math_comb test_math_comb)
 
 add_executable(test_transform
-    test_transform.cc
+        test_transform.cc
 )
 
-message(STATUS "thor_scsi::core gtpsa cpp include dir ${gtpsa_cpp_INCLUDE_DIR}")
+
 target_include_directories(test_transform
-    PUBLIC
-    "$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/../../>"
-    "$<BUILD_INTERFACE:${gtpsa_cpp_INCLUDE_DIR}>"
-    "$<BUILD_INTERFACE:${gtpsa_mad_ng_INCLUDE_DIR}>"
-    "$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>"
+        PUBLIC
+        "$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/../../>"
+        "$<BUILD_INTERFACE:${gtpsa_cpp_INCLUDE_DIR}>"
+        "$<BUILD_INTERFACE:${gtpsa_mad_ng_INCLUDE_DIR}>"
+        "$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>"
 )
 
 target_link_libraries(test_transform
-    ${Boost_PRG_EXEC_MONITOR_LIBRARY}
-    ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}
+        ${Boost_PRG_EXEC_MONITOR_LIBRARY}
+        ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}
 )
 
 add_test(transform test_transform)
@@ -38,12 +54,8 @@ target_include_directories(test_two_dimensional_multipoles
     "$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>"
 )
 target_link_libraries(test_two_dimensional_multipoles
-  # PRIVATE Boost::boost
-  #  quadmath
-  # gsl
-  # gslcblas
-  gtpsa::c++
-  gtpsa
-    ${Boost_PRG_EXEC_MONITOR_LIBRARY}
-    ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}
+        gtpsa::gtpsa-c++
+        gtpsa::gtpsa
+        ${Boost_PRG_EXEC_MONITOR_LIBRARY}
+        ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}
 )
diff --git a/src/thor_scsi/core/math_comb.cc b/src/thor_scsi/core/math_comb.cc
index 82db350f..cb06581f 100644
--- a/src/thor_scsi/core/math_comb.cc
+++ b/src/thor_scsi/core/math_comb.cc
@@ -1,56 +1,60 @@
 #include <thor_scsi/core/math_comb.h>
 #include <limits.h>
-#include <cassert>
+#include <stdexcept>
+#include <sstream>
 
 namespace tsc = thor_scsi;
 
 unsigned long thor_scsi::core::binom(unsigned long n, unsigned long k) {
-	unsigned long c = 1, i;
-
-	if (k > n-k) // take advantage of symmetry
-		k = n-k;
-
-	for (i = 1; i <= k; i++, n--) {
-		if (c/i >= ULONG_MAX/n){
-			// raise an exception
-			assert(0);
-			return 0;
-		}
-		// split c * n / i into (c / i * i + c % i) * n / i
-		c = c / i * n + c % i * n / i;
-	}
-
-	return c;
+        unsigned long c = 1, i;
+
+        if (k > n-k) // take advantage of symmetry
+                k = n-k;
+
+        for (i = 1; i <= k; i++, n--) {
+                if (c/i >= ULONG_MAX/n){
+                        std::stringstream strm;
+                        strm << "(c/i = "<< c/i << ") >= (ULONG_MAX/n " << ULONG_MAX/n << ")."
+			     << " n = " << n << ", k = " << k << ", c = " << c
+                                ;
+                        throw std::overflow_error(strm.str());
+                        return 0;
+                }
+                // split c * n / i into (c / i * i + c % i) * n / i
+                c = c / i * n + c % i * n / i;
+        }
+
+        return c;
 }
 
 #if 0
 double thor_scsi::binom (int const n, int const k)
 {
-	int i;
-	double num, den;
+        int i;
+        double num, den;
 
-	if (k<0  || n<0 || k>n){
-		return -1.0;
-	}
+        if (k<0  || n<0 || k>n){
+                return -1.0;
+        }
 
-	if (k==0 || k==n){
-		return 1.0;
-	}
+        if (k==0 || k==n){
+                return 1.0;
+        }
 
-	if (k==1 || k==n-1){
-		return n;
-	}
+        if (k==1 || k==n-1){
+                return n;
+        }
 
-	num = den = 1.0;
+        num = den = 1.0;
 
-	for(i=n-k+1; i<=n; i++){
-		num *= i;
-	}
-	for(i=1;     i<=k; i++){
-		den *= i;
-	}
+        for(i=n-k+1; i<=n; i++){
+                num *= i;
+        }
+        for(i=1;     i<=k; i++){
+                den *= i;
+        }
 
-	return num/den;
+        return num/den;
 }
 #endif
 /*
diff --git a/src/thor_scsi/core/math_comb.h b/src/thor_scsi/core/math_comb.h
index a80d75ec..60daa762 100644
--- a/src/thor_scsi/core/math_comb.h
+++ b/src/thor_scsi/core/math_comb.h
@@ -1,12 +1,12 @@
 #ifndef _THOR_SCSI_MATH_COMB_H_
 #define _THOR_SCSI_MATH_COMB_H_
 namespace thor_scsi::core {
-	/**
-	 * @brief calculate binomial coefficient
-	 *
-	 * @ todo: check if in standard library or in GSL
-	 */
-	unsigned long binom (unsigned long const n, unsigned long const k);
+        /**
+         * @brief calculate binomial coefficient
+         *
+         * @todo: check if in standard library or in GSL
+         */
+        unsigned long binom (unsigned long const n, unsigned long const k);
 }
 #endif /* THOR_SCSI_MATH_COMB_H */
 /*
diff --git a/src/thor_scsi/core/test_math_comb.cc b/src/thor_scsi/core/test_math_comb.cc
new file mode 100644
index 00000000..a8b6ed78
--- /dev/null
+++ b/src/thor_scsi/core/test_math_comb.cc
@@ -0,0 +1,57 @@
+#define BOOST_TEST_MODULE transform
+#define BOOST_TEST_DYN_LINK
+
+#include <boost/test/unit_test.hpp>
+#include <thor_scsi/core/math_comb.h>
+
+auto binom = thor_scsi::core::binom;
+
+BOOST_AUTO_TEST_CASE(test10_binom_pascal_triangle)
+{
+    //BOOST_TEST(binom(0, 1) == 0);
+    //BOOST_TEST(binom(0, 2) == 0);
+    //BOOST_TEST(binom(0, 3) == 0);
+
+    // the top
+    BOOST_TEST(binom(1, 1) == 1);
+
+    // quadratic expansion
+    BOOST_TEST(binom(2, 0) == 1);
+    BOOST_TEST(binom(2, 1) == 2);
+    BOOST_TEST(binom(2, 2) == 1);
+
+    // cubic expansion
+    BOOST_TEST(binom(3, 0) == 1);
+    BOOST_TEST(binom(3, 1) == 3);
+    BOOST_TEST(binom(3, 2) == 3);
+    BOOST_TEST(binom(3, 3) == 1);
+
+    // 4th order  expansion
+    BOOST_TEST(binom(4, 0) == 1);
+    BOOST_TEST(binom(4, 1) == 4);
+    BOOST_TEST(binom(4, 2) == 6);
+    BOOST_TEST(binom(4, 3) == 4);
+    BOOST_TEST(binom(4, 4) == 1);
+
+    // 5th order  expansion
+    BOOST_TEST(binom(5, 0) == 1);
+    BOOST_TEST(binom(5, 1) == 5);
+    BOOST_TEST(binom(5, 2) == 10);
+    BOOST_TEST(binom(5, 3) == 10);
+    BOOST_TEST(binom(5, 4) == 5);
+    BOOST_TEST(binom(5, 5) == 1);
+
+    // 6th order  expansion
+    BOOST_TEST(binom(6, 0) == 1);
+    BOOST_TEST(binom(6, 1) == 6);
+    BOOST_TEST(binom(6, 2) == 15);
+    BOOST_TEST(binom(6, 3) == 20);
+    BOOST_TEST(binom(6, 4) == 15);
+    BOOST_TEST(binom(6, 5) == 6);
+    BOOST_TEST(binom(6, 6) == 1);
+}
+
+BOOST_AUTO_TEST_CASE(test10_binom_too_large)
+{
+    BOOST_CHECK_THROW(binom(2000 * 1000, 3000), std::overflow_error);
+}
\ No newline at end of file
diff --git a/src/thor_scsi/core/test_two_dimensional_multipoles.cc b/src/thor_scsi/core/test_two_dimensional_multipoles.cc
index 07fed58a..66c57ef0 100644
--- a/src/thor_scsi/core/test_two_dimensional_multipoles.cc
+++ b/src/thor_scsi/core/test_two_dimensional_multipoles.cc
@@ -3,7 +3,6 @@
 #include <boost/test/unit_test.hpp>
 #include <boost/test/tools/output_test_stream.hpp>
 #include <thor_scsi/core/multipoles.h>
-#include <gtpsa/utils_tps.hpp>
 #include <cmath>
 #include <array>
 
@@ -266,18 +265,6 @@ BOOST_AUTO_TEST_CASE(test20_set_harmonic_quadrupole)
 			  << "By\n " << By << "\n";
 
 	}
-	{
-
-		auto x = tps(0, 0+1);
-		auto y = tps(0, 2+1);
-		tps Bx, By;
-
-		std::cout << "x\n " << x << "\n"
-			  << "y\n " << y << "\n";
-		h.field(x, y, &Bx, &By);
-		std::cout << "Bx\n " << Bx << "\n"
-			  << "By\n " << By << "\n";
-	}
 	{
 		auto x = t_ref.clone(), y = t_ref.clone(), Bx = t_ref.clone(), By = t_ref.clone();
 		x = 1;