QPMS
/
qpms
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Code reuse

This commit is contained in:
Marek Nečada 2020-07-23 10:00:39 +03:00
parent 868e603f1c
commit 5246229b2d
1 changed files with 50 additions and 78 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

124
qpms/scatsystem.c
View File

@ -2043,6 +2043,38 @@ ccart3_t qpms_scatsysw_scattered_E(const qpms_scatsys_at_omega_t *ssw,
cvf, where); 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 // Alternative implementation, using translation operator and regular dipole waves at zero
ccart3_t qpms_scatsys_scattered_E__alt(const qpms_scatsys_t *ss, ccart3_t qpms_scatsys_scattered_E__alt(const qpms_scatsys_t *ss,
qpms_bessel_t btyp, 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 = {0,0,0};
ccart3_t res_kc = {0,0,0}; // kahan sum compensation ccart3_t res_kc = {0,0,0}; // kahan sum compensation
static const int dipspecn = 3; // We have three basis vectors ccart3_t regdipoles_0[DIPSPECN];
// bspec containing only electric dipoles const qpms_vswf_set_spec_t dipspec = qpms_fill_regdipoles_0(regdipoles_0, ss->c->normalisation);
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);
}
complex double *s; // Translation matrix complex double *s; // Translation matrix
QPMS_CRASHING_MALLOC(s, ss->max_bspecn * sizeof(*s) * dipspec.n); 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}; const cart3_t origin_cart = {0, 0, 0};
QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_lc3p( 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 i = 0; i < bspec->n; ++i)
for(size_t j = 0; j < dipspecn; ++j){ for(size_t j = 0; j < dipspec.n; ++j){
ccart3_t summand = ccart3_scale(particle_cv[i] * s[dipspecn*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.x), &(res_kc.x), summand.x);
ckahanadd(&(res.y), &(res_kc.y), summand.y); ckahanadd(&(res.y), &(res_kc.y), summand.y);
ckahanadd(&(res.z), &(res_kc.z), summand.z); ckahanadd(&(res.z), &(res_kc.z), summand.z);
@ -2125,28 +2137,8 @@ ccart3_t qpms_scatsyswk_scattered_E(const qpms_scatsys_at_omega_k_t *sswk,
ccart3_t res = {0,0,0}; ccart3_t res = {0,0,0};
ccart3_t res_kc = {0,0,0}; // kahan sum compensation ccart3_t res_kc = {0,0,0}; // kahan sum compensation
static const int dipspecn = 3; // We have three basis vectors ccart3_t regdipoles_0[DIPSPECN];
// bspec containing only electric dipoles const qpms_vswf_set_spec_t dipspec = qpms_fill_regdipoles_0(regdipoles_0, ss->c->normalisation);
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);
}
complex double *s; // Translation matrix complex double *s; // Translation matrix
QPMS_CRASHING_MALLOC(s, ss->max_bspecn * sizeof(*s) * dipspec.n); 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( QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_e32(
ss->c, s, NULL, ss->c, s, NULL,
&dipspec, 1, bspec, dipspecn, &dipspec, 1, bspec, dipspec.n,
sswk->eta, sswk->ssw->wavenumber, sswk->eta, sswk->ssw->wavenumber,
cart3xy2cart2(ss->per.lattice_basis[0]), cart3xy2cart2(ss->per.lattice_basis[1]), cart3xy2cart2(ss->per.lattice_basis[0]), cart3xy2cart2(ss->per.lattice_basis[1]),
cart2_from_double_array(sswk->k), cart3_substract(where, particle_pos) /*CHECKSIGN*/, cart2_from_double_array(sswk->k), cart3_substract(where, particle_pos) /*CHECKSIGN*/,
maxR, maxK)); maxR, maxK));
for(size_t i = 0; i < bspec->n; ++i) 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){
ccart3_t summand = ccart3_scale(particle_cv[i] * s[dipspecn*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.x), &(res_kc.x), summand.x);
ckahanadd(&(res.y), &(res_kc.y), summand.y); ckahanadd(&(res.y), &(res_kc.y), summand.y);
ckahanadd(&(res.z), &(res_kc.z), summand.z); ckahanadd(&(res.z), &(res_kc.z), summand.z);
@ -2201,28 +2193,8 @@ qpms_errno_t qpms_scatsyswk_scattered_field_basis(
//ccart3_t res = {0,0,0}; //ccart3_t res = {0,0,0};
//ccart3_t res_kc = {0,0,0}; // kahan sum compensation //ccart3_t res_kc = {0,0,0}; // kahan sum compensation
static const int dipspecn = 3; // We have three basis vectors ccart3_t regdipoles_0[DIPSPECN];
// bspec containing only electric dipoles const qpms_vswf_set_spec_t dipspec = qpms_fill_regdipoles_0(regdipoles_0, ss->c->normalisation);
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);
}
complex double *s; // Translation matrix complex double *s; // Translation matrix
QPMS_CRASHING_MALLOC(s, ss->max_bspecn * sizeof(*s) * dipspec.n); 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( QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_e32(
ss->c, s, NULL, ss->c, s, NULL,
&dipspec, 1, bspec, dipspecn, &dipspec, 1, bspec, dipspec.n,
sswk->eta, sswk->ssw->wavenumber, sswk->eta, sswk->ssw->wavenumber,
cart3xy2cart2(ss->per.lattice_basis[0]), cart3xy2cart2(ss->per.lattice_basis[1]), cart3xy2cart2(ss->per.lattice_basis[0]), cart3xy2cart2(ss->per.lattice_basis[1]),
cart2_from_double_array(sswk->k), cart3_substract(where, particle_pos) /*CHECKSIGN*/, cart2_from_double_array(sswk->k), cart3_substract(where, particle_pos) /*CHECKSIGN*/,
maxR, maxK)); maxR, maxK));
for(size_t i = 0; i < bspec->n; ++i) 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], target[ss->fecv_pstarts[pi] + i] = ccart3_add(target[ss->fecv_pstarts[pi] + i],
ccart3_scale(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[dipspecn*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.x), &(res_kc.x), summand.x);
//ckahanadd(&(res.y), &(res_kc.y), summand.y); //ckahanadd(&(res.y), &(res_kc.y), summand.y);
//ckahanadd(&(res.z), &(res_kc.z), summand.z); //ckahanadd(&(res.z), &(res_kc.z), summand.z);