From e06cdd31711c8ec9a48c91872b4d1f3cc8a2f35c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= <marek@necada.org>
Date: Wed, 14 Mar 2018 17:17:52 +0200
Subject: [PATCH] =?UTF-8?q?Dudom.=20Implementace=20besselov=C3=BDch=20fc?=
 =?UTF-8?q?=C3=AD=20v=20chodu.?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Former-commit-id: 63827799543f410dfb1a83cab7ca10287d76f7c2
---
 qpms/bessel.c        | 16 ++++++++++++++++
 qpms/qpms_specfunc.h | 14 ++++++++++++++
 2 files changed, 30 insertions(+)

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"                 *
  ******************************************************************************/