diff --git a/dipdip-dirty/lrhankel_recspace_dirty.c b/dipdip-dirty/lrhankel_recspace_dirty.c index f1b2644..ac06bdb 100644 --- a/dipdip-dirty/lrhankel_recspace_dirty.c +++ b/dipdip-dirty/lrhankel_recspace_dirty.c @@ -18,14 +18,17 @@ typedef complex double (*lrhankelspec)(double, double, double, 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) { + // FIXME 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) { + // FIXME 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) { + const complex double *a, const complex double *b, const complex double *d, const complex double *e) { double t = 2/(k*k); + // FIXME 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]) @@ -41,6 +44,7 @@ complex double fk5q2n0l(double c, double k0, double k, } 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) { + // FIXME return ( b[0]*a[0] - 5 *b[1]*a[1] +10 *b[2]*a[2] @@ -60,6 +64,51 @@ complex double fk5q2n2l(double c, double k0, double k, ) / (k*k*k0*k0); } +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 (e[0]-5*e[1]+10*e[2]-10*e[3]+5*e[4]-e[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 (-d[0]+5*d[1]-10*d[2]+10*d[3]-5*d[4]+d[5])/(k0*k); +} +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) { + 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 fk5q2n1s(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 fk5q2n2s(double c, double k0, double k, + const complex double *a, const complex double *b, const complex double *d, const complex double *e) { + // FIXME + 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) { @@ -120,4 +169,33 @@ void lrhankel_recpart_fill(complex double *target, } } +#include +int main() { + double k0 = 0.7; + double c = 0.1324; + double kmin = 0.000; + double kmax = 20; + double kstep = 0.001; + size_t kappa = 5; + for (double k = kmin; k <= kmax; k += kstep) { + printf("%f ", k); + 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]; + } + for (size_t qm = 0; qm <= MAXQM; ++qm) + for (size_t n = 0; n <= MAXN; ++n) + if (!((qm==1)&&(n==0))){ // skip q==2, n=0 for now +// complex double fun(double c, double k0, double k, ccd *a, ccd *b, ccd *d, ccd *e) + complex double result = + transfuns_f[kappa][qm][n](c,k0,k,a,b,d,e); + printf("%.16e %.16e ", creal(result), cimag(result)); + } + printf("\n"); + } + return 0; +}