[dirty dipoles] Hankel transform of the long-range part of Hankel

functions, dimensionful

TODO more comments, dimensionless everything


Former-commit-id: f05f55c3d444bb633b0fe4877120ea001dd34aa7

dipdip-dirty/bessels.h

@@ -6,27 +6,35 @@
 
 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){
+trindex_cd(complex double *arr, size_t n){
-  return hankelcoefftable + n*(n+1)/2;
+  return arr + n*(n+1)/2;
 }
 
+// general, gives the offset such that result[ql] is 
+// the coefficient corresponding to the e**(I * x) * x**(-ql-1)
+// term of the n-th Hankel function; no boundary checks!
+static inline complex double *
+hankelcoeffs_get(complex double *hankelcoefftable, size_t n){
+  return trindex_cd(hankelcoefftable, n);
+}
 
 // general; target_longrange and target_shortrange are of size (maxn+1)
 // if target_longrange is NULL, only the short-range part is calculated
 void hankelparts_fill(complex double *target_longrange, complex double *target_shortrange,
 		size_t maxn, size_t longrange_order_cutoff, // x**(-(order+1)-1) terms go completely to short-range part
 		complex double *hankelcoefftable,
-		unsigned kappa, double c, double x); // x = k0 * r
+		unsigned kappa, double vc, double x); // x = k0 * r
 
 
 // this declaration is general; however, the implementation 
 // is so far only for kappa == ???, maxn == ??? TODO
-void lrhankel_recpart_fill(complex  double *target_longrange_kspace,
+void lrhankel_recpart_fill(complex  double *target_longrange_kspace /*Must be of size maxn*(maxn+1)/2*/,
-	       size_t maxn, size_t longrange_k_cutoff, 
+	       size_t maxp, size_t longrange_k_cutoff /* terms e**(I x)/x**(k+1), k>= longrange_k_cutoff go
+							 completely to the shortrange part
+							 index with hankelcoeffs_get(target,p)l[delta_m] */, 
 	       complex double  *hankelcoefftable,
 	       unsigned kappa, double c, double k0, double k);
 

dipdip-dirty/lrhankel_recspace_dirty.c

@@ -1,18 +1,38 @@
-#include "bessels.h"
+/*
-//#include "mdefs.h"
+ * This is a dirty implementation of lrhankel_recpart_fill() that calculates
-#include <complex.h>
+ * the (cylindrical) Hankel transforms of the regularised part of the spherical Hankel
-#include <string.h>
+ * functions that are to be summed in the reciprocal space.
-#define SQ(x) ((x)*(x))
+ *
+ * For now, only the regularisation with κ == 5 && q <= 2 && n <= 5 is implemented
+ * by writing down the explicit formula for each q,n pair and k>k0 vs k<k0 case,
+ * only with the help of some macros to make the whole thing shorter.
+ * 
+ * N.B. the results for very small k/k0 differ significantly (sometimes even in the first
+ * digit for n >= 3, probably due to catastrophic cancellation (hopefully not due
+ * to an error in the formula!). On the other hand,
+ * these numbers are tiny in their absolute value, so their contribution to the
+ * lattice sum should be negligible.
+ */
 
 #define MAXQM 1
 #define MAXN 5
 #define MAXKAPPA 5
 
-#define FF (-1)
+#include "bessels.h"
+//#include "mdefs.h"
+#include <complex.h>
+#include <string.h>
+#include <assert.h>
+#define SQ(x) ((x)*(x))
 
 #define P4(x) (((x)*(x))*((x)*(x)))
 
 
+/* 
+ * General form of the κ == 5 transforms. One usually has to put a (-1) factor
+ * to some part of the zeroth term of one of the cases k < k0 or k > k0 in order
+ * to stay on the correct branch of complex square root...
+ */
 #define KAPPA5SUM(form) (\
                 (form(0, 1))  \
             - 5*(form(1, 1)) \
@@ -31,11 +51,16 @@
             -   (form(5, 1)) \
 )
 
-#define LRHANKELDEF(fname) complex double fname(double c, double k0, double k, \
+/* 
+ * Prototype for the individual (per q,n) Bessel transform calculating functions.
+ * a, b, d, e, ash are recurring pre-calculated intermediate results, see the definition
+ * of lrhankel_recpart_fill() below to see their meaning
+ */
+#define LRHANKELDEF(fname) complex double fname(const double c, const double k0, const double k, \
 		const complex double *a, const complex double *b, const complex double *d, \
 		const complex double *e, const complex double *ash)
 
-typedef complex double (*lrhankelspec)(double, double, double, 
+typedef complex double (*lrhankelspec)(const double, const double, const double, 
 		const complex double *,
 		const complex double *,
 		const complex double *,
@@ -185,7 +210,6 @@ LRHANKELDEF(fk5q2n5s){
 #undef FORMK5Q2N5
 
 
-
 static const lrhankelspec transfuns_f[MAXKAPPA+1][MAXQM+1][MAXN+1] = {
 	{{NULL,NULL,NULL,NULL,NULL,NULL},{NULL,NULL,NULL,NULL,NULL,NULL}},
 	{{NULL,NULL,NULL,NULL,NULL,NULL},{NULL,NULL,NULL,NULL,NULL,NULL}},
@@ -195,7 +219,7 @@ static const lrhankelspec transfuns_f[MAXKAPPA+1][MAXQM+1][MAXN+1] = {
 	{{fk5q1n0l,fk5q1n1l,fk5q1n2l,fk5q1n3l,fk5q1n4l,fk5q1n5l},{fk5q2n0,fk5q2n1l,fk5q2n2l,fk5q2n3l,fk5q2n4l,fk5q2n5l}}
 };
 
-static lrhankelspec transfuns_n[MAXKAPPA+1][MAXQM+1][MAXN+1] = {
+static const lrhankelspec transfuns_n[MAXKAPPA+1][MAXQM+1][MAXN+1] = {
 	{{NULL,NULL,NULL,NULL,NULL,NULL},{NULL,NULL,NULL,NULL,NULL,NULL}},
 	{{NULL,NULL,NULL,NULL,NULL,NULL},{NULL,NULL,NULL,NULL,NULL,NULL}},
 	{{NULL,NULL,NULL,NULL,NULL,NULL},{NULL,NULL,NULL,NULL,NULL,NULL}},
@@ -203,19 +227,33 @@ static lrhankelspec transfuns_n[MAXKAPPA+1][MAXQM+1][MAXN+1] = {
 	{{NULL,NULL,NULL,NULL,NULL,NULL},{NULL,NULL,NULL,NULL,NULL,NULL}},
 	{{fk5q1n0s,fk5q1n1s,fk5q1n2s,fk5q1n3s,fk5q1n4s,fk5q1n5s},{fk5q2n0,fk5q2n1s,fk5q2n2s,fk5q2n3s,fk5q2n4s,fk5q2n5s}}
 };
+
 void lrhankel_recpart_fill(complex  double *target,
-               size_t maxn, size_t lrk_cutoff,
+               size_t maxp /*max. degree of transformed spherical Hankel fun, 
+	        also the max. order of the Hankel transform */, 
+	       size_t lrk_cutoff,
                complex double  *hct,
                unsigned kappa, double c, double k0, double k)
 {
-	memset(target, 0, (maxn+1)*sizeof(complex double));
+	assert(5 == kappa); // Only kappa == 5 implemented so far
-	complex double a[kappa+1], b[kappa+1], d[kappa+1], e[kappa+1];
+	assert(maxp <= 5); // only n <= implemented so far
+	// assert(lrk_cutoff <= TODO);
+	const lrhankelspec (*funarr)[MAXQM+1][MAXN+1] = (k>k0) ? transfuns_f : transfuns_n;
+	memset(target, 0, maxp*(maxp+1)/2*sizeof(complex double));
+	complex double a[kappa+1], b[kappa+1], d[kappa+1], e[kappa+1], ash[kappa+1];
 	for (size_t sigma = 0; sigma <= kappa; ++sigma) {
 		a[sigma] = (sigma * c - I * k0);
 		b[sigma] = csqrt(1+k*k/(a[sigma]*a[sigma]));
 		d[sigma] = 1/b[sigma];
 		e[sigma] = d[sigma] / a[sigma];
+		ash[sigma] = casinh(a[sigma]/k);
 	}
+	for (size_t ql = 0; (ql <= maxp) && (ql < lrk_cutoff); ++ql)  // ql is q-1, i.e. corresponds to the hankel term power
+		for (size_t deltam = 0; deltam <= maxp; ++deltam){
+			complex double result = funarr[kappa][ql][deltam](c,k0,k,a,b,d,e,ash);
+			for (size_t p = 0; p <= maxp; ++p)
+				trindex_cd(target,p)[deltam] += result * hankelcoeffs_get(hct,p)[ql];
+		}
 }
 
 #ifdef TESTING