From 8001b67603d9df49e110a4d4e2e2038effe00edb Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= <marek@necada.org>
Date: Thu, 18 Jan 2018 11:22:27 +0000
Subject: [PATCH] [dirty dipoles] bessel function coefficient table

Former-commit-id: 99679b765621c7237bce927169cf83331fe61c99
---
 dipdip-dirty/bessels.c | 40 ++++++++++++++++++++++++++++++++++++++++
 dipdip-dirty/bessels.h | 18 ++++++++++++++++++
 2 files changed, 58 insertions(+)
 create mode 100644 dipdip-dirty/bessels.c
 create mode 100644 dipdip-dirty/bessels.h

diff --git a/dipdip-dirty/bessels.c b/dipdip-dirty/bessels.c
new file mode 100644
index 0000000..cd49680
--- /dev/null
+++ b/dipdip-dirty/bessels.c
@@ -0,0 +1,40 @@
+#include "bessels.h"
+#include <stdlib.h>
+#include <math.h>
+
+static const double ln2 = 0.69314718055994531;
+
+#include <stdio.h>
+
+// general; gives an array of size 
+complex double * hankelcoefftable_init(size_t maxn) {
+	complex double *hct = malloc((maxn+1)*(maxn+2)/2 * sizeof(complex double));
+	for(size_t n = 0; n <= maxn; ++n) {
+		complex double *hcs = hankelcoeffs_get(hct,n);
+		for (size_t k = 0; k <= n; ++k) {
+			double lcoeff = lgamma(n+k+1) - lgamma(n-k+1) - lgamma(k+1) - k*ln2;
+			printf("%f, %.16f\n", lcoeff, exp(lcoeff));
+			// for some reason, casting k-n to double does not work,so
+			// cpow (I, k-n-1) cannot be used...
+			complex double ifactor;
+			switch ((n+1-k) % 4) {
+				case 0:
+					ifactor = 1;
+					break;
+				case 1:
+					ifactor = -I;
+					break;
+				case 2:
+					ifactor = -1;
+					break;
+				case 3:
+					ifactor = I;
+					break;
+			}
+			// the result should be integer, so round to remove inaccuracies
+			hcs[k] = round(exp(lcoeff)) * ifactor;
+		}
+	}
+	return hct;
+}
+
diff --git a/dipdip-dirty/bessels.h b/dipdip-dirty/bessels.h
new file mode 100644
index 0000000..0553b17
--- /dev/null
+++ b/dipdip-dirty/bessels.h
@@ -0,0 +1,18 @@
+#ifndef BESSELS_H
+#define BESSELS_H
+
+#include <stddef.h>
+#include <complex.h>
+
+complex double *hankelcoefftable_init(size_t maxn);
+
+// general, gives the offset such that result[k] is 
+// the coefficient corresponding to the e**(I * x) * x**(-k-1)
+// term of the Hankel function; no boundary checks!
+static inline complex double *
+hankelcoeffs_get(complex double *hankelcoefftable, size_t n){
+        return hankelcoefftable + n*(n+1)/2;
+}
+
+
+#endif //BESSELS_H