From f943cc0cdb70539dc9326471f8850eed6ee027b1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Sun, 6 May 2018 17:34:35 +0000 Subject: [PATCH] Found a working set of parameters for the VSWF translation. Former-commit-id: 3b36597e152bd89e6251e34d549fe1db3ecfd955 --- BUGS.rst | 26 +++ .../test_planewave_decomposition.c | 1 + qpms/{ => tests}/test_vswf_translations.c | 1 + qpms/tests/vswf_translation_test_rays.c | 196 +----------------- ...vswf_translation_test_rays_combocompile.sh | 43 ++++ qpms/translations.c | 62 +++++- 6 files changed, 137 insertions(+), 192 deletions(-) rename qpms/{ => tests}/test_planewave_decomposition.c (98%) rename qpms/{ => tests}/test_vswf_translations.c (98%) create mode 100644 qpms/tests/vswf_translation_test_rays_combocompile.sh diff --git a/BUGS.rst b/BUGS.rst index 93eb87a..2d84484 100644 --- a/BUGS.rst +++ b/BUGS.rst @@ -8,6 +8,32 @@ B(1,0,n,n)/B(1,0,n,-n) == -(2n)! at (x,y,z) = (x,0,0) (expected plus or minus 1). A-coefficients seem to behave correctly. +Translation coefficients inconsistent +------------------------------------- +The translation coefficients currently do not work correctly except for certain +combinations of additional i-factors (defined by the [AB][NMF][123] macros +in translations.c) and only for certain normalisations. +QPMS_NORMALISATION_KRISTENSSON_CS does not work at all +QPMS_NORMALISATION_NONE_CS does not work at all +QPMS_NORMALISATION_TAYLOR_CS works for the following macros defined: + AN1 AM0 BN1 BM0 BF0 + AN1 AM0 BN3 BM2 BF2 + AN3 AM2 BN1 BM0 BF0 + AN3 AM2 BN3 BM2 BF2 +QPMS_NORMALISATION_TAYLOR works for the following macros defined: + AN1 AM2 BN1 BM3 BF0 + AN1 AM2 BN3 BM0 BF2 + AN3 AM0 BN1 BM2 BF0 + AN3 AM0 BN3 BM0 BF2 + +The default behaviour is now that the QPMS_NORMALISATION_TAYLOR_CS works. + +Longitudinal waves +------------------ +Plane wave decompositions gives wrong value on the longitudinal part. +The implementation of the L coefficients OR the longitudinal waves +is thus probably wrong. + Scattering result asymmetry --------------------------- The lattice scattering code (such as finitesqlatzsym-scattery.py) produces diff --git a/qpms/test_planewave_decomposition.c b/qpms/tests/test_planewave_decomposition.c similarity index 98% rename from qpms/test_planewave_decomposition.c rename to qpms/tests/test_planewave_decomposition.c index 3e2bd7e..f4a8b92 100644 --- a/qpms/test_planewave_decomposition.c +++ b/qpms/tests/test_planewave_decomposition.c @@ -1,3 +1,4 @@ +// c99 -o test_planewave_decomposition -ggdb -I .. test_planewave_decomposition.c ../vswf.c -lgsl -lm -lblas ../legendre.c ../bessel.c #include #include #include diff --git a/qpms/test_vswf_translations.c b/qpms/tests/test_vswf_translations.c similarity index 98% rename from qpms/test_vswf_translations.c rename to qpms/tests/test_vswf_translations.c index 29f578d..49ea7e7 100644 --- a/qpms/test_vswf_translations.c +++ b/qpms/tests/test_vswf_translations.c @@ -1,3 +1,4 @@ +//c99 -o test_vswf_translations -ggdb -I .. test_vswf_translations.c ../translations.c ../gaunt.c -lgsl -lm -lblas ../vecprint.c ../vswf.c ../legendre.c #include "translations.h" #include "vswf.h" #include diff --git a/qpms/tests/vswf_translation_test_rays.c b/qpms/tests/vswf_translation_test_rays.c index 8351a0e..47ca019 100644 --- a/qpms/tests/vswf_translation_test_rays.c +++ b/qpms/tests/vswf_translation_test_rays.c @@ -1,3 +1,5 @@ +// c99 -o ../../tests/raytests/n_nu0 -ggdb -I .. vswf_translation_test_rays.c ../translations.c ../vswf.c ../gaunt.c ../legendre.c -lgsl -lm -lblas + #include "translations.h" #include "vswf.h" #include @@ -58,7 +60,7 @@ int main(int argc, char **argv) { qpms_y_t y1 = qpms_mn2y(m1, l1); //qpms_l_t viewlMax = 2; - int npoints = 10; + int npoints = 3; double sigma = 4.; //double shiftsigma = 0.01; @@ -66,9 +68,9 @@ int main(int argc, char **argv) { // o2minuso1.x = gsl_ran_gaussian(rng, shiftsigma); // o2minuso1.y = gsl_ran_gaussian(rng, shiftsigma); // o2minuso1.z = gsl_ran_gaussian(rng, shiftsigma); - o2minuso1_max.x = 6; - o2minuso1_max.y = 0; - o2minuso1_max.z = 0; + o2minuso1_max.x = 1; + o2minuso1_max.y = 2; + o2minuso1_max.z = 5; int nraysteps = 512; cart3_t o2mo1_step = cart3_scale(1./nraysteps, o2minuso1_max); @@ -235,189 +237,3 @@ int main(int argc, char **argv) { return 0; } -#if 0 -int test_sphwave_translation(const qpms_trans_calculator *c, qpms_bessel_t wavetype, - cart3_t sc, int npoints, cart3_t *points) { - puts("=============================================================="); - printf("Test translation o2-o1 = %fx̂ + %fŷ + %fẑ", sc.x, sc.y, sc.z); - sph_t ss = cart2sph(sc); - printf(" = %fr̂ @ o1, θ = %f, φ = %f\n", ss.r, ss.theta, ss.phi); - printf("lMax = %d, norm: %s, csphase = %d\n", - (int)c->lMax, normstr(c->normalisation), qpms_normalisation_t_csphase(c->normalisation)); - printf("wave type J = %d\n", wavetype); - - qpms_l_t lMax = c->lMax; - qpms_y_t nelem = c->nelem; - csphvec_t N1[nelem], /* N2[nelem], */ M1[nelem] /*, M2[nelem]*/; - - for (int i = 0; i < npoints; i++) { - printf("-------- Point %d --------\n", i); - cart3_t w1c = points[i]; - // cart3_t w2c = cart3_add(w1c, cart3_scale(-1, sc)); - cart3_t w2c = cart3_add(w1c, sc); - sph_t w1s = cart2sph(w1c); - sph_t w2s = cart2sph(w2c); - printf(" = %fx̂ + %fŷ + %fẑ @o1\n", w1c.x, w1c.y, w1c.z); - printf(" = %fx̂ + %fŷ + %fẑ @o2\n", w2c.x, w2c.y, w2c.z); - printf(" = %fr̂ @ o1, θ = %f, φ = %f\n", w1s.r, w1s.theta, w1s.phi); - printf(" = %fr̂ @ o2, θ = %f, φ = %f\n", w2s.r, w2s.theta, w2s.phi); - printf("Outside the sphere centered in o1 intersecting o2: %s; by %f\n", (w1s.r > ss.r) ? "true" : "false", - w1s.r - ss.r); - if(QPMS_SUCCESS != qpms_vswf_fill(NULL, M1, N1, lMax, w1s, wavetype, c->normalisation)) - abort(); // original wave set - - for(qpms_y_t y1 = 0; y1 < nelem; ++y1) { //index of the wave originating in o1 that will be reconstructed in o2 - qpms_m_t m1; - qpms_l_t l1; - qpms_y2mn_p(y1, &m1, &l1); - printf("*** wave l = %d, m = %d ***\n", l1, m1); - - complex double A[nelem], B[nelem]; - for(qpms_y_t y2 = 0; y2 < nelem; ++y2){ - qpms_m_t m2; qpms_l_t l2; - qpms_y2mn_p(y2, &m2, &l2); - if(qpms_trans_calculator_get_AB_p(c, A+y2, B+y2, m2, l2, m1, l1, ss, (w1s.r > ss.r), wavetype)) - abort(); - } - - printf("M = "); - print_csphvec(M1[y1]); - printf(" @ o1\n = "); - ccart3_t M1c = csphvec2ccart(M1[y1], w1s); - print_ccart3(M1c); - printf("\n = "); - csphvec_t M1s2 = ccart2csphvec(M1c, w2s); - print_csphvec(M1s2); - printf(" @ o2\n"); - csphvec_t M2s2 = qpms_eval_vswf(w2s, NULL, A, B, lMax, wavetype, c->normalisation); - printf("Mr= "); - print_csphvec(M2s2); - printf(" @ o2\n"); - - printf("N = "); - print_csphvec(N1[y1]); - printf(" @ o1\n = "); - ccart3_t N1c = csphvec2ccart(N1[y1], w1s); - print_ccart3(N1c); - printf("\n = "); - csphvec_t N1s2 = ccart2csphvec(N1c, w2s); - print_csphvec(N1s2); - printf(" @o2\nNr= "); - csphvec_t N2s2 = qpms_eval_vswf(w2s, NULL, B, A, lMax, wavetype, c->normalisation); - print_csphvec(N2s2); - printf(" @o2\n"); - } - } - - return 0; // FIXME something more meaningful here... -} -#endif - -#if 0 -int test_planewave_decomposition(cart3_t k, ccart3_t E, qpms_l_t lMax, qpms_normalisation_t norm, int npoints, cart3_t *points){ - qpms_y_t nelem = qpms_lMax2nelem(lMax); - complex double lc[nelem], mc[nelem], ec[nelem]; - if (QPMS_SUCCESS != - qpms_planewave2vswf_fill_cart(k, E, lc, mc, ec, lMax, norm)) { - printf("Error\n"); - return -1; - } - printf("==============================================================\n"); - printf("Test wave k = %fx̂ + %fŷ + %fẑ", k.x, k.y, k.z); - sph_t k_sph = cart2sph(k); - printf(" = %fr̂ @ θ = %f, φ = %f\n", k_sph.r, k_sph.theta, k_sph.phi); - printf(" E_0 = (%f+%fj)x̂ + (%f+%fj)ŷ + (%f+%fj)ẑ", - creal(E.x),cimag(E.x), - creal(E.y),cimag(E.y), - creal(E.z),cimag(E.z)); - csphvec_t E_s = ccart2csphvec(E, k_sph); - printf(" = (%f+%fj)r̂ + (%f+%fj)θ̂ + (%f+%fj)φ̂ @ k\n", - creal(E_s.rc), cimag(E_s.rc), creal(E_s.thetac), cimag(E_s.thetac), - creal(E_s.phic), cimag(E_s.phic)); - printf(" lMax = %d, norm: %s, csphase = %d\n", - (int)lMax, normstr(norm), qpms_normalisation_t_csphase(norm)); - printf("a_L: "); - for(qpms_y_t y = 0; y < nelem; ++y) printf("%g+%gj ", creal(lc[y]), cimag(lc[y])); - printf("\na_M: "); - for(qpms_y_t y = 0; y < nelem; ++y) printf("%g+%gj ", creal(mc[y]), cimag(mc[y])); - printf("\na_N: "); - for(qpms_y_t y = 0; y < nelem; ++y) printf("%g+%gj ", creal(ec[y]), cimag(ec[y])); - printf("\n"); - for (int i = 0; i < npoints; i++) { - cart3_t w = points[i]; - sph_t w_sph = cart2sph(w); - printf("Point %d: x = %f, y = %f, z = %f,\n", i, w.x, w.y, w.z); - printf(" |r| = %f, θ = %f, φ = %f:\n", w_sph.r, w_sph.theta, w_sph.phi); - double phase = cart3_dot(k,w); - printf(" k.r = %f\n", phase); - complex double phfac = cexp(phase * I); - ccart3_t Ew = ccart3_scale(phfac, E); - printf(" pw E(r) = (%f+%fj)x̂ + (%f+%fj)ŷ + (%f+%fj)ẑ", - creal(Ew.x),cimag(Ew.x), - creal(Ew.y),cimag(Ew.y), - creal(Ew.z),cimag(Ew.z)); - csphvec_t Ew_s = ccart2csphvec(Ew, w_sph); - printf(" = (%f+%fj)r̂ + (%f+%fj)θ̂ + (%f+%fj)φ̂ @ r\n", - creal(Ew_s.rc), cimag(Ew_s.rc), - creal(Ew_s.thetac), cimag(Ew_s.thetac), - creal(Ew_s.phic), cimag(Ew_s.phic)); - w_sph.r *= k_sph.r; /// NEVER FORGET THIS!!! - csphvec_t Ew_s_recomp = qpms_eval_vswf(w_sph, - lc, mc, ec, lMax, QPMS_BESSEL_REGULAR, norm); - ccart3_t Ew_recomp = csphvec2ccart(Ew_s_recomp, w_sph); - printf(" rec E(r) = (%f+%fj)x̂ + (%f+%fj)ŷ + (%f+%fj)ẑ", - creal(Ew_recomp.x),cimag(Ew_recomp.x), - creal(Ew_recomp.y),cimag(Ew_recomp.y), - creal(Ew_recomp.z),cimag(Ew_recomp.z)); - printf(" = (%f+%fj)r̂ + (%f+%fj)θ̂ + (%f+%fj)φ̂ @ r\n", - creal(Ew_s_recomp.rc), cimag(Ew_s_recomp.rc), - creal(Ew_s_recomp.thetac), cimag(Ew_s_recomp.thetac), - creal(Ew_s_recomp.phic), cimag(Ew_s_recomp.phic)); - double relerrfac = 2/(cabs(Ew_s_recomp.rc) + cabs(Ew_s.rc) - +cabs(Ew_s_recomp.thetac) + cabs(Ew_s.thetac) - +cabs(Ew_s_recomp.phic) + cabs(Ew_s.phic)); - printf(" rel. err. magnitude: %g @ r̂, %g @ θ̂, %g @ φ̂\n", - cabs(Ew_s_recomp.rc - Ew_s.rc) * relerrfac, - cabs(Ew_s_recomp.thetac - Ew_s.thetac) * relerrfac, - cabs(Ew_s_recomp.phic - Ew_s.phic) * relerrfac - ); - } - return 0; -} - -int test_planewave_decomposition_silent(cart3_t k, ccart3_t E, qpms_l_t lMax, qpms_normalisation_t norm, int npoints, cart3_t *points, double relerrthreshold, double *relerrs) { - qpms_y_t nelem = qpms_lMax2nelem(lMax); - int failcount = 0; - complex double lc[nelem], mc[nelem], ec[nelem]; - if (QPMS_SUCCESS != - qpms_planewave2vswf_fill_cart(k, E, lc, mc, ec, lMax, norm)) { - printf("Error\n"); - return -1; - } - sph_t k_sph = cart2sph(k); - csphvec_t E_s = ccart2csphvec(E, k_sph); - for (int i = 0; i < npoints; i++) { - cart3_t w = points[i]; - sph_t w_sph = cart2sph(w); - w_sph.r *= k_sph.r; - double phase = cart3_dot(k,w); - complex double phfac = cexp(phase * I); - ccart3_t Ew = ccart3_scale(phfac, E); - csphvec_t Ew_s = ccart2csphvec(Ew, w_sph); - csphvec_t Ew_s_recomp = qpms_eval_vswf(w_sph, - lc, mc, ec, lMax, QPMS_BESSEL_REGULAR, norm); - ccart3_t Ew_recomp = csphvec2ccart(Ew_s_recomp, w_sph); - double relerrfac = 2/(cabs(Ew_s_recomp.rc) + cabs(Ew_s.rc) - +cabs(Ew_s_recomp.thetac) + cabs(Ew_s.thetac) - +cabs(Ew_s_recomp.phic) + cabs(Ew_s.phic)); - - double relerr = (cabs(Ew_s_recomp.rc - Ew_s.rc) - + cabs(Ew_s_recomp.thetac - Ew_s.thetac) - + cabs(Ew_s_recomp.phic - Ew_s.phic) - ) * relerrfac; - if(relerrs) relerrs[i] = relerr; - if(relerr > relerrthreshold) ++failcount; - } - return failcount; -} -#endif diff --git a/qpms/tests/vswf_translation_test_rays_combocompile.sh b/qpms/tests/vswf_translation_test_rays_combocompile.sh new file mode 100644 index 0000000..514c5b4 --- /dev/null +++ b/qpms/tests/vswf_translation_test_rays_combocompile.sh @@ -0,0 +1,43 @@ +for an in 0 1 2 3 ; do + for am in 0 2 ; do + for bn in 0 1 2 3 ; do + for bm in 0 2 ; do + for bf in 0 1 2 3 ; do + c99 -o ../../tests/raytests/z/AN${an}M${am}BN${bn}M${bm}F${bf} -DAN${an} -DAM${am} -DBN${bn} -DBM${bm} -DBF${bf} -ggdb -O2 -I .. vswf_translation_test_rays.c ../translations.c ../vswf.c ../gaunt.c ../legendre.c -lgsl -lm -lblas + done + done + done + done +done + +export NORM="129" +for an in 0 1 2 3 ; do + for am in 0 2 ; do + for bn in 0 1 2 3 ; do + for bm in 0 2 ; do + for bf in 0 1 2 3 ; do + ./AN${an}M${am}BN${bn}M${bm}F${bf} an${an}m${am}bn${bn}m${bm}f${bf}_${NORM} 1 0 ${NORM} + sed -e 's/\./,/g' an${an}m${am}bn${bn}m${bm}f${bf}_${NORM}.0 > an${an}m${am}bn${bn}m${bm}f${bf}_${NORM}.0.tsv + sed -e 's/\./,/g' an${an}m${am}bn${bn}m${bm}f${bf}_${NORM}.1 > an${an}m${am}bn${bn}m${bm}f${bf}_${NORM}.1.tsv + sed -e 's/\./,/g' an${an}m${am}bn${bn}m${bm}f${bf}_${NORM}.2 > an${an}m${am}bn${bn}m${bm}f${bf}_${NORM}.2.tsv +# for i in an${an}m${am}bn${bn}m${bm}_130.? ; do +# sed -e 's/\./,/g' $i >${i}.tsv +# done + done + done + done + done +done + + +export NORM="129" +export smer="125" +export an=1 am=0 bn=1 bm=0 bf=0 + c99 -o ../../tests/raytests/$smer/AN${an}M${am}BN${bn}M${bm}F${bf} -DAN${an} -DAM${am} -DBN${bn} -DBM${bm} -DBF${bf} -ggdb -I .. vswf_translation_test_rays.c ../translations.c ../vswf.c ../gaunt.c ../legendre.c -lgsl -lm -lblas + cd ../../tests/raytests/$smer + ./AN${an}M${am}BN${bn}M${bm}F${bf} an${an}m${am}bn${bn}m${bm}f${bf}_${NORM} 1 0 ${NORM} + sed -e 's/\./,/g' an${an}m${am}bn${bn}m${bm}f${bf}_${NORM}.0 > an${an}m${am}bn${bn}m${bm}f${bf}_${NORM}.0.tsv + sed -e 's/\./,/g' an${an}m${am}bn${bn}m${bm}f${bf}_${NORM}.1 > an${an}m${am}bn${bn}m${bm}f${bf}_${NORM}.1.tsv + sed -e 's/\./,/g' an${an}m${am}bn${bn}m${bm}f${bf}_${NORM}.2 > an${an}m${am}bn${bn}m${bm}f${bf}_${NORM}.2.tsv + cd - + \ No newline at end of file diff --git a/qpms/translations.c b/qpms/translations.c index 90b6826..c2dd817 100644 --- a/qpms/translations.c +++ b/qpms/translations.c @@ -11,6 +11,35 @@ #include //abort() #include + +/* + * Define macros with additional factors that "should not be there" according + * to the "original" formulae but are needed to work with my vswfs. + * (actually, I don't know whether the error is in using "wrong" implementation + * of vswfs, "wrong" implementation of Xu's translation coefficient formulae, + * error/inconsintency in Xu's paper or something else) + * Anyway, the zeroes give the correct _numerical_ values according to Xu's + * paper tables (without Xu's typos, of course), while + * the predefined macros give the correct translations of the VSWFs. + */ +#if !(defined AN0 || defined AN1 || defined AN2 || defined AN3) +#pragma message "using AN1 macro as default" +#define AN1 +#endif +//#if !(defined AM0 || defined AM2) +//#define AM1 +//#endif +#if !(defined BN0 || defined BN1 || defined BN2 || defined BN3) +#pragma message "using BN1 macro as default" +#define BN1 +#endif +//#if !(defined BM0 || defined BM2) +//#define BM1 +//#endif +//#if !(defined BF0 || defined BF1 || defined BF2 || defined BF3) +//#define BF1 +//#endif + // if defined, the pointer B_multipliers[y] corresponds to the q = 1 element; // otherwise, it corresponds to the q = 0 element, which should be identically zero #ifdef QPMS_PACKED_B_MULTIPLIERS @@ -483,7 +512,18 @@ static void qpms_trans_calculator_multipliers_A_general( 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 /* * ipow(n-nu) */; + dest[q] = presum * summandfac * Ppfac +#ifdef AN1 + * ipow(n-nu) +#elif defined AN2 + * min1pow(-n+nu) +#elif defined AN3 + * ipow (nu - n) +#endif +#ifdef AM2 + * min1pow(-m+mu) +#endif //NNU + ; // FIXME I might not need complex here } } @@ -533,7 +573,25 @@ void qpms_trans_calculator_multipliers_B_general( * (isq(n+nu+1)-isq(p+1)) ); dest[q-BQ_OFFSET] = presum * t * Ppfac_ - * cruzan_bfactor(-m,n,mu,nu,p) * ipow(p+1) /* * ipow(n-nu) */; + * cruzan_bfactor(-m,n,mu,nu,p) * ipow(p+1) +#ifdef BN1 + * ipow(n-nu) +#elif defined BN2 + * min1pow(-n+nu) +#elif defined BN3 + * ipow (nu - n) +#endif +#ifdef BM2 + * min1pow(-m+mu) +#endif +#ifdef BF1 + * I +#elif defined BF2 + * (-1) +#elif defined BF3 + * (-I) +#endif + ;// NNU } }