[dirty dipoles] regularised hankel functions

Former-commit-id: 58b4d91b718a0d864555961bb199761d1774313d

dipdip-dirty/bessels.c

@@ -1,19 +1,19 @@
 #include "bessels.h"
 #include <stdlib.h>
 #include <math.h>
+#include <stdio.h>
+#include <string.h>
 
 static const double ln2 = 0.69314718055994531;
 
-#include <stdio.h>
 
-// general; gives an array of size 
+// general; gives an array of size xxx with TODODESC
 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;
@@ -38,3 +38,57 @@ complex double * hankelcoefftable_init(size_t maxn) {
 	return hct;
 }
 
+void hankelLR2_fill(complex double *target, size_t maxn, complex double *hct, 
+								unsigned kappa, double c, double x) {
+	memset(target, 0, (maxn+1)*sizeof(complex double));
+	double regularisator = pow(1. - exp(-c * x), (double) kappa);
+	complex double expix = cexp(I * x);
+	double xfrac = 1.; // x ** (-1-k)
+	for (size_t k = 0; k < 2; ++k) {
+					xfrac /= x;
+					for(size_t n = k; n <= maxn; ++n) 
+						target[n] += regularisator * xfrac * hankelcoeffs_get(hct,n)[k];
+	}
+	for(size_t n = 0; n <= maxn; ++n)
+					target[n] *= expix;
+}
+
+void hankelSR2_fill(complex double *target, size_t maxn, complex double *hct, 
+								unsigned kappa, double c, double x) {
+	memset(target, 0, (maxn+1)*sizeof(complex double));
+	double antiregularisator = 1 - pow(1. - exp(-c * x), (double) kappa);
+	complex double expix = cexp(I * x);
+	double xfrac = 1.; // x ** (-1-k)
+	for (size_t k = 0; k < maxn; ++k) {
+					xfrac /= x;
+					for(size_t n = k; n <= maxn; ++n) 
+						target[n] += ((k<2) ? antiregularisator : 1) 
+										* xfrac * hankelcoeffs_get(hct,n)[k];
+	}
+	for(size_t n = 0; n <= maxn; ++n)
+					target[n] *= expix;
+}
+
+void hankelparts_fill(complex double *lrt, complex double *srt, size_t maxn,
+								size_t lrk_cutoff, complex double *hct,
+								unsigned kappa, double c, double x) {
+	memset(lrt, 0, (maxn+1)*sizeof(complex double));
+	memset(srt, 0, (maxn+1)*sizeof(complex double));
+	double regularisator = pow(1. - exp(-c * x), (double) kappa);
+	double antiregularisator = 1. - regularisator;
+	double xfrac = 1.; // x ** (-1-k)
+	for (size_t k = 0; k <= maxn; ++k) {
+	xfrac /= x;
+	  for(size_t n = k; n <= maxn; ++n) 
+	  	srt[n] += ((k<lrk_cutoff) ? antiregularisator : 1) 
+		  				* xfrac * hankelcoeffs_get(hct,n)[k];
+	  if (k < lrk_cutoff) for (size_t n = k; n <= maxn; ++n)
+	  	lrt[n] += regularisator * xfrac * hankelcoeffs_get(hct,n)[k];
+	}
+
+	complex double expix = cexp(I * x);
+	for(size_t n = 0; n <= maxn; ++n){
+					srt[n] *= expix;
+					lrt[n] *= expix;
+	}
+}

dipdip-dirty/bessels.h

@@ -11,8 +11,14 @@ complex double *hankelcoefftable_init(size_t maxn);
 // 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;
+  return hankelcoefftable + n*(n+1)/2;
 }
 
 
+// general; target_longrange and target_shortrange are of size (maxn+1)
+void hankelparts_fill(complex double *target_longrange, complex double *target_shortrange,
+		size_t maxn, size_t longrange_k_cutoff, // x**(-(k+1)-1) terms go completely to short-range part
+		complex double *hankelcoefftable,
+		unsigned kappa, double c, double x);
+
 #endif //BESSELS_H

dipdip-dirty/testmain.c

@@ -0,0 +1,47 @@
+#include "bessels.h"
+#include <stdio.h>
+#include <stdlib.h>
+#if 0
+int main() {
+	size_t maxn = 5;
+	complex double *hct = hankelcoefftable_init(maxn);
+	for (size_t n = 0; n <= maxn; ++n) {
+		printf("n = %zd\n", n);
+		for(size_t k = 0; k<=n; ++k)
+			printf("%p: %f + %fj,\n", hankelcoeffs_get(hct,n) + k, creal(hankelcoeffs_get(hct,n)[k]),
+					cimag(hankelcoeffs_get(hct,n)[k]));
+		printf("\n");
+	}
+	printf("%f+%fj\n",creal(cpow(I,(ptrdiff_t)-1)), cimag(cpow(I,(ptrdiff_t)-1)));
+	return 0;
+}
+#endif
+
+//#if 0
+int main() {
+	size_t maxn = 6;
+	size_t lrk_cutoff = 2;
+	size_t kappa = 4;
+	double c = 0.1324;
+
+	double xmin = 0.;
+	double xstep = 0.001;
+	double xmax = 20;
+
+	complex double *hct = hankelcoefftable_init(maxn);
+	complex double srhankel[maxn+1];
+	complex double lrhankel[maxn+1];
+
+	for(double x = xmin; x <= xmax; x += xstep) {
+		hankelparts_fill(lrhankel, srhankel, maxn, lrk_cutoff, hct, kappa, c, x);
+		printf("%f ", x);
+		for(size_t n = 0; n <= maxn; ++n)
+			printf("%.16e %.16e %.16e %.16e ", creal(lrhankel[n]), cimag(lrhankel[n]), creal(srhankel[n]), cimag(srhankel[n]));
+		printf("\n");
+	}
+	
+	free(hct);
+	return 0;
+}
+//#endif
+