Allow different Bessel kind in scattered_E

Former-commit-id: afbc59cf542daec6e89d4a48f01acfb7e515819c

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

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":

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, 
+    const csphvec_t E_sph = qpms_eval_uvswf(bspec, particle_cv, kr, btyp);
-        QPMS_HANKEL_PLUS);
     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,
+          ss->c, s, &dipspec, 1, bspec, dipspecn, k, particle_pos, where, btyp));
-          QPMS_HANKEL_PLUS));
 
     for(size_t i = 0; i < bspec->n; ++i)
-      for(size_t j = 0; j < 3; ++j){
+      for(size_t j = 0; j < dipspecn; ++j){
-        ccart3_t summand = ccart3_scale(particle_cv[i] * s[3*i+j], regdipoles_0[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) {
+    qpms_bessel_t btyp, const complex double *cvf,  const cart3_t where) {
-  return qpms_scatsys_scattered_E__alt(ssw->ss, ssw->wavenumber,
+  return qpms_scatsys_scattered_E__alt(ssw->ss, btyp, ssw->wavenumber,
       cvf, where);
 }
 

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.
 		);