From c0f23fce558271867183538dc4702cf3a45792c3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Mon, 14 May 2018 06:52:32 +0300 Subject: [PATCH] ./dipdip-dirty stuff moved to ./qpms + comments, asserts and consts Former-commit-id: 6005fcf91ca007b7e6d7093d1e82a0a305f28ffd --- dipdip-dirty/bessels.h | 41 -------------- {dipdip-dirty => qpms}/bessels.c | 4 +- qpms/bessels.h | 54 +++++++++++++++++++ .../lrhankel_recspace_dirty.c | 9 ++-- 4 files changed, 62 insertions(+), 46 deletions(-) delete mode 100644 dipdip-dirty/bessels.h rename {dipdip-dirty => qpms}/bessels.c (93%) create mode 100644 qpms/bessels.h rename {dipdip-dirty => qpms}/lrhankel_recspace_dirty.c (97%) diff --git a/dipdip-dirty/bessels.h b/dipdip-dirty/bessels.h deleted file mode 100644 index 774b9c5..0000000 --- a/dipdip-dirty/bessels.h +++ /dev/null @@ -1,41 +0,0 @@ -#ifndef BESSELS_H -#define BESSELS_H - -#include -#include - -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 diff --git a/dipdip-dirty/bessels.c b/qpms/bessels.c similarity index 93% rename from dipdip-dirty/bessels.c rename to qpms/bessels.c index 3d8b778..e1aaf5e 100644 --- a/dipdip-dirty/bessels.c +++ b/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 +#include + +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 diff --git a/dipdip-dirty/lrhankel_recspace_dirty.c b/qpms/lrhankel_recspace_dirty.c similarity index 97% rename from dipdip-dirty/lrhankel_recspace_dirty.c rename to qpms/lrhankel_recspace_dirty.c index d286e54..47d2b9f 100644 --- a/dipdip-dirty/lrhankel_recspace_dirty.c +++ b/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(lrk_cutoff <= TODO); + assert(maxp <= MAXN); // only n <= 5 implemented so far + 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];