diff --git a/dipdip-dirty/bessels.c b/dipdip-dirty/bessels.c index cd49680..80df3ff 100644 --- a/dipdip-dirty/bessels.c +++ b/dipdip-dirty/bessels.c @@ -1,19 +1,19 @@ #include "bessels.h" #include #include +#include +#include static const double ln2 = 0.69314718055994531; -#include -// 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 +#include +#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 +