diff --git a/qpms/test_translations.c b/qpms/test_translations.c index 1ce2c29..03ed04e 100644 --- a/qpms/test_translations.c +++ b/qpms/test_translations.c @@ -4,6 +4,7 @@ #include typedef struct { + qpms_normalisation_t norm; int m, n, mu, nu; sph_t kdlj; qpms_bessel_t J; diff --git a/qpms/test_translations2.c b/qpms/test_translations2.c index 5d24293..e0c67e4 100644 --- a/qpms/test_translations2.c +++ b/qpms/test_translations2.c @@ -4,6 +4,7 @@ #include typedef struct { + qpms_normalisation_t norm; int m, n, mu, nu; sph_t kdlj; qpms_bessel_t J; @@ -17,14 +18,19 @@ testcase_single_trans_t testcases_Taylor[] = { int lMax=10; int main() { - qpms_trans_calculator *c = qpms_trans_calculator_init(lMax, QPMS_NORMALIZATION_TAYLOR); + qpms_trans_calculator *c = qpms_trans_calculator_init(lMax, QPMS_NORMALISATION_TAYLOR); for(testcase_single_trans_t *tc = testcases_Taylor; tc->J != QPMS_BESSEL_UNDEF; tc++) { if (!tc->n || !tc->nu || tc->n > lMax || tc->nu > lMax ) continue; printf("m=%d, n=%d, mu=%d, nu=%d,\n", tc->m,tc->n,tc->mu,tc->nu); +#if 0 complex double A = qpms_trans_single_A_Taylor(tc->m, tc->n, tc->mu, tc->nu, tc->kdlj, true, tc->J); complex double B = qpms_trans_single_B_Taylor(tc->m, tc->n, tc->mu, tc->nu, tc->kdlj, true, tc->J); +#else + complex double A = qpms_trans_single_A(QPMS_NORMALISATION_TAYLOR, tc->m, tc->n, tc->mu, tc->nu, tc->kdlj, true, tc->J); + complex double B = qpms_trans_single_B(QPMS_NORMALISATION_TAYLOR,tc->m, tc->n, tc->mu, tc->nu, tc->kdlj, true, tc->J); +#endif complex double A2 = qpms_trans_calculator_get_A(c, tc->m, tc->n, tc->mu, tc->nu, tc->kdlj, true, tc->J); complex double B2 = qpms_trans_calculator_get_B(c, tc->m, tc->n, tc->mu, tc->nu, tc->kdlj, true, tc->J); printf("A = %.16f+%.16fj, relerr=%.16f, J=%d\n", diff --git a/qpms/test_translations_xu.c b/qpms/test_translations_xu.c new file mode 100644 index 0000000..f23d8c1 --- /dev/null +++ b/qpms/test_translations_xu.c @@ -0,0 +1,61 @@ +#include "translations.h" +#include +//#include +#include + +typedef struct { + qpms_normalisation_t norm; + int m, n, mu, nu; + sph_t kdlj; + qpms_bessel_t J; + complex double result_A, result_B; +} testcase_single_trans_t; + +testcase_single_trans_t testcases_xu[] = { +#include "testcases_translations_Xu" +}; + +int lMax=10; + +int main() { + qpms_trans_calculator *c = qpms_trans_calculator_init(lMax, QPMS_NORMALISATION_XU); + + for(testcase_single_trans_t *tc = testcases_xu; tc->J != QPMS_BESSEL_UNDEF; tc++) { + if (!tc->n || !tc->nu || tc->n > lMax || tc->nu > lMax ) continue; + + printf("m=%d, n=%d, mu=%d, nu=%d,\n", tc->m,tc->n,tc->mu,tc->nu); + complex double A = qpms_trans_single_A(QPMS_NORMALISATION_XU,tc->m, tc->n, tc->mu, tc->nu, tc->kdlj, true, tc->J); + complex double B = qpms_trans_single_B(QPMS_NORMALISATION_XU,tc->m, tc->n, tc->mu, tc->nu, tc->kdlj, true, tc->J); + complex double A2 = qpms_trans_calculator_get_A(c, tc->m, tc->n, tc->mu, tc->nu, tc->kdlj, true, tc->J); + complex double B2 = qpms_trans_calculator_get_B(c, tc->m, tc->n, tc->mu, tc->nu, tc->kdlj, true, tc->J); + printf("A = %.16f+%.16fj, relerr=%.16f, J=%d\n", + creal(A), cimag(A), (0 == cabs(tc->result_A - A)) ? 0 : + cabs(tc->result_A - A)/((cabs(A) < cabs(tc->result_A)) ? cabs(A) : cabs(tc->result_A)), + tc->J); + printf("A' = %.16f+%.16fj, relerr=%.16f, relerr2=%.3e\n", + creal(A2), cimag(A2), (0 == cabs(tc->result_A - A2)) ? 0 : + cabs(tc->result_A - A2)/((cabs(A2) < cabs(tc->result_A)) ? cabs(A2) : cabs(tc->result_A)), + (0 == cabs(A - A2)) ? 0 : + cabs(A - A2)/((cabs(A2) < cabs(A)) ? cabs(A2) : cabs(A)) + ); + printf("B = %.16f+%.16fj, relerr=%.16f, J=%d\n", + creal(B), cimag(B), (0 == cabs(tc->result_B - B)) ? 0 : + cabs(tc->result_B - B)/((cabs(B) < cabs(tc->result_B)) ? cabs(B) : cabs(tc->result_B)), + tc->J); + printf("B' = %.16f+%.16fj, relerr=%.16f, relerr2=%.3e\n", + creal(B2), cimag(B2), (0 == cabs(tc->result_B - B2)) ? 0 : + cabs(tc->result_B - B2)/((cabs(B2) < cabs(tc->result_B)) ? cabs(B2) : cabs(tc->result_B)), + (0 == cabs(B - B2)) ? 0 : + cabs(B - B2)/((cabs(B2) < cabs(B)) ? cabs(B2) : cabs(B)) + ); + } + complex double A,B; + // Test of zero R + sph_t kdlj = {0, 1, 2}; + int m = -1, n = 1, mu = -1, nu = 1; + qpms_trans_calculator_get_AB_p(c,&A,&B,m,n,mu,nu,kdlj,false,3); + printf("A = %.6e+%.6ej, B = %.6e+%.6ej\n", creal(A),cimag(A),creal(B),cimag(B)); + qpms_trans_calculator_free(c); +} + + diff --git a/qpms/testcases_taylor.REMOVED.git-id b/qpms/testcases_taylor.REMOVED.git-id index ba0ba20..a714f6b 100644 --- a/qpms/testcases_taylor.REMOVED.git-id +++ b/qpms/testcases_taylor.REMOVED.git-id @@ -1 +1 @@ -2b61a18b6426306ba80a50ee30339651313fc44f \ No newline at end of file +7f001489805a48cf44317137517740bc40512657 \ No newline at end of file diff --git a/qpms/testcases_translations_Xu b/qpms/testcases_translations_Xu new file mode 100644 index 0000000..1923bcf --- /dev/null +++ b/qpms/testcases_translations_Xu @@ -0,0 +1,33 @@ +//{norm, m,n,mu,nu, {kd.r, kd.theta, kd.phi}, J, areal + aimag * I, breal + bimag* I}, +// Yu-lin Xu, Journal of computational physics 127, 285–298 (1996), Table II +{QPMS_NORMALISATION_XU, 1,2,1,2, {2., 0.5, 0.5}, QPMS_HANKEL_PLUS, .40952764e0 + -.37687246e-1 * I, -.95439506e0 + .14571311e0* I}, +{QPMS_NORMALISATION_XU, 1,5,2,6, {2., 0.5, 0.5}, QPMS_HANKEL_PLUS, .42246528e6 + .23079387e6 * I, .24466070e5 + -.44784830e5* I}, +{QPMS_NORMALISATION_XU, -2,10,7,15, {2., 0.5, 0.5}, QPMS_HANKEL_PLUS, .13425435e14 + .62258200e14 * I, .20835862e13 + -.44930709e12* I}, +{QPMS_NORMALISATION_XU, 10,20,3,22, {2., 0.5, 0.5}, QPMS_HANKEL_PLUS, .12643610e58 + .33753591e58 * I, .18943589e57 + -.70959965e56* I}, +{QPMS_NORMALISATION_XU, -28,30,29,32, {2., 0.5, 0.5}, QPMS_HANKEL_PLUS, -.13928339e-10 + .60669708-10 * I, -.31805959e-13 + -.73019004e-14* I}, +{QPMS_NORMALISATION_XU, 22,40,-39,45, {2., 0.5, 0.5}, QPMS_HANKEL_PLUS, .34997554e216 + .45574651e216 * I, -.23949371e214 + .18391132e214* I}, +{QPMS_NORMALISATION_XU, 5,50,-22,61, {2., 0.5, 0.5}, QPMS_HANKEL_PLUS, .89567769e226 + -.12101307e227 * I, .50442227e224 + .37334790e224* I}, +// Yu-lin Xu, Journal of computational physics 139, 137–165 (1998), Table 12 +{QPMS_NORMALISATION_XU, 8,10,-9,12, {2., 0.5, 0.5}, QPMS_HANKEL_PLUS, .3663964990e+35 + -.2762412192e+35 * I, -.8370892023e+32 + -.1110285257e+32* I}, +{QPMS_NORMALISATION_XU, 0,10,0,10, {2., 0.5, 0.5}, QPMS_HANKEL_PLUS, .2969682019e+00 + -.1928601440e+18 * I, .0000000000e+00 + .0000000000e+00* I}, +{QPMS_NORMALISATION_XU, -2,11,3,9, {2., 0.5, 0.5}, QPMS_HANKEL_PLUS, .7726121583e+12 + .1034255820e+13 * I, .1222239141e+11 + -.9130398908e+10* I}, +{QPMS_NORMALISATION_XU, -12,13,10,15, {2., 0.5, 0.5}, QPMS_HANKEL_PLUS, .3290937356e+01 + .1456483748e-01 * I, -.1763167849e-03 + .3983892680e-01* I}, +{QPMS_NORMALISATION_XU, -15,16,17,18, {2., 0.5, 0.5}, QPMS_HANKEL_PLUS, .3793897303e-08 + -.1261972860e-07 * I, -.3042702016e-11 + -.9147343290e-12* I}, +{QPMS_NORMALISATION_XU, -5,20,5,20, {2., 0.5, 0.5}, QPMS_HANKEL_PLUS, .4040625669e+34 + -.1195269260e+34 * I, .0000000000e+00 + .0000000000e+00* I}, +{QPMS_NORMALISATION_XU, 10,18,15,22, {2., 0.5, 0.5}, QPMS_HANKEL_PLUS, -.6206840651e+36 + -.8308775621e+36 * I, -.3610252125e+35 + .2696938836e+35* I}, +{QPMS_NORMALISATION_XU, 10,30,-10,30, {2., 0.5, 0.5}, QPMS_HANKEL_PLUS, .1807705110e+110 + .2788115866e+110 * I, .0000000000e+00 + .0000000000e+00* I}, +{QPMS_NORMALISATION_XU, 18,33,20,38, {2., 0.5, 0.5}, QPMS_HANKEL_PLUS, .3343492687e+92 + .5207181338e+92 * I, .1759309957e+91 + -.1129639932e+91* I}, +{QPMS_NORMALISATION_XU, -35,36,11,12, {2., 0.5, 0.5}, QPMS_HANKEL_PLUS, -.1901528547e-15 + .1197320691e-15 * I, -.1618572254e-18 + -.2570540515e-18* I}, +{QPMS_NORMALISATION_XU, 36,36,-38,38, {2., 0.5, 0.5}, QPMS_HANKEL_PLUS, .4146334728e+191 + -.4931584782e+191 * I, .0000000000e+00 + .0000000000e+00* I}, +{QPMS_NORMALISATION_XU, -35,40,35,40, {2., 0.5, 0.5}, QPMS_HANKEL_PLUS, -.6514262216e-05 + .1374854333e-04 * I, .0000000000e+00 + .0000000000e+00* I}, +{QPMS_NORMALISATION_XU, 32,35,-43,45, {2., 0.5, 0.5}, QPMS_HANKEL_PLUS, .2762232925e+212 + -.1368895313e+213 * I, .8373862584e+209 + .1689724460e+209* I}, +{QPMS_NORMALISATION_XU, 38,42,-39,45, {2., 0.5, 0.5}, QPMS_HANKEL_PLUS, -.2298689786e+235 + .2371029493e+235 * I, .1277697908e+232 + .1238711556e+232* I}, +{QPMS_NORMALISATION_XU, -42,42,45,45, {2., 0.5, 0.5}, QPMS_HANKEL_PLUS, .3488835702e-28 + -.1826524477e-28 * I, .0000000000e+00 + .0000000000e+00* I}, +{QPMS_NORMALISATION_XU, -43,45,41,42, {2., 0.5, 0.5}, QPMS_HANKEL_PLUS, .5178100899e-22 + .1186503822e-21 * I, .3274958627e-25 + -.1429246656e-25* I}, +{QPMS_NORMALISATION_XU, 48,50,-30,49, {2., 0.5, 0.5}, QPMS_HANKEL_PLUS, .3393827523e+267 + -.1226717423e+268 * I, -.5718637033e+265 + -.1582113973e+265* I}, +{QPMS_NORMALISATION_XU, -72,72,1,3, {2., 0.5, 0.5}, QPMS_HANKEL_PLUS, .6946365327e-42 + -.1782022552e-41 * I, .1833377882e-43 + .7146549596e-44* I}, +{QPMS_NORMALISATION_XU, 42,52,9,81, {2., 0.5, 0.5}, QPMS_HANKEL_PLUS, .3656934399e+271 + .3705813223e+271 * I, -.4499925012e+269 + -.4440572037e+269* I}, +{QPMS_NORMALISATION_XU, 18,100,-5,45, {2., 0.5, 0.5}, QPMS_HANKEL_PLUS, .4118769973e+293 + .7460688240e+293 * I, .5914795871e+291 + -.3265339985e+291* I}, +// THE END +{QPMS_NORMALISATION_UNDEF, 0,0,0,0, {0, 0, 0}, QPMS_BESSEL_UNDEF, 0 + 0 * I, 0 + 0* I} + diff --git a/qpms/testcases_translations_Xu.fods b/qpms/testcases_translations_Xu.fods new file mode 100644 index 0000000..2f94c40 --- /dev/null +++ b/qpms/testcases_translations_Xu.fods @@ -0,0 +1,2759 @@ + + + + 2018-02-07T18:03:11.6448020812018-02-07T19:45:24.778557362PT29M30S6LibreOffice/5.1.6.2$Linux_X86_64 LibreOffice_project/10m0$Build-2 + + + 0 + 0 + 36119 + 14449 + + + view1 + + + 10 + 9 + 0 + 0 + 0 + 0 + 2 + 0 + 0 + 0 + 0 + 0 + 100 + 60 + true + + + Taulukko1 + 1853 + 0 + 100 + 60 + false + true + true + true + 12632256 + true + true + true + true + false + false + false + 1000 + 1000 + 1 + 1 + true + + + + + false + false + true + true + false + false + false + 1000 + 1000 + true + true + true + true + true + true + true + false + 12632256 + false + + 3 + 1 + true + 1 + true + + true + 0 + + + + + + + + + + + + + + + + + + + + + + + + + + + € + + + + - + + + € + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + ??? + + + + Sivu 1 + + + + + + + ??? (???) + + + 00.00.0000, 00:00:00 + + + + + Sivu 1 / 99 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + // + + + { + + + norm + + + , + + + m + + + , + + + n + + + , + + + mu + + + , + + + nu + + + , { + + + kd.r + + + , + + + kd.theta + + + , + + + kd.phi + + + }, + + + J + + + , + + + areal + + + + + + + aimag + + + * I, + + + breal + + + + + + + bimag + + + * I}, + + + + + // Yu-lin Xu, Journal of computational physics 127, 285–298 (1996), Table II + + + + + + + + + + + + + + + { + + + QPMS_NORMALISATION_XU + + + , + + + 1 + + + , + + + 2 + + + , + + + 1 + + + , + + + 2 + + + , { + + + 2. + + + , + + + 0.5 + + + , + + + 0.5 + + + }, + + + QPMS_HANKEL_PLUS + + + , + + + .40952764e0 + + + + + + + -.37687246e-1 + + + * I, + + + -.95439506e0 + + + + + + + .14571311e0 + + + * I}, + + + + + + { + + + QPMS_NORMALISATION_XU + + + , + + + 1 + + + , + + + 5 + + + , + + + 2 + + + , + + + 6 + + + , { + + + 2. + + + , + + + 0.5 + + + , + + + 0.5 + + + }, + + + QPMS_HANKEL_PLUS + + + , + + + .42246528e6 + + + + + + + .23079387e6 + + + * I, + + + .24466070e5 + + + + + + + -.44784830e5 + + + * I}, + + + + + + { + + + QPMS_NORMALISATION_XU + + + , + + + -2 + + + , + + + 10 + + + , + + + 7 + + + , + + + 15 + + + , { + + + 2. + + + , + + + 0.5 + + + , + + + 0.5 + + + }, + + + QPMS_HANKEL_PLUS + + + , + + + .13425435e14 + + + + + + + .62258200e14 + + + * I, + + + .20835862e13 + + + + + + + -.44930709e12 + + + * I}, + + + + + + { + + + QPMS_NORMALISATION_XU + + + , + + + 10 + + + , + + + 20 + + + , + + + 3 + + + , + + + 22 + + + , { + + + 2. + + + , + + + 0.5 + + + , + + + 0.5 + + + }, + + + QPMS_HANKEL_PLUS + + + , + + + .12643610e58 + + + + + + + .33753591e58 + + + * I, + + + .18943589e57 + + + + + + + -.70959965e56 + + + * I}, + + + + + + { + + + QPMS_NORMALISATION_XU + + + , + + + -28 + + + , + + + 30 + + + , + + + 29 + + + , + + + 32 + + + , { + + + 2. + + + , + + + 0.5 + + + , + + + 0.5 + + + }, + + + QPMS_HANKEL_PLUS + + + , + + + -.13928339e-10 + + + + + + + .60669708-10 + + + * I, + + + -.31805959e-13 + + + + + + + -.73019004e-14 + + + * I}, + + + + + + { + + + QPMS_NORMALISATION_XU + + + , + + + 22 + + + , + + + 40 + + + , + + + -39 + + + , + + + 45 + + + , { + + + 2. + + + , + + + 0.5 + + + , + + + 0.5 + + + }, + + + QPMS_HANKEL_PLUS + + + , + + + .34997554e216 + + + + + + + .45574651e216 + + + * I, + + + -.23949371e214 + + + + + + + .18391132e214 + + + * I}, + + + + + + { + + + QPMS_NORMALISATION_XU + + + , + + + 5 + + + , + + + 50 + + + , + + + -22 + + + , + + + 61 + + + , { + + + 2. + + + , + + + 0.5 + + + , + + + 0.5 + + + }, + + + QPMS_HANKEL_PLUS + + + , + + + .89567769e226 + + + + + + + -.12101307e227 + + + * I, + + + .50442227e224 + + + + + + + .37334790e224 + + + * I}, + + + + + // Yu-lin Xu, Journal of computational physics 139, 137–165 (1998), Table 12 + + + + + + + + + + + + + + + { + + + QPMS_NORMALISATION_XU + + + , + + + 8 + + + , + + + 10 + + + , + + + -9 + + + , + + + 12 + + + , { + + + 2. + + + , + + + 0.5 + + + , + + + 0.5 + + + }, + + + QPMS_HANKEL_PLUS + + + , + + + .3663964990e+35 + + + + + + + -.2762412192e+35 + + + * I, + + + -.8370892023e+32 + + + + + + + -.1110285257e+32 + + + * I}, + + + + + + { + + + QPMS_NORMALISATION_XU + + + , + + + 0 + + + , + + + 10 + + + , + + + 0 + + + , + + + 10 + + + , { + + + 2. + + + , + + + 0.5 + + + , + + + 0.5 + + + }, + + + QPMS_HANKEL_PLUS + + + , + + + .2969682019e+00 + + + + + + + -.1928601440e+18 + + + * I, + + + .0000000000e+00 + + + + + + + .0000000000e+00 + + + * I}, + + + + + + { + + + QPMS_NORMALISATION_XU + + + , + + + -2 + + + , + + + 11 + + + , + + + 3 + + + , + + + 9 + + + , { + + + 2. + + + , + + + 0.5 + + + , + + + 0.5 + + + }, + + + QPMS_HANKEL_PLUS + + + , + + + .7726121583e+12 + + + + + + + .1034255820e+13 + + + * I, + + + .1222239141e+11 + + + + + + + -.9130398908e+10 + + + * I}, + + + + + + { + + + QPMS_NORMALISATION_XU + + + , + + + -12 + + + , + + + 13 + + + , + + + 10 + + + , + + + 15 + + + , { + + + 2. + + + , + + + 0.5 + + + , + + + 0.5 + + + }, + + + QPMS_HANKEL_PLUS + + + , + + + .3290937356e+01 + + + + + + + .1456483748e-01 + + + * I, + + + -.1763167849e-03 + + + + + + + .3983892680e-01 + + + * I}, + + + + + + { + + + QPMS_NORMALISATION_XU + + + , + + + -15 + + + , + + + 16 + + + , + + + 17 + + + , + + + 18 + + + , { + + + 2. + + + , + + + 0.5 + + + , + + + 0.5 + + + }, + + + QPMS_HANKEL_PLUS + + + , + + + .3793897303e-08 + + + + + + + -.1261972860e-07 + + + * I, + + + -.3042702016e-11 + + + + + + + -.9147343290e-12 + + + * I}, + + + + + + { + + + QPMS_NORMALISATION_XU + + + , + + + -5 + + + , + + + 20 + + + , + + + 5 + + + , + + + 20 + + + , { + + + 2. + + + , + + + 0.5 + + + , + + + 0.5 + + + }, + + + QPMS_HANKEL_PLUS + + + , + + + .4040625669e+34 + + + + + + + -.1195269260e+34 + + + * I, + + + .0000000000e+00 + + + + + + + .0000000000e+00 + + + * I}, + + + + + + { + + + QPMS_NORMALISATION_XU + + + , + + + 10 + + + , + + + 18 + + + , + + + 15 + + + , + + + 22 + + + , { + + + 2. + + + , + + + 0.5 + + + , + + + 0.5 + + + }, + + + QPMS_HANKEL_PLUS + + + , + + + -.6206840651e+36 + + + + + + + -.8308775621e+36 + + + * I, + + + -.3610252125e+35 + + + + + + + .2696938836e+35 + + + * I}, + + + + + + { + + + QPMS_NORMALISATION_XU + + + , + + + 10 + + + , + + + 30 + + + , + + + -10 + + + , + + + 30 + + + , { + + + 2. + + + , + + + 0.5 + + + , + + + 0.5 + + + }, + + + QPMS_HANKEL_PLUS + + + , + + + .1807705110e+110 + + + + + + + .2788115866e+110 + + + * I, + + + .0000000000e+00 + + + + + + + .0000000000e+00 + + + * I}, + + + + + + { + + + QPMS_NORMALISATION_XU + + + , + + + 18 + + + , + + + 33 + + + , + + + 20 + + + , + + + 38 + + + , { + + + 2. + + + , + + + 0.5 + + + , + + + 0.5 + + + }, + + + QPMS_HANKEL_PLUS + + + , + + + .3343492687e+92 + + + + + + + .5207181338e+92 + + + * I, + + + .1759309957e+91 + + + + + + + -.1129639932e+91 + + + * I}, + + + + + + { + + + QPMS_NORMALISATION_XU + + + , + + + -35 + + + , + + + 36 + + + , + + + 11 + + + , + + + 12 + + + , { + + + 2. + + + , + + + 0.5 + + + , + + + 0.5 + + + }, + + + QPMS_HANKEL_PLUS + + + , + + + -.1901528547e-15 + + + + + + + .1197320691e-15 + + + * I, + + + -.1618572254e-18 + + + + + + + -.2570540515e-18 + + + * I}, + + + + + + { + + + QPMS_NORMALISATION_XU + + + , + + + 36 + + + , + + + 36 + + + , + + + -38 + + + , + + + 38 + + + , { + + + 2. + + + , + + + 0.5 + + + , + + + 0.5 + + + }, + + + QPMS_HANKEL_PLUS + + + , + + + .4146334728e+191 + + + + + + + -.4931584782e+191 + + + * I, + + + .0000000000e+00 + + + + + + + .0000000000e+00 + + + * I}, + + + + + + { + + + QPMS_NORMALISATION_XU + + + , + + + -35 + + + , + + + 40 + + + , + + + 35 + + + , + + + 40 + + + , { + + + 2. + + + , + + + 0.5 + + + , + + + 0.5 + + + }, + + + QPMS_HANKEL_PLUS + + + , + + + -.6514262216e-05 + + + + + + + .1374854333e-04 + + + * I, + + + .0000000000e+00 + + + + + + + .0000000000e+00 + + + * I}, + + + + + + { + + + QPMS_NORMALISATION_XU + + + , + + + 32 + + + , + + + 35 + + + , + + + -43 + + + , + + + 45 + + + , { + + + 2. + + + , + + + 0.5 + + + , + + + 0.5 + + + }, + + + QPMS_HANKEL_PLUS + + + , + + + .2762232925e+212 + + + + + + + -.1368895313e+213 + + + * I, + + + .8373862584e+209 + + + + + + + .1689724460e+209 + + + * I}, + + + + + + { + + + QPMS_NORMALISATION_XU + + + , + + + 38 + + + , + + + 42 + + + , + + + -39 + + + , + + + 45 + + + , { + + + 2. + + + , + + + 0.5 + + + , + + + 0.5 + + + }, + + + QPMS_HANKEL_PLUS + + + , + + + -.2298689786e+235 + + + + + + + .2371029493e+235 + + + * I, + + + .1277697908e+232 + + + + + + + .1238711556e+232 + + + * I}, + + + + + + { + + + QPMS_NORMALISATION_XU + + + , + + + -42 + + + , + + + 42 + + + , + + + 45 + + + , + + + 45 + + + , { + + + 2. + + + , + + + 0.5 + + + , + + + 0.5 + + + }, + + + QPMS_HANKEL_PLUS + + + , + + + .3488835702e-28 + + + + + + + -.1826524477e-28 + + + * I, + + + .0000000000e+00 + + + + + + + .0000000000e+00 + + + * I}, + + + + + + { + + + QPMS_NORMALISATION_XU + + + , + + + -43 + + + , + + + 45 + + + , + + + 41 + + + , + + + 42 + + + , { + + + 2. + + + , + + + 0.5 + + + , + + + 0.5 + + + }, + + + QPMS_HANKEL_PLUS + + + , + + + .5178100899e-22 + + + + + + + .1186503822e-21 + + + * I, + + + .3274958627e-25 + + + + + + + -.1429246656e-25 + + + * I}, + + + + + + { + + + QPMS_NORMALISATION_XU + + + , + + + 48 + + + , + + + 50 + + + , + + + -30 + + + , + + + 49 + + + , { + + + 2. + + + , + + + 0.5 + + + , + + + 0.5 + + + }, + + + QPMS_HANKEL_PLUS + + + , + + + .3393827523e+267 + + + + + + + -.1226717423e+268 + + + * I, + + + -.5718637033e+265 + + + + + + + -.1582113973e+265 + + + * I}, + + + + + + { + + + QPMS_NORMALISATION_XU + + + , + + + -72 + + + , + + + 72 + + + , + + + 1 + + + , + + + 3 + + + , { + + + 2. + + + , + + + 0.5 + + + , + + + 0.5 + + + }, + + + QPMS_HANKEL_PLUS + + + , + + + .6946365327e-42 + + + + + + + -.1782022552e-41 + + + * I, + + + .1833377882e-43 + + + + + + + .7146549596e-44 + + + * I}, + + + + + + { + + + QPMS_NORMALISATION_XU + + + , + + + 42 + + + , + + + 52 + + + , + + + 9 + + + , + + + 81 + + + , { + + + 2. + + + , + + + 0.5 + + + , + + + 0.5 + + + }, + + + QPMS_HANKEL_PLUS + + + , + + + .3656934399e+271 + + + + + + + .3705813223e+271 + + + * I, + + + -.4499925012e+269 + + + + + + + -.4440572037e+269 + + + * I}, + + + + + + { + + + QPMS_NORMALISATION_XU + + + , + + + 18 + + + , + + + 100 + + + , + + + -5 + + + , + + + 45 + + + , { + + + 2. + + + , + + + 0.5 + + + , + + + 0.5 + + + }, + + + QPMS_HANKEL_PLUS + + + , + + + .4118769973e+293 + + + + + + + .7460688240e+293 + + + * I, + + + .5914795871e+291 + + + + + + + -.3265339985e+291 + + + * I}, + + + + + // THE END + + + + + + + + + + + + + + + { + + + QPMS_NORMALISATION_UNDEF + + + , + + + 0 + + + , + + + 0 + + + , + + + 0 + + + , + + + 0 + + + , { + + + 0 + + + , + + + 0 + + + , + + + 0 + + + }, + + + QPMS_BESSEL_UNDEF + + + , + + + 0 + + + + + + + 0 + + + * I, + + + 0 + + + + + + + 0 + + + * I} + + + + + + + \ No newline at end of file diff --git a/qpms/translations.c b/qpms/translations.c index 6ecdc16..39a1f36 100644 --- a/qpms/translations.c +++ b/qpms/translations.c @@ -14,6 +14,11 @@ * [1] Yu-Lin Xu, Journal of Computational Physics 127, 285–298 (1996) */ +/* + * GENERAL TODO: use normalised Legendre functions for Kristensson and Taylor conventions directly + * instead of normalising them here (the same applies for csphase). + */ + static const double sqrtpi = 1.7724538509055160272981674833411451827975494561223871; //static const double ln2 = 0.693147180559945309417232121458176568075500134360255254120; @@ -68,10 +73,46 @@ int qpms_sph_bessel_fill(qpms_bessel_t typ, int lmax, double x, complex double * assert(0); } -// [1], eq. (82) -complex double qpms_trans_single_A_Xu(int m, int n, int mu, int nu, sph_t kdlj, +static inline double qpms_trans_normlogfac(qpms_normalisation_t norm, + int m, int n, int mu, int nu) { + //int csphase = qpms_normalisation_t csphase(norm); // probably not needed here + norm = qpms_normalisation_t_normonly(norm); + switch(norm) { + case QPMS_NORMALISATION_KRISTENSSON: + case QPMS_NORMALISATION_TAYLOR: + return -0.5*(lgamma(n+m+1)-lgamma(n-m+1)+lgamma(nu-mu+1)-lgamma(nu+mu+1)); + break; + case QPMS_NORMALISATION_XU: + return 0; + break; + default: + abort(); + } +} + +static inline double qpms_trans_normfac(qpms_normalisation_t norm, + int m, int n, int mu, int nu) { + int csphase = qpms_normalisation_t_csphase(norm); // FIXME USEME TODO + norm = qpms_normalisation_t_normonly(norm); + double normfac = 1.; + switch(norm) { + case QPMS_NORMALISATION_KRISTENSSON: + normfac *= sqrt((nu*(nu+1.))/(n*(n+1.))); + case QPMS_NORMALISATION_TAYLOR: + normfac *= sqrt((2.*n+1)/(2.*nu+1)); + break; + case QPMS_NORMALISATION_XU: + break; + default: + abort(); + } + + return normfac; +} + +complex double qpms_trans_single_A(qpms_normalisation_t norm, + int m, int n, int mu, int nu, sph_t kdlj, bool r_ge_d, qpms_bessel_t J) { - abort(); // FIXME, THIS IS STILL TAYLOR if(r_ge_d) J = QPMS_BESSEL_REGULAR; double costheta = cos(kdlj.theta); @@ -108,59 +149,16 @@ complex double qpms_trans_single_A_Xu(int m, int n, int mu, int nu, sph_t kdlj, complex double presum = exp(exponent); presum *= cexp(I*(mu-m)*kdlj.phi) * min1pow(m) * ipow(nu+n) / (4*n); - complex double prenormratio = ipow(nu-n) * sqrt(((2.*nu+1)/(2.*n+1))* exp( - lgamma(n+m+1)-lgamma(n-m+1)+lgamma(nu-mu+1)-lgamma(nu+mu+1))); - return (presum / prenormratio) * sum; + double normlogfac = qpms_trans_normlogfac(norm,m,n,mu,nu); + double normfac = qpms_trans_normfac(norm,m,n,mu,nu); + + // int csphase = qpms_normalisation_t_csphase(norm); FIXME EITHER TO NORMFAC OR USE HERE + + + presum *= ipow(n-nu) * (normfac * exp(normlogfac)); + return presum * sum; } - - -complex double qpms_trans_single_A_Kristensson(int m, int n, int mu, int nu, sph_t kdlj, - bool r_ge_d, qpms_bessel_t J) { - abort();// FIXME, THIS IS STILL TAYLOR - if(r_ge_d) J = QPMS_BESSEL_REGULAR; - - double costheta = cos(kdlj.theta); - - int qmax = gaunt_q_max(-m,n,mu,nu); // nemá tu být +m? - // N.B. -m !!!!!! - double a1q[qmax+1]; - int err; - gaunt_xu(-m,n,mu,nu,qmax,a1q,&err); - double a1q0 = a1q[0]; - if (err) abort(); - - double leg[gsl_sf_legendre_array_n(n+nu)]; - if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu,costheta,-1,leg)) abort(); - complex double bes[n+nu+1]; - if (qpms_sph_bessel_fill(J, n+nu, kdlj.r, bes)) abort(); - complex double sum = 0; - for(int q = 0; q <= qmax; ++q) { - int p = n+nu-2*q; - int Pp_order = mu-m; - //if(p < abs(Pp_order)) continue; // FIXME raději nastav lépe meze - assert(p >= abs(Pp_order)); - double a1q_n = a1q[q] / a1q0; - double Pp = leg[gsl_sf_legendre_array_index(p, abs(Pp_order))]; - if (Pp_order < 0) Pp *= min1pow(mu-m) * exp(lgamma(1+p+Pp_order)-lgamma(1+p-Pp_order)); - complex double zp = bes[p]; - complex double summandq = (n*(n+1) + nu*(nu+1) - p*(p+1)) * min1pow(q) * a1q_n * zp * Pp; - sum += summandq; - } - - double exponent=(lgamma(2*n+1)-lgamma(n+2)+lgamma(2*nu+3)-lgamma(nu+2) - +lgamma(n+nu+m-mu+1)-lgamma(n-m+1)-lgamma(nu+mu+1) - +lgamma(n+nu+1) - lgamma(2*(n+nu)+1)); - complex double presum = exp(exponent); - presum *= cexp(I*(mu-m)*kdlj.phi) * min1pow(m) * ipow(nu+n) / (4*n); - - complex double prenormratio = ipow(nu-n) * sqrt(((2.*nu+1)/(2.*n+1))* exp( - lgamma(n+m+1)-lgamma(n-m+1)+lgamma(nu-mu+1)-lgamma(nu+mu+1))); - return (presum / prenormratio) * sum; -} - - - complex double qpms_trans_single_A_Taylor(int m, int n, int mu, int nu, sph_t kdlj, bool r_ge_d, qpms_bessel_t J) { if(r_ge_d) J = QPMS_BESSEL_REGULAR; @@ -207,7 +205,6 @@ complex double qpms_trans_single_A_Taylor(int m, int n, int mu, int nu, sph_t kd // [1], eq. (83) complex double qpms_trans_single_B_Xu(int m, int n, int mu, int nu, sph_t kdlj, bool r_ge_d, qpms_bessel_t J) { - abort(); // FIXME, this is still Taylor if(r_ge_d) J = QPMS_BESSEL_REGULAR; double costheta = cos(kdlj.theta); @@ -257,12 +254,71 @@ complex double qpms_trans_single_B_Xu(int m, int n, int mu, int nu, sph_t kdlj, (4*n)*(n+1)*(n+m+1)); // Taylor normalisation v2, proven to be equivalent - complex double prenormratio = ipow(nu-n) * sqrt(((2.*nu+1)/(2.*n+1))* exp( - lgamma(n+m+1)-lgamma(n-m+1)+lgamma(nu-mu+1)-lgamma(nu+mu+1))); + complex double prenormratio = ipow(nu-n); return (presum / prenormratio) * sum; } +complex double qpms_trans_single_B(qpms_normalisation_t norm, + int m, int n, int mu, int nu, sph_t kdlj, + bool r_ge_d, qpms_bessel_t J) { + if(r_ge_d) J = QPMS_BESSEL_REGULAR; + double costheta = cos(kdlj.theta); + + int q2max = gaunt_q_max(-m-1,n+1,mu+1,nu); + int Qmax = gaunt_q_max(-m,n+1,mu,nu); + double a2q[q2max+1], a3q[Qmax+1], a2q0, a3q0; + int err; + if (mu == nu) { + for (int q = 0; q <= q2max; ++q) + a2q[q] = 0; + a2q0 = 1; + } + else { + gaunt_xu(-m-1,n+1,mu+1,nu,q2max,a2q,&err); if (err) abort(); + a2q0 = a2q[0]; + } + gaunt_xu(-m,n+1,mu,nu,Qmax,a3q,&err); if (err) abort(); + a3q0 = a3q[0]; + + double leg[gsl_sf_legendre_array_n(n+nu+1)]; + if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1,costheta,-1,leg)) abort(); + complex double bes[n+nu+2]; + if (qpms_sph_bessel_fill(J, n+nu+1, kdlj.r, bes)) abort(); + + complex double sum = 0; + for (int q = 0; q <= Qmax; ++q) { + int p = n+nu-2*q; + double a2q_n = a2q[q]/a2q0; + double a3q_n = a3q[q]/a3q0; + complex double zp_ = bes[p+1]; + int Pp_order_ = mu-m; + //if(p+1 < abs(Pp_order_)) continue; // FIXME raději nastav lépe meze + assert(p+1 >= abs(Pp_order_)); + double Pp_ = leg[gsl_sf_legendre_array_index(p+1, abs(Pp_order_))]; + if (Pp_order_ < 0) Pp_ *= min1pow(mu-m) * exp(lgamma(1+1+p+Pp_order_)-lgamma(1+1+p-Pp_order_)); + complex double summandq = ((2*(n+1)*(nu-mu)*a2q_n + -(-nu*(nu+1) - n*(n+3) - 2*mu*(n+1)+p*(p+3))* a3q_n) + *min1pow(q) * zp_ * Pp_); + sum += summandq; + } + + double exponent=(lgamma(2*n+3)-lgamma(n+2)+lgamma(2*nu+3)-lgamma(nu+2) + +lgamma(n+nu+m-mu+2)-lgamma(n-m+1)-lgamma(nu+mu+1) + +lgamma(n+nu+2) - lgamma(2*(n+nu)+3)); + complex double presum = exp(exponent); + presum *= cexp(I*(mu-m)*kdlj.phi) * min1pow(m) * ipow(nu+n+1) / ( + (4*n)*(n+1)*(n+m+1)); + + double normlogfac = qpms_trans_normlogfac(norm,m,n,mu,nu); + double normfac = qpms_trans_normfac(norm,m,n,mu,nu); + + // int csphase = qpms_normalisation_t_csphase(norm); FIXME EITHER TO NORMFAC OR USE HERE + + presum *= ipow(n-nu) * (exp(normlogfac) * normfac); + + return presum * sum; +} complex double qpms_trans_single_B_Taylor(int m, int n, int mu, int nu, sph_t kdlj, bool r_ge_d, qpms_bessel_t J) { @@ -354,6 +410,105 @@ static inline size_t qpms_trans_calculator_index_yyu(const qpms_trans_calculator #define SQ(x) ((x)*(x)) + +static void qpms_trans_calculator_multipliers_A_general( + qpms_normalisation_t norm, + complex double *dest, int m, int n, int mu, int nu, int qmax) { + assert(qmax == gaunt_q_max(-m,n,mu,nu)); + double a1q[qmax+1]; + int err; + gaunt_xu(-m,n,mu,nu,qmax,a1q,&err); + if (err) abort(); + double a1q0 = a1q[0]; + + int csphase = qpms_normalisation_t_csphase(norm); //TODO FIXME use this + norm = qpms_normalisation_t_normonly(norm); + double normlogfac = qpms_trans_normlogfac(norm,m,n,mu,nu); + double normfac = qpms_trans_normfac(norm,m,n,mu,nu); + + // TODO use csphase to modify normfac here!!!! + // normfac = xxx ? -normfac : normfac; + normfac *= min1pow(m+n); + double exponent=(lgamma(2*n+1)-lgamma(n+2)+lgamma(2*nu+3)-lgamma(nu+2) + +lgamma(n+nu+m-mu+1)-lgamma(n-m+1)-lgamma(nu+mu+1) + +lgamma(n+nu+1) - lgamma(2*(n+nu)+1)) + + normlogfac; + double presum = exp(exponent); + presum *= normfac / (4.*n); + + for(int q = 0; q <= qmax; q++) { + int p = n+nu-2*q; + int Pp_order = mu - m; + assert(p >= abs(Pp_order)); + double a1q_n = a1q[q] / a1q0; + // Assuming non_normalized legendre polynomials! + double Ppfac = (Pp_order >= 0) ? 1 : + min1pow(mu-m) * exp(lgamma(1+p+Pp_order)-lgamma(1+p-Pp_order)); + double summandfac = (n*(n+1) + nu*(nu+1) - p*(p+1)) * min1pow(q) * a1q_n; + dest[q] = presum * summandfac * Ppfac; + // FIXME I might not need complex here + } +} + + +static void qpms_trans_calculator_multipliers_B_general( + qpms_normalisation_t norm, + complex double *dest, int m, int n, int mu, int nu, int Qmax) { + assert(Qmax == gaunt_q_max(-m,n+1,mu,nu)); + int q2max = gaunt_q_max(-m-1,n+1,mu+1,nu); + assert(Qmax == q2max); + // FIXME remove the q2max variable altogether, as it is probably equal + // to Qmax + double a2q[q2max+1], a3q[Qmax+1], a2q0, a3q0; + int err; + if (mu == nu) { + for (int q = 0; q <= q2max; ++q) + a2q[q] = 0; + a2q0 = 1; + } + else { + gaunt_xu(-m-1,n+1,mu+1,nu,q2max,a2q,&err); if (err) abort(); + a2q0 = a2q[0]; + } + gaunt_xu(-m,n+1,mu,nu,Qmax,a3q,&err); if (err) abort(); + a3q0 = a3q[0]; + + + int csphase = qpms_normalisation_t_csphase(norm); //TODO FIXME use this + norm = qpms_normalisation_t_normonly(norm); + double normlogfac= qpms_trans_normlogfac(norm,m,n,mu,nu); + double normfac = qpms_trans_normfac(norm,m,n,mu,nu); + // TODO use csphase to modify normfac here!!!! + // normfac = xxx ? -normfac : normfac; + normfac *= min1pow(m+n); + + + + double exponent=(lgamma(2*n+3)-lgamma(n+2)+lgamma(2*nu+3)-lgamma(nu+2) + +lgamma(n+nu+m-mu+2)-lgamma(n-m+1)-lgamma(nu+mu+1) + +lgamma(n+nu+2) - lgamma(2*(n+nu)+3)) + +normlogfac; + complex double presum = exp(exponent); + presum *= I * normfac / ( + (4*n)*(n+1)*(n+m+1)); + + for (int q = 0; q <= Qmax; ++q) { + int p = n+nu-2*q; + double a2q_n = a2q[q]/a2q0; + double a3q_n = a3q[q]/a3q0; + int Pp_order_ = mu-m; + //if(p+1 < abs(Pp_order_)) continue; // FIXME raději nastav lépe meze + assert(p+1 >= abs(Pp_order_)); + double Ppfac = (Pp_order_ >= 0) ? 1 : + + min1pow(mu-m) * exp(lgamma(1+1+p+Pp_order_)-lgamma(1+1+p-Pp_order_)); + double summandq = ((2*(n+1)*(nu-mu)*a2q_n + -(-nu*(nu+1) - n*(n+3) - 2*mu*(n+1)+p*(p+3))* a3q_n) + *min1pow(q)); + dest[q] = Ppfac * summandq * presum; + } +} + //#if 0 static void qpms_trans_calculator_multipliers_A_Taylor( complex double *dest, int m, int n, int mu, int nu, int qmax) { @@ -472,9 +627,16 @@ static void qpms_trans_calculator_multipliers_B_Taylor( int qpms_trans_calculator_multipliers_A(qpms_normalisation_t norm, complex double *dest, int m, int n, int mu, int nu, int qmax) { switch (norm) { case QPMS_NORMALISATION_TAYLOR: +#ifdef USE_SEPARATE_TAYLOR qpms_trans_calculator_multipliers_A_Taylor(dest,m,n,mu,nu,qmax); return 0; break; +#endif + case QPMS_NORMALISATION_XU: + case QPMS_NORMALISATION_KRISTENSSON: + qpms_trans_calculator_multipliers_A_general(norm, dest, m, n, mu, nu, qmax); + return 0; + break; default: abort(); } @@ -484,9 +646,16 @@ int qpms_trans_calculator_multipliers_A(qpms_normalisation_t norm, complex doubl int qpms_trans_calculator_multipliers_B(qpms_normalisation_t norm, complex double *dest, int m, int n, int mu, int nu, int qmax) { switch (norm) { case QPMS_NORMALISATION_TAYLOR: +#ifdef USE_SEPARATE_TAYLOR qpms_trans_calculator_multipliers_B_Taylor(dest,m,n,mu,nu,qmax); return 0; break; +#endif + case QPMS_NORMALISATION_XU: + case QPMS_NORMALISATION_KRISTENSSON: + qpms_trans_calculator_multipliers_B_general(norm, dest, m, n, mu, nu, qmax); + return 0; + break; default: abort(); } @@ -586,7 +755,10 @@ complex double qpms_trans_calculator_get_A_buf(const qpms_trans_calculator *c, // TODO warn? return NAN+I*NAN; switch(c->normalisation) { + // TODO use normalised legendre functions for Taylor and Kristensson case QPMS_NORMALISATION_TAYLOR: + case QPMS_NORMALISATION_KRISTENSSON: + case QPMS_NORMALISATION_XU: { double costheta = cos(kdlj.theta); if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu, @@ -632,6 +804,8 @@ complex double qpms_trans_calculator_get_B_buf(const qpms_trans_calculator *c, return NAN+I*NAN; switch(c->normalisation) { case QPMS_NORMALISATION_TAYLOR: + case QPMS_NORMALISATION_KRISTENSSON: + case QPMS_NORMALISATION_XU: { double costheta = cos(kdlj.theta); if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1, @@ -661,6 +835,8 @@ int qpms_trans_calculator_get_AB_buf_p(const qpms_trans_calculator *c, } switch(c->normalisation) { case QPMS_NORMALISATION_TAYLOR: + case QPMS_NORMALISATION_KRISTENSSON: + case QPMS_NORMALISATION_XU: { double costheta = cos(kdlj.theta); if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1, diff --git a/qpms/translations.h b/qpms/translations.h index 4904c61..b63593d 100644 --- a/qpms/translations.h +++ b/qpms/translations.h @@ -7,7 +7,7 @@ #include // TODO replace the xplicit "Taylor" functions with general, -// taking qpms_bessel_t argument. +// taking qpms_normalisation_t argument. complex double qpms_trans_single_A_Taylor(qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, sph_t kdlj, bool r_ge_d, qpms_bessel_t J); @@ -20,6 +20,12 @@ complex double qpms_trans_single_A_Taylor_ext(qpms_m_t m, qpms_l_t n, qpms_m_t m complex double qpms_trans_single_B_Taylor_ext(qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, double kdlj_r, double kdlj_th, double kdlj_phi, int r_ge_d, int J); +complex double qpms_trans_single_A(qpms_normalisation_t norm, qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, sph_t kdlj, + bool r_ge_d, qpms_bessel_t J); + +complex double qpms_trans_single_B(qpms_normalisation_t norm, qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, sph_t kdlj, + bool r_ge_d, qpms_bessel_t J); + typedef struct qpms_trans_calculator { qpms_normalisation_t normalisation; qpms_l_t lMax;