Plane wave decomposition now works for the transversal part (not yet longitudinal)

Former-commit-id: 83dadf35d24a66c7576395f6f156f5e7586551cd
This commit is contained in:
Marek Nečada 2018-02-04 11:43:42 +00:00
parent 6198a992f9
commit dbc118e712
4 changed files with 129 additions and 19 deletions

View File

@ -22,30 +22,64 @@ char *normstr(qpms_normalisation_t norm) {
} }
} }
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(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() { int main() {
gsl_rng *rng = gsl_rng_alloc(gsl_rng_ranlxs0); gsl_rng *rng = gsl_rng_alloc(gsl_rng_ranlxs0);
gsl_rng_set(rng, 666); gsl_rng_set(rng, 666);
qpms_l_t lMax = 20; qpms_l_t lMax = 50;
int npoints = 10; int npoints = 10;
double sigma = 1.3; double sigma = 1.3;
cart3_t points[npoints]; cart3_t points[npoints];
double relerrs[npoints];
memset(points, 0, npoints * sizeof(cart3_t)); memset(points, 0, npoints * sizeof(cart3_t));
points[0].x = points[1].y = points[2].z = 1.; points[1].x = points[2].y = points[3].z = 1.;
for (unsigned i = 3; i < npoints; ++i) { double relerrthreshold = 1e-11;
for (unsigned i = 4; i < npoints; ++i) {
cart3_t *w = points+i; cart3_t *w = points+i;
w->x = gsl_ran_gaussian(rng, sigma); w->x = gsl_ran_gaussian(rng, sigma);
w->y = gsl_ran_gaussian(rng, sigma); w->y = gsl_ran_gaussian(rng, sigma);
w->z = gsl_ran_gaussian(rng, sigma); w->z = gsl_ran_gaussian(rng, sigma);
} }
double ksigma = 15;
for(int i = 1; i < 2; ++i) {
cart3_t k = {0, 0, 2};
ccart3_t E = {0., 1.1+I, 0.};
if(i){
k.x = gsl_ran_gaussian(rng, ksigma);
k.y = gsl_ran_gaussian(rng, ksigma);
k.z = gsl_ran_gaussian(rng, ksigma);
E.x = gsl_ran_gaussian(rng, ksigma);
E.x += I* gsl_ran_gaussian(rng, ksigma);
E.y = gsl_ran_gaussian(rng, ksigma);
E.y += I* gsl_ran_gaussian(rng, ksigma);
E.z = gsl_ran_gaussian(rng, ksigma);
E.z += I* gsl_ran_gaussian(rng, ksigma);
}
#if 0
printf("Test wave k = %gx̂ + %gŷ + %gẑ", k.x, k.y, k.z);
printf(", E_0 = (%g+%gj)x̂ + (%g+%gj)ŷ + (%g+%gj)ẑ\n",
creal(E.x),cimag(E.x),
creal(E.y),cimag(E.y),
creal(E.z),cimag(E.z));
printf("%d points, sigma = %g, rel. err. threshold = %g\n", npoints, sigma, relerrthreshold);
for(int i = 1; i <= 3; ++i){
int res = test_planewave_decomposition_silent(k, E, lMax, i, npoints, points, relerrthreshold, relerrs);
printf("\n%s %d/%d %s %s\n", res ? "!!" : "OK", npoints-res, npoints, normstr(i), "");
for(int p = 0; p < npoints; ++p) printf("%.3g ", relerrs[p]);
res = test_planewave_decomposition_silent(k, E, lMax, i |QPMS_NORMALISATION_T_CSBIT, npoints, points, relerrthreshold, relerrs);
printf("\n%s %d/%d %s %s\n", res ? "!!" : "OK", npoints-res, npoints, normstr(i), "with C.-S. phase");
for(int p = 0; p < npoints; ++p) printf("%.3g ", relerrs[p]);
}
#endif
for(int i = 1; i <= 3; ++i){
test_planewave_decomposition(k, E, lMax, i, npoints, points);
test_planewave_decomposition(k, E, lMax, i | QPMS_NORMALISATION_T_CSBIT, npoints, points);
}
cart3_t k = {0, 0, 1};
ccart3_t E = {0., 1., 0.};
for(int i = 1; i <= 3; ++i){
test_planewave_decomposition(k, E, lMax, i, npoints, points);
test_planewave_decomposition(k, E, lMax, i | QPMS_NORMALISATION_T_CSBIT, npoints, points);
} }
gsl_rng_free(rng); gsl_rng_free(rng);
} }
@ -84,9 +118,9 @@ int test_planewave_decomposition(cart3_t k, ccart3_t E, qpms_l_t lMax, qpms_norm
sph_t w_sph = cart2sph(w); sph_t w_sph = cart2sph(w);
printf("Point %d: x = %f, y = %f, z = %f,\n", i, w.x, w.y, w.z); 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); printf(" |r| = %f, θ = %f, φ = %f:\n", w_sph.r, w_sph.theta, w_sph.phi);
double ph_2pi = cart3_dot(k,w); double phase = cart3_dot(k,w);
printf(" k.r = %f\n", ph_2pi); printf(" k.r = %f\n", phase);
double phfac = cexp(2 * M_PI * ph_2pi * I); complex double phfac = cexp(phase * I);
ccart3_t Ew = ccart3_scale(phfac, E); ccart3_t Ew = ccart3_scale(phfac, E);
printf(" pw E(r) = (%f+%fj)x̂ + (%f+%fj)ŷ + (%f+%fj)ẑ", printf(" pw E(r) = (%f+%fj)x̂ + (%f+%fj)ŷ + (%f+%fj)ẑ",
creal(Ew.x),cimag(Ew.x), creal(Ew.x),cimag(Ew.x),
@ -97,6 +131,7 @@ int test_planewave_decomposition(cart3_t k, ccart3_t E, qpms_l_t lMax, qpms_norm
creal(Ew_s.rc), cimag(Ew_s.rc), creal(Ew_s.rc), cimag(Ew_s.rc),
creal(Ew_s.thetac), cimag(Ew_s.thetac), creal(Ew_s.thetac), cimag(Ew_s.thetac),
creal(Ew_s.phic), cimag(Ew_s.phic)); 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, csphvec_t Ew_s_recomp = qpms_eval_vswf(w_sph,
lc, mc, ec, lMax, QPMS_BESSEL_REGULAR, norm); lc, mc, ec, lMax, QPMS_BESSEL_REGULAR, norm);
ccart3_t Ew_recomp = csphvec2ccart(Ew_s_recomp, w_sph); ccart3_t Ew_recomp = csphvec2ccart(Ew_s_recomp, w_sph);
@ -108,13 +143,50 @@ int test_planewave_decomposition(cart3_t k, ccart3_t E, qpms_l_t lMax, qpms_norm
creal(Ew_s_recomp.rc), cimag(Ew_s_recomp.rc), creal(Ew_s_recomp.rc), cimag(Ew_s_recomp.rc),
creal(Ew_s_recomp.thetac), cimag(Ew_s_recomp.thetac), creal(Ew_s_recomp.thetac), cimag(Ew_s_recomp.thetac),
creal(Ew_s_recomp.phic), cimag(Ew_s_recomp.phic)); creal(Ew_s_recomp.phic), cimag(Ew_s_recomp.phic));
printf(" rel. err. magnitude: %e @ r̂, %e @ θ̂, %e @ φ̂\n", double relerrfac = 2/(cabs(Ew_s_recomp.rc) + cabs(Ew_s.rc)
2. * cabs(Ew_s_recomp.rc - Ew_s.rc) +cabs(Ew_s_recomp.thetac) + cabs(Ew_s.thetac)
/(cabs(Ew_s_recomp.rc) + cabs(Ew_s.rc)), +cabs(Ew_s_recomp.phic) + cabs(Ew_s.phic));
2. * cabs(Ew_s_recomp.thetac - Ew_s.thetac) printf(" rel. err. magnitude: %g @ r̂, %g @ θ̂, %g @ φ̂\n",
/(cabs(Ew_s_recomp.thetac) + cabs(Ew_s.thetac)), cabs(Ew_s_recomp.rc - Ew_s.rc) * relerrfac,
2. * cabs(Ew_s_recomp.phic - Ew_s.phic) cabs(Ew_s_recomp.thetac - Ew_s.thetac) * relerrfac,
/(cabs(Ew_s_recomp.phic) + cabs(Ew_s.phic))); cabs(Ew_s_recomp.phic - Ew_s.phic) * relerrfac
);
} }
return 0; 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;
}

33
qpms/vecprint.c Normal file
View File

@ -0,0 +1,33 @@
#include <assert.h>
#include "vectors.h"
#include "complex.h"
#include <stdio.h>
void print_csphvec(csphvec_t v)
{
printf("(%f+%fj)r̂ + (%f+%fj)θ̂ + (%f+%fj)φ̂",
creal(v.rc), cimag(v.rc),
creal(v.thetac), cimag(v.thetac),
creal(v.phic), cimag(v.phic)
);
}
void print_cart3(cart3_t v)
{
printf("%fx̂ + %fŷ + %fẑ", v.x, v.y, v.z);
}
void print_ccart3(ccart3_t v)
{
printf("(%f+%fj)x̂ + (%f+%fj)ŷ + (%f+%fj)ẑ",
creal(v.x), cimag(v.x),
creal(v.y), cimag(v.y),
creal(v.z), cimag(v.z)
);
}
void print_sph(sph_t r)
{
printf("(r=%g, θ=%g, φ=%g)", r.r, r.theta, r.phi);
}

View File

@ -107,4 +107,9 @@ static inline csphvec_t ccart2csphvec(const ccart3_t cartvec, const sph_t at) {
return res; return res;
} }
void print_csphvec(csphvec_t);
void print_ccart3(ccart3_t);
void print_cart3(cart3_t);
void print_sph(sph_t);
#endif //VECTORS_H #endif //VECTORS_H

View File

@ -374,8 +374,8 @@ qpms_errno_t qpms_vecspharm_dual_fill(csphvec_t *const a1target, csphvec_t *cons
p3->phic = 0; p3->phic = 0;
++p3; ++p3;
} }
++pleg; ++ppi; ++ptau;
} }
++pleg; ++ppi; ++ptau;
} }
qpms_pitau_free(pt); qpms_pitau_free(pt);
return QPMS_SUCCESS; return QPMS_SUCCESS;