From 85b881b1d6c3560e6741f8fb0ca2b3059200ecee Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Thu, 18 Jan 2018 17:55:44 +0000 Subject: [PATCH] [dirty dipoles] few long-range hankel transforms rewriten from mathematica Former-commit-id: 2554e2f670bd9e405210d3e2f62596b2295b0b7d --- dipdip-dirty/lrhankel_recspace_dirty.c | 123 +++++++++++++++++++++++++ 1 file changed, 123 insertions(+) create mode 100644 dipdip-dirty/lrhankel_recspace_dirty.c diff --git a/dipdip-dirty/lrhankel_recspace_dirty.c b/dipdip-dirty/lrhankel_recspace_dirty.c new file mode 100644 index 0000000..f1b2644 --- /dev/null +++ b/dipdip-dirty/lrhankel_recspace_dirty.c @@ -0,0 +1,123 @@ +#include "bessels.h" +//#include "mdefs.h" +#include +#include +#define SQ(x) ((x)*(x)) + +#define MAXQM 1 +#define MAXN 2 +#define MAXKAPPA 5 + +typedef complex double (*lrhankelspec)(double, double, double, + const complex double *, + const complex double *, + const complex double *, + const complex double *); +// complex double fun(double c, double k0, double k, ccd *a, ccd *b, ccd *d, ccd *e) + + +complex double fk5q1n0l(double c, double k0, double k, + const complex double *a, const complex double *b, const complex double *d, const complex double *e) { + return (e[0]-5*e[1]+10*e[2]-10*e[3]+5*e[4]-e[5])/k0; +} +complex double fk5q1n1l(double c, double k0, double k, + const complex double *a, const complex double *b, const complex double *d, const complex double *e) { + return (-d[0]+5*d[1]-10*d[2]+10*d[3]-5*d[4]+d[5])/(k0*k); +} +complex double fk5q1n2l(double c, double k0, double k, + const complex double *a, const complex double *b, const complex double *d, const complex double *e) { + double t = 2/(k*k); + return ( (e[0] - t*a[0] + t*d[0]*a[0]) + -5 * (e[1] - t*a[1] + t*d[1]*a[1]) + +10 *(e[2] - t*a[2] + t*d[2]*a[2]) + -10 *(e[3] - t*a[3] + t*d[3]*a[3]) + +5 * (e[4] - t*a[4] + t*d[4]*a[4]) + - (e[5] - t*a[5] + t*d[5]*a[5]) + )/k0; +} +complex double fk5q2n0l(double c, double k0, double k, + const complex double *a, const complex double *b, const complex double *d, const complex double *e) { + return 0; // FIXME +} +complex double fk5q2n1l(double c, double k0, double k, + const complex double *a, const complex double *b, const complex double *d, const complex double *e) { + return ( b[0]*a[0] + - 5 *b[1]*a[1] + +10 *b[2]*a[2] + -10 *b[3]*a[3] + + 5 *b[4]*a[4] + - b[5]*a[5] + )/(k*k0*k0); +} +complex double fk5q2n2l(double c, double k0, double k, + const complex double *a, const complex double *b, const complex double *d, const complex double *e) { + return ( b[0]*a[0]*a[0] + + 5 * b[1]*a[1]*a[1] + -10 * b[2]*a[2]*a[2] + +10 * b[3]*a[3]*a[3] + - 5 * b[4]*a[4]*a[4] + + b[5]*a[5]*a[5] + ) / (k*k*k0*k0); +} + +#if 0 +complex double fk5q1n0s(double c, double k0, double k, + const complex double *a, const complex double *b, const complex double *d, const complex double *e) { + return (p[0]-5*p[1]+10*p[2]-10*p[3]+5*p[4]-p[5])/k0; +} +complex double fk5q1n1s(double c, double k0, double k, + const complex double *a, const complex double *b, const complex double *d, const complex double *e) { + return ; +} +complex double fk5q1n2s(double c, double k0, double k, + const complex double *a, const complex double *b, const complex double *d, const complex double *e) { + return ; +} +complex double fk5q2n0s(double c, double k0, double k, + const complex double *a, const complex double *b, const complex double *d, const complex double *e) { + return ; +} +complex double fk5q2n1s(double c, double k0, double k, + const complex double *a, const complex double *b, const complex double *d, const complex double *e) { + return ; +} +complex double fk5q2n2s(double c, double k0, double k, + const complex double *a, const complex double *b, const complex double *d, const complex double *e) { + return ; +} + +static 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}, + {NULL,NULL,NULL},{NULL,NULL,NULL},{NULL,NULL,NULL}, + {TODO,TODO,TODO},{TODO,TODO,TODO},{TODO,TODO,TODO} +}; +#endif + +static 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}}, + {{NULL,NULL,NULL},{NULL,NULL,NULL}}, + {{fk5q1n0l,fk5q1n1l,fk5q1n2l},{fk5q2n0l/*FIXME*/,fk5q2n1l,fk5q2n2l}} +}; + +void lrhankel_recpart_fill(complex double *target, + size_t maxn, size_t lrk_cutoff, + complex double *hct, + unsigned kappa, double c, double k0, double k) +{ + memset(target, 0, (maxn+1)*sizeof(complex double)); + complex double a[kappa+1], b[kappa+1], d[kappa+1], e[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]; + } +} + +