diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx index 4e6da2c..a9c97d2 100644 --- a/qpms/qpms_c.pyx +++ b/qpms/qpms_c.pyx @@ -11,7 +11,7 @@ from .qpms_cdefs cimport * from .cyquaternions cimport IRot3, CQuat from .cybspec cimport BaseSpec from .cycommon cimport make_c_string -from .cycommon import string_c2py, PointGroupClass +from .cycommon import string_c2py, PointGroupClass, BesselType from .cytmatrices cimport CTMatrix, TMatrixFunction, TMatrixGenerator, TMatrixInterpolator from .cymaterials cimport EpsMuGenerator, EpsMu from libc.stdlib cimport malloc, free, calloc @@ -710,7 +710,8 @@ cdef class ScatteringSystem: return retdict - def scattered_E(self, cdouble wavenumber, scatcoeffvector_full, evalpos, bint alt=False): + def scattered_E(self, cdouble wavenumber, scatcoeffvector_full, evalpos, bint alt=False, btyp=BesselType.HANKEL_PLUS): + cdef qpms_bessel_t btyp_c = BesselType(btyp) evalpos = np.array(evalpos, dtype=float, copy=False) if evalpos.shape[-1] != 3: raise ValueError("Last dimension of evalpos has to be 3") @@ -726,9 +727,9 @@ cdef class ScatteringSystem: pos.y = evalpos_a[i,1] pos.z = evalpos_a[i,2] if alt: - res = qpms_scatsys_scattered_E__alt(self.s, wavenumber, &scv_view[0], pos) + res = qpms_scatsys_scattered_E__alt(self.s, btyp_c, wavenumber, &scv_view[0], pos) else: - res = qpms_scatsys_scattered_E(self.s, wavenumber, &scv_view[0], pos) + res = qpms_scatsys_scattered_E(self.s, btyp_c, wavenumber, &scv_view[0], pos) results[i,0] = res.x results[i,1] = res.y results[i,2] = res.z @@ -867,7 +868,8 @@ cdef class _ScatteringSystemAtOmega: def translation_matrix_full(self, blochvector = None): return self.ss_pyref.translation_matrix_full(wavenumber=self.wavenumber, blochvector=blochvector) - def scattered_E(self, scatcoeffvector_full, evalpos, bint alt=False): + def scattered_E(self, scatcoeffvector_full, evalpos, bint alt=False, btyp=QPMS_HANKEL_PLUS): + cdef qpms_bessel_t btyp_c = BesselType(btyp) evalpos = np.array(evalpos, dtype=float, copy=False) if evalpos.shape[-1] != 3: raise ValueError("Last dimension of evalpos has to be 3") @@ -883,9 +885,9 @@ cdef class _ScatteringSystemAtOmega: pos.y = evalpos_a[i,1] pos.z = evalpos_a[i,2] if alt: - res = qpms_scatsysw_scattered_E__alt(self.ssw, &scv_view[0], pos) + res = qpms_scatsysw_scattered_E__alt(self.ssw, btyp_c, &scv_view[0], pos) else: - res = qpms_scatsysw_scattered_E(self.ssw, &scv_view[0], pos) + res = qpms_scatsysw_scattered_E(self.ssw, btyp_c, &scv_view[0], pos) results[i,0] = res.x results[i,1] = res.y results[i,2] = res.z diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index c6d806a..a08f41d 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -625,13 +625,13 @@ cdef extern from "scatsystem.h": cdouble omega_centre, double omega_rr, double omega_ri, size_t contour_npoints, double rank_tol, size_t rank_sel_min, double res_tol) const qpms_vswf_set_spec_t *qpms_ss_bspec_pi(const qpms_scatsys_t *ss, qpms_ss_pi_t pi) - ccart3_t qpms_scatsys_scattered_E(const qpms_scatsys_t *ss, cdouble wavenumber, + ccart3_t qpms_scatsys_scattered_E(const qpms_scatsys_t *ss, qpms_bessel_t btyp, cdouble wavenumber, const cdouble *f_excitation_vector_full, cart3_t where) - ccart3_t qpms_scatsysw_scattered_E(const qpms_scatsys_at_omega_t *ssw, + ccart3_t qpms_scatsysw_scattered_E(const qpms_scatsys_at_omega_t *ssw, qpms_bessel_t btyp, const cdouble *f_excitation_vector_full, cart3_t where) - ccart3_t qpms_scatsys_scattered_E__alt(const qpms_scatsys_t *ss, cdouble wavenumber, + ccart3_t qpms_scatsys_scattered_E__alt(const qpms_scatsys_t *ss, qpms_bessel_t btyp, cdouble wavenumber, const cdouble *f_excitation_vector_full, cart3_t where) - ccart3_t qpms_scatsysw_scattered_E__alt(const qpms_scatsys_at_omega_t *ssw, + ccart3_t qpms_scatsysw_scattered_E__alt(const qpms_scatsys_at_omega_t *ssw, qpms_bessel_t btyp, const cdouble *f_excitation_vector_full, cart3_t where) cdef extern from "ewald.h": diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index 8f0844d..17a5e47 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -1987,6 +1987,7 @@ complex double *qpms_scatsysw_apply_Tmatrices_full( } ccart3_t qpms_scatsys_scattered_E(const qpms_scatsys_t *ss, + qpms_bessel_t btyp, const complex double k, const complex double *cvf, const cart3_t where @@ -2002,8 +2003,7 @@ ccart3_t qpms_scatsys_scattered_E(const qpms_scatsys_t *ss, const csph_t kr = sph_cscale(k, cart2sph( cart3_substract(where, particle_pos))); - const csphvec_t E_sph = qpms_eval_uvswf(bspec, particle_cv, kr, - QPMS_HANKEL_PLUS); + const csphvec_t E_sph = qpms_eval_uvswf(bspec, particle_cv, kr, btyp); const ccart3_t E_cart = csphvec2ccart_csph(E_sph, kr); ckahanadd(&(res.x), &(res_kc.x), E_cart.x); ckahanadd(&(res.y), &(res_kc.y), E_cart.y); @@ -2013,13 +2013,15 @@ ccart3_t qpms_scatsys_scattered_E(const qpms_scatsys_t *ss, } ccart3_t qpms_scatsysw_scattered_E(const qpms_scatsys_at_omega_t *ssw, + qpms_bessel_t btyp, const complex double *cvf, const cart3_t where) { - return qpms_scatsys_scattered_E(ssw->ss, ssw->wavenumber, + return qpms_scatsys_scattered_E(ssw->ss, btyp, ssw->wavenumber, cvf, where); } // 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, const complex double k, const complex double *cvf, const cart3_t where @@ -2028,25 +2030,26 @@ 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 = 3, + .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, + .lMax=1, .lMax_M=1, .lMax_N=1, .lMax_L=-1, .capacity=0, .norm = ss->c->normalisation, }; - ccart3_t regdipoles_0[3]; { + 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[3]; + 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 < 3; ++i) + for(int i = 0; i < dipspecn; ++i) regdipoles_0[i] = csphvec2ccart(regdipoles_0_sph[i], origin_sph); } @@ -2061,12 +2064,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, 3, k, particle_pos, where, - QPMS_HANKEL_PLUS)); + ss->c, s, &dipspec, 1, bspec, dipspecn, k, particle_pos, where, btyp)); for(size_t i = 0; i < bspec->n; ++i) - for(size_t j = 0; j < 3; ++j){ - ccart3_t summand = ccart3_scale(particle_cv[i] * s[3*i+j], regdipoles_0[j]); + for(size_t j = 0; j < dipspecn; ++j){ + ccart3_t summand = ccart3_scale(particle_cv[i] * s[dipspecn*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); @@ -2077,8 +2079,8 @@ ccart3_t qpms_scatsys_scattered_E__alt(const qpms_scatsys_t *ss, } ccart3_t qpms_scatsysw_scattered_E__alt(const qpms_scatsys_at_omega_t *ssw, - const complex double *cvf, const cart3_t where) { - return qpms_scatsys_scattered_E__alt(ssw->ss, ssw->wavenumber, + qpms_bessel_t btyp, const complex double *cvf, const cart3_t where) { + return qpms_scatsys_scattered_E__alt(ssw->ss, btyp, ssw->wavenumber, cvf, where); } diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index d7260a0..b40e347 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -705,6 +705,7 @@ complex double *qpms_scatsys_incident_field_vector_irrep_packed( */ ccart3_t qpms_scatsys_scattered_E( const qpms_scatsys_t *ss, + qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS). complex double wavenumber, ///< Wavenumber of the background medium. const complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$. cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated. @@ -723,6 +724,7 @@ ccart3_t qpms_scatsys_scattered_E( */ ccart3_t qpms_scatsysw_scattered_E( const qpms_scatsys_at_omega_t *ssw, + qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS). const complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$. cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated. ); @@ -738,6 +740,7 @@ ccart3_t qpms_scatsysw_scattered_E( */ ccart3_t qpms_scatsys_scattered_E__alt( const qpms_scatsys_t *ss, + qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS). complex double wavenumber, ///< Wavenumber of the background medium. const complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$. cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated. @@ -754,6 +757,7 @@ ccart3_t qpms_scatsys_scattered_E__alt( */ ccart3_t qpms_scatsysw_scattered_E__alt( const qpms_scatsys_at_omega_t *ssw, + qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS). const complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$. cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated. );