diff --git a/qpms/test_vswf_translations.c b/qpms/test_vswf_translations.c new file mode 100644 index 0000000..7f22e25 --- /dev/null +++ b/qpms/test_vswf_translations.c @@ -0,0 +1,258 @@ +#include "translations.h" +#include "vswf.h" +#include +#include +#include +#include +#include +#include "vectors.h" +#include "string.h" +#include "indexing.h" + +char *normstr(qpms_normalisation_t norm) { + //int csphase = qpms_normalisation_t_csphase(norm); + norm = qpms_normalisation_t_normonly(norm); + switch (norm) { + case QPMS_NORMALISATION_NONE: + return "none"; + case QPMS_NORMALISATION_SPHARM: + return "spharm"; + case QPMS_NORMALISATION_POWER: + return "power"; + default: + return "!!!undef!!!"; + } +} + +int test_sphwave_translation(const qpms_trans_calculator *c, qpms_bessel_t wavetype, + cart3_t o2minuso1, int npoints, cart3_t *o1points); +//int test_planewave_decomposition(cart3_t k, ccart3_t E, qpms_l_t lMax, qpms_normalisation_t norm, int npoints, cart3_t *points); +//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); + +int main() { + gsl_rng *rng = gsl_rng_alloc(gsl_rng_ranlxs0); + gsl_rng_set(rng, 666); + + qpms_l_t lMax = 8; + //qpms_l_t viewlMax = 2; + int npoints = 10; + double sigma = 1.3; + + cart3_t o2minuso1; + o2minuso1.x = gsl_ran_gaussian(rng, sigma); + o2minuso1.y = gsl_ran_gaussian(rng, sigma); + o2minuso1.z = gsl_ran_gaussian(rng, sigma); + + cart3_t points[npoints]; + double relerrs[npoints]; + memset(points, 0, npoints * sizeof(cart3_t)); + points[0].x = points[1].y = points[2].z = 1.; + double relerrthreshold = 1e-11; + for (unsigned i = 3; i < npoints; ++i) { + cart3_t *w = points+i; + w->x = gsl_ran_gaussian(rng, sigma); + w->y = gsl_ran_gaussian(rng, sigma); + w->z = gsl_ran_gaussian(rng, sigma); + } + + for(int use_csbit = 0; use_csbit <= 1; ++use_csbit) { + for(int i = 1; i <= 3; ++i){ + qpms_normalisation_t norm = i | (use_csbit ? QPMS_NORMALISATION_T_CSBIT : 0); + qpms_trans_calculator *c = qpms_trans_calculator_init(lMax, norm); + for(int J = 1; J <= 4; ++J) + test_sphwave_translation(c, J, o2minuso1, npoints, points); + qpms_trans_calculator_free(c); + } + + } + gsl_rng_free(rng); +} + +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)); + 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... +} + + + + + +#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/translations.c b/qpms/translations.c index 2328b1e..442ce8e 100644 --- a/qpms/translations.c +++ b/qpms/translations.c @@ -627,7 +627,7 @@ 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) { + switch (qpms_normalisation_t_normonly(norm)) { case QPMS_NORMALISATION_TAYLOR: #ifdef USE_SEPARATE_TAYLOR qpms_trans_calculator_multipliers_A_Taylor(dest,m,n,mu,nu,qmax); @@ -646,7 +646,7 @@ 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) { + switch (qpms_normalisation_t_normonly(norm)) { case QPMS_NORMALISATION_TAYLOR: #ifdef USE_SEPARATE_TAYLOR qpms_trans_calculator_multipliers_B_Taylor(dest,m,n,mu,nu,qmax); @@ -837,7 +837,7 @@ int qpms_trans_calculator_get_AB_buf_p(const qpms_trans_calculator *c, // TODO warn? different return value? return 0; } - switch(c->normalisation) { + switch(qpms_normalisation_t_normonly(c->normalisation)) { case QPMS_NORMALISATION_TAYLOR: case QPMS_NORMALISATION_KRISTENSSON: case QPMS_NORMALISATION_XU: diff --git a/qpms/vecprint.c b/qpms/vecprint.c index d7f87cb..8b33506 100644 --- a/qpms/vecprint.c +++ b/qpms/vecprint.c @@ -5,7 +5,7 @@ void print_csphvec(csphvec_t v) { - printf("(%f+%fj)r̂ + (%f+%fj)θ̂ + (%f+%fj)φ̂", + printf("(%g+%gj)r̂ + (%g+%gj)θ̂ + (%g+%gj)φ̂", creal(v.rc), cimag(v.rc), creal(v.thetac), cimag(v.thetac), creal(v.phic), cimag(v.phic) @@ -15,12 +15,12 @@ void print_csphvec(csphvec_t v) void print_cart3(cart3_t v) { - printf("%fx̂ + %fŷ + %fẑ", v.x, v.y, v.z); + printf("%gx̂ + %gŷ + %gẑ", v.x, v.y, v.z); } void print_ccart3(ccart3_t v) { - printf("(%f+%fj)x̂ + (%f+%fj)ŷ + (%f+%fj)ẑ", + printf("(%g+%gj)x̂ + (%g+%gj)ŷ + (%g+%gj)ẑ", creal(v.x), cimag(v.x), creal(v.y), cimag(v.y), creal(v.z), cimag(v.z)