./dipdip-dirty stuff moved to ./qpms + comments, asserts and consts

Former-commit-id: 6005fcf91ca007b7e6d7093d1e82a0a305f28ffd

dipdip-dirty/bessels.h

@@ -1,41 +0,0 @@
-#ifndef BESSELS_H
-#define BESSELS_H
-
-#include <stddef.h>
-#include <complex.h>
-
-complex double *hankelcoefftable_init(size_t maxn);
-
-
-
-static inline complex double *
-trindex_cd(complex double *arr, size_t n){
-  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 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 /*Must be of size maxn*(maxn+1)/2*/,
-	       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);
-
-#endif //BESSELS_H

qpms/bessels.c

@@ -39,7 +39,7 @@ complex double * hankelcoefftable_init(size_t maxn) {
 }
 
 void hankelparts_fill(complex double *lrt, complex double *srt, size_t maxn,
-								size_t lrk_cutoff, complex double *hct,
+								size_t lrk_cutoff, complex double const * const hct,
 								unsigned kappa, double c, double x) {
 	if (lrt) memset(lrt, 0, (maxn+1)*sizeof(complex double));
 	memset(srt, 0, (maxn+1)*sizeof(complex double));
@@ -48,7 +48,7 @@ void hankelparts_fill(complex double *lrt, complex double *srt, size_t maxn,
 	double xfrac = 1.; // x ** (-1-k)
 	for (size_t k = 0; k <= maxn; ++k) {
 	xfrac /= x;
-	  for(size_t n = k; n <= maxn; ++n) 
+	  for(size_t n = k; n <= maxn; ++n)  // TODO Kahan sums here
 	  	srt[n] += ((k<lrk_cutoff) ? antiregularisator : 1) 
 		  				* xfrac * hankelcoeffs_get(hct,n)[k];
 	  if (lrt && k < lrk_cutoff) for (size_t n = k; n <= maxn; ++n)

qpms/bessels.h

@@ -0,0 +1,54 @@
+#ifndef BESSELS_H
+#define BESSELS_H
+/* Short- and long-range parts of spherical Hankel functions
+ * and (cylindrical) Hankel transforms of the long-range parts.
+ * Currently, the implementation lies in bessels.c and 
+ * lrhankel_recspace_dirty.c. The latter contains the implementation
+ * of the Hankel transforms, but currenty only for a pretty limited
+ * set of parameters. The general implementation is a BIG TODO here.
+ */
+
+#include <stddef.h>
+#include <complex.h>
+
+complex double *hankelcoefftable_init(size_t maxn);
+
+
+// For navigating in the coefficients, maybe not for public use
+static inline complex double *
+trindex_cd(complex double const * const arr, size_t n){
+  return (complex double *)(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 const * const 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, /* terms e**(I x)/x**(k+1), 
+							       k>= longrange_order_cutoff go 
+							       completely to short-range part */
+		complex double const * const hankelcoefftable,
+		unsigned kappa, double vc, double x); // x = k0 * r
+
+
+
+/* Hankel transforms of the long-range parts of the spherical Hankel functions */
+// this declaration is general; however, the implementation 
+// is so far only for kappa == 5, maxp == 5 TODO
+void lrhankel_recpart_fill(complex  double *target_longrange_kspace /*Must be of size maxn*(maxn+1)/2*/,
+	       size_t maxp /* Max. degree of transformed spherical Hankel function,
+			      also the max. order of the Hankel transform */,
+	       size_t longrange_order_cutoff /* terms e**(I x)/x**(k+1), k>= longrange_order_cutoff go
+					        completely to the shortrange part
+					        index with hankelcoeffs_get(target,p)l[delta_m] */, 
+	       complex double const * const hankelcoefftable,
+	       unsigned kappa, double c, double k0, double k);
+
+#endif //BESSELS_H

qpms/lrhankel_recspace_dirty.c

@@ -12,6 +12,9 @@
  * 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.
+ *
+ * Therefore TODO use kahan summation.
+ *
  */
 
 #define MAXQM 1
@@ -232,12 +235,12 @@ void lrhankel_recpart_fill(complex  double *target,
                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,
+               complex double const *const hct,
                unsigned kappa, double c, double k0, double k)
 {
 	assert(5 == kappa); // Only kappa == 5 implemented so far
-	assert(maxp <= 5); // only n <= implemented so far
+	assert(maxp <= MAXN); // only n <= 5 implemented so far
-	// assert(lrk_cutoff <= TODO);
+	assert(lrk_cutoff <= MAXQM); // only q <= 2 implemented so far
 	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];