diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index f36c0de..3fd27fb 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -2043,6 +2043,38 @@ ccart3_t qpms_scatsysw_scattered_E(const qpms_scatsys_at_omega_t *ssw, cvf, where); } + +static const int DIPSPECN = 3; +// Evaluates the regular electric dipole waves in the origin. The returned +// value is not to be freed as in the usual case. +static inline const qpms_vswf_set_spec_t qpms_fill_regdipoles_0( + ccart3_t regdipoles_0[DIPSPECN], qpms_normalisation_t normalisation) { + static const int dipspecn = DIPSPECN; // We have three basis vectors + // bspec containing only electric dipoles + const qpms_vswf_set_spec_t dipspec = { + .n = dipspecn, + .ilist = (qpms_uvswfi_t[]){ + qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, -1, 1), + qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, 0, 1), + qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, +1, 1), + }, + .lMax=1, .lMax_M=0, .lMax_N=1, .lMax_L=-1, + .capacity=0, + .norm = normalisation, + }; + + const sph_t origin_sph = {.r = 0, .theta = M_PI_2, .phi=0}; // Should work with any theta/phi (TESTWORTHY) + csphvec_t regdipoles_0_sph[dipspecn]; + QPMS_ENSURE_SUCCESS(qpms_uvswf_fill(regdipoles_0_sph, &dipspec, + sph2csph(origin_sph), QPMS_BESSEL_REGULAR)); + for(int i = 0; i < dipspecn; ++i) + regdipoles_0[i] = csphvec2ccart(regdipoles_0_sph[i], origin_sph); + + return dipspec; +} + + + // Alternative implementation, using translation operator and regular dipole waves at zero ccart3_t qpms_scatsys_scattered_E__alt(const qpms_scatsys_t *ss, qpms_bessel_t btyp, @@ -2055,28 +2087,8 @@ ccart3_t qpms_scatsys_scattered_E__alt(const qpms_scatsys_t *ss, ccart3_t res = {0,0,0}; ccart3_t res_kc = {0,0,0}; // kahan sum compensation - static const int dipspecn = 3; // We have three basis vectors - // bspec containing only electric dipoles - const qpms_vswf_set_spec_t dipspec = { - .n = dipspecn, - .ilist = (qpms_uvswfi_t[]){ - qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, -1, 1), - qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, 0, 1), - qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, +1, 1), - }, - .lMax=1, .lMax_M=0, .lMax_N=1, .lMax_L=-1, - .capacity=0, - .norm = ss->c->normalisation, - }; - - ccart3_t regdipoles_0[dipspecn]; { - const sph_t origin_sph = {.r = 0, .theta = M_PI_2, .phi=0}; // Should work with any theta/phi (TESTWORTHY) - csphvec_t regdipoles_0_sph[dipspecn]; - QPMS_ENSURE_SUCCESS(qpms_uvswf_fill(regdipoles_0_sph, &dipspec, - sph2csph(origin_sph), QPMS_BESSEL_REGULAR)); - for(int i = 0; i < dipspecn; ++i) - regdipoles_0[i] = csphvec2ccart(regdipoles_0_sph[i], origin_sph); - } + ccart3_t regdipoles_0[DIPSPECN]; + const qpms_vswf_set_spec_t dipspec = qpms_fill_regdipoles_0(regdipoles_0, ss->c->normalisation); complex double *s; // Translation matrix QPMS_CRASHING_MALLOC(s, ss->max_bspecn * sizeof(*s) * dipspec.n); @@ -2089,11 +2101,11 @@ ccart3_t qpms_scatsys_scattered_E__alt(const qpms_scatsys_t *ss, const cart3_t origin_cart = {0, 0, 0}; QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_lc3p( - ss->c, s, &dipspec, 1, bspec, dipspecn, k, particle_pos, where, btyp)); + ss->c, s, &dipspec, 1, bspec, dipspec.n, k, particle_pos, where, btyp)); for(size_t i = 0; i < bspec->n; ++i) - for(size_t j = 0; j < dipspecn; ++j){ - ccart3_t summand = ccart3_scale(particle_cv[i] * s[dipspecn*i+j], regdipoles_0[j]); + for(size_t j = 0; j < dipspec.n; ++j){ + ccart3_t summand = ccart3_scale(particle_cv[i] * s[dipspec.n*i+j], regdipoles_0[j]); ckahanadd(&(res.x), &(res_kc.x), summand.x); ckahanadd(&(res.y), &(res_kc.y), summand.y); ckahanadd(&(res.z), &(res_kc.z), summand.z); @@ -2124,29 +2136,9 @@ ccart3_t qpms_scatsyswk_scattered_E(const qpms_scatsys_at_omega_k_t *sswk, QPMS_NOT_IMPLEMENTED("Only 2D-periodic lattices implemented"); ccart3_t res = {0,0,0}; ccart3_t res_kc = {0,0,0}; // kahan sum compensation - - static const int dipspecn = 3; // We have three basis vectors - // bspec containing only electric dipoles - const qpms_vswf_set_spec_t dipspec = { - .n = dipspecn, - .ilist = (qpms_uvswfi_t[]){ - qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, -1, 1), - qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, 0, 1), - qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, +1, 1), - }, - .lMax=1, .lMax_M=0, .lMax_N=1, .lMax_L=-1, - .capacity=0, - .norm = ss->c->normalisation, - }; - - ccart3_t regdipoles_0[dipspecn]; { - const sph_t origin_sph = {.r = 0, .theta = M_PI_2, .phi=0}; // Should work with any theta/phi (TESTWORTHY) - csphvec_t regdipoles_0_sph[dipspecn]; - QPMS_ENSURE_SUCCESS(qpms_uvswf_fill(regdipoles_0_sph, &dipspec, - sph2csph(origin_sph), QPMS_BESSEL_REGULAR)); - for(int i = 0; i < dipspecn; ++i) - regdipoles_0[i] = csphvec2ccart(regdipoles_0_sph[i], origin_sph); - } + + ccart3_t regdipoles_0[DIPSPECN]; + const qpms_vswf_set_spec_t dipspec = qpms_fill_regdipoles_0(regdipoles_0, ss->c->normalisation); complex double *s; // Translation matrix QPMS_CRASHING_MALLOC(s, ss->max_bspecn * sizeof(*s) * dipspec.n); @@ -2168,15 +2160,15 @@ ccart3_t qpms_scatsyswk_scattered_E(const qpms_scatsys_at_omega_k_t *sswk, QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_e32( ss->c, s, NULL, - &dipspec, 1, bspec, dipspecn, + &dipspec, 1, bspec, dipspec.n, sswk->eta, sswk->ssw->wavenumber, cart3xy2cart2(ss->per.lattice_basis[0]), cart3xy2cart2(ss->per.lattice_basis[1]), cart2_from_double_array(sswk->k), cart3_substract(where, particle_pos) /*CHECKSIGN*/, maxR, maxK)); for(size_t i = 0; i < bspec->n; ++i) - for(size_t j = 0; j < dipspecn; ++j){ - ccart3_t summand = ccart3_scale(particle_cv[i] * s[dipspecn*i+j], regdipoles_0[j]); + for(size_t j = 0; j < dipspec.n; ++j){ + ccart3_t summand = ccart3_scale(particle_cv[i] * s[dipspec.n*i+j], regdipoles_0[j]); ckahanadd(&(res.x), &(res_kc.x), summand.x); ckahanadd(&(res.y), &(res_kc.y), summand.y); ckahanadd(&(res.z), &(res_kc.z), summand.z); @@ -2200,29 +2192,9 @@ qpms_errno_t qpms_scatsyswk_scattered_field_basis( QPMS_NOT_IMPLEMENTED("Only 2D-periodic lattices implemented"); //ccart3_t res = {0,0,0}; //ccart3_t res_kc = {0,0,0}; // kahan sum compensation - - static const int dipspecn = 3; // We have three basis vectors - // bspec containing only electric dipoles - const qpms_vswf_set_spec_t dipspec = { - .n = dipspecn, - .ilist = (qpms_uvswfi_t[]){ - qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, -1, 1), - qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, 0, 1), - qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, +1, 1), - }, - .lMax=1, .lMax_M=0, .lMax_N=1, .lMax_L=-1, - .capacity=0, - .norm = ss->c->normalisation, - }; - - ccart3_t regdipoles_0[dipspecn]; { - const sph_t origin_sph = {.r = 0, .theta = M_PI_2, .phi=0}; // Should work with any theta/phi (TESTWORTHY) - csphvec_t regdipoles_0_sph[dipspecn]; - QPMS_ENSURE_SUCCESS(qpms_uvswf_fill(regdipoles_0_sph, &dipspec, - sph2csph(origin_sph), QPMS_BESSEL_REGULAR)); - for(int i = 0; i < dipspecn; ++i) - regdipoles_0[i] = csphvec2ccart(regdipoles_0_sph[i], origin_sph); - } + + ccart3_t regdipoles_0[DIPSPECN]; + const qpms_vswf_set_spec_t dipspec = qpms_fill_regdipoles_0(regdipoles_0, ss->c->normalisation); complex double *s; // Translation matrix QPMS_CRASHING_MALLOC(s, ss->max_bspecn * sizeof(*s) * dipspec.n); @@ -2246,17 +2218,17 @@ qpms_errno_t qpms_scatsyswk_scattered_field_basis( QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_e32( ss->c, s, NULL, - &dipspec, 1, bspec, dipspecn, + &dipspec, 1, bspec, dipspec.n, sswk->eta, sswk->ssw->wavenumber, cart3xy2cart2(ss->per.lattice_basis[0]), cart3xy2cart2(ss->per.lattice_basis[1]), cart2_from_double_array(sswk->k), cart3_substract(where, particle_pos) /*CHECKSIGN*/, maxR, maxK)); for(size_t i = 0; i < bspec->n; ++i) - for(size_t j = 0; j < dipspecn; ++j){ + for(size_t j = 0; j < dipspec.n; ++j){ target[ss->fecv_pstarts[pi] + i] = ccart3_add(target[ss->fecv_pstarts[pi] + i], - ccart3_scale(s[dipspecn*i+j], regdipoles_0[j])); - //ccart3_t summand = ccart3_scale(particle_cv[i] * s[dipspecn*i+j], regdipoles_0[j]); + ccart3_scale(s[dipspec.n*i+j], regdipoles_0[j])); + //ccart3_t summand = ccart3_scale(particle_cv[i] * s[dipspec.n*i+j], regdipoles_0[j]); //ckahanadd(&(res.x), &(res_kc.x), summand.x); //ckahanadd(&(res.y), &(res_kc.y), summand.y); //ckahanadd(&(res.z), &(res_kc.z), summand.z);