diff --git a/qpms/bessel.c b/qpms/bessel.c
index 162a0f2..468f2f0 100644
--- a/qpms/bessel.c
+++ b/qpms/bessel.c
@@ -1,8 +1,11 @@
 #include <assert.h>
 #include "qpms_specfunc.h"
 #include <stdlib.h>
+#include <stddef.h>
+#include <kahansum.h>
 #include <gsl/gsl_sf_bessel.h>
 
+// There is a big issue with gsl's precision of spherical bessel function; these have to be implemented differently
 qpms_errno_t qpms_sph_bessel_fill(qpms_bessel_t typ, qpms_l_t lmax, double x, complex double *result_array) {
 	int retval;
 	double tmparr[lmax+1];
@@ -36,3 +39,16 @@ qpms_errno_t qpms_sph_bessel_fill(qpms_bessel_t typ, qpms_l_t lmax, double x, co
 	assert(0);
 }
 
+static inline ptrdiff_t akn_index(qpms_l_t n, qpms_l_t k) {
+	assert(k <= n);
+	return ((ptrdiff_t) n + 1) * n / 2 + k;
+}
+static inline ptrdiff_t bkn_index(qpms_l_t n, qpms_l_t k) {
+	assert(k <= n+1);
+	return ((ptrdiff_t) n + 2) * (n + 1) / 2 - 1 + k;
+}
+
+
+qpms_bessel_calculator_t *qpms_bessel_calculator_init(...){
+	TODO dudom;
+}
diff --git a/qpms/qpms_specfunc.h b/qpms/qpms_specfunc.h
index 1f50262..75a5ca7 100644
--- a/qpms/qpms_specfunc.h
+++ b/qpms/qpms_specfunc.h
@@ -11,6 +11,20 @@
 qpms_errno_t qpms_sph_bessel_fill(qpms_bessel_t typ, qpms_l_t lmax, double x, complex double *result_array);
 
 
+
+typedef struct qpms_bessel_calculator_t {
+        qpms_l_t lMax;
+        double *akn; // coefficients as in DLMF 10.49.1
+	complex double *bkn; // coefficients of the derivatives
+};
+
+qpms_bessel_calculator_t *qpms_bessel_calculator_init(qpms_l_t lMax);
+void qpms_bessel_calculator_free(qpms_bessel_calculator_t *c);
+
+qpms_errno_t qpms_sph_bessel_calc_fill(qpms_bessel_calculator_t *c, qpms_bessel_t typ, qpms_l_t lmax,
+	double x, complex double *result_array);
+
+
 /******************************************************************************
  *         Legendre functions and their "angular derivatives"                 *
  ******************************************************************************/