"Basis fields" for finite systems.

There seem to be race conditions in the prange'd cython parts.

qpms/qpms_c.pyx

@@ -1289,8 +1289,9 @@ cdef class _ScatteringSystemAtOmega:
         return results.reshape(evalpos.shape)
 
     @boundscheck(False)
-    def scattered_field_basis(self, evalpos, blochvector=None, btyp=QPMS_HANKEL_PLUS):
+    def scattered_field_basis(self, evalpos, blochvector=None, particle=None, btyp=QPMS_HANKEL_PLUS):
         # TODO examples
+        # FIXME periodic case not implemented
         """Evaluate scattered field "basis" 
 
         This function enables the evaluation of "scattered" fields
@@ -1304,25 +1305,31 @@ cdef class _ScatteringSystemAtOmega:
             Evaluation points in cartesian coordinates.
         blochvector: array_like or None
             Bloch vector, must be supplied (non-None) for periodic systems, else None.
+        particle_index: int or None (default), optional
+            A valid particle index; if specified, only the contribution of the given particle
+            is evaluated. Not applicable for periodic arrays.
         btyp: BesselType, optional
             Kind of the waves. Defaults to BesselType.HANKEL_PLUS.
 
         Returns
         ------- 
-        ndarray of complex, with the shape `evalpos.shape[:-1] + (self.fecv_size, 3)`
+        ndarray of complex, with the shape `evalpos.shape[:-1] + (n, 3)`
             "Basis" fields at the positions given in `evalpos`, in cartesian coordinates.
+            `n` here stays either for `self.fecv_size` (if `particle==None`)
+            or `len(self.bspec_pi(particle))`.
         """
-        if(btyp != QPMS_HANKEL_PLUS):
+        #if(btyp != QPMS_HANKEL_PLUS): # #TODO, IIRC not supported only for periodic systems
-            raise NotImplementedError("Only first kind Bessel function-based fields are supported")
+        #    raise NotImplementedError("Only first kind Bessel function-based fields are supported")
         cdef qpms_bessel_t btyp_c = BesselType(btyp)
-        cdef Py_ssize_t fecv_size = self.fecv_size
+        cdef Py_ssize_t basissize = self.fecv_size if particle is None else len(self.bspec_pi(particle))
+        cdef qpms_ss_pi_t pi = particle
         evalpos = np.array(evalpos, dtype=float, copy=False)
         if evalpos.shape[-1] != 3:
             raise ValueError("Last dimension of evalpos has to be 3")
         cdef np.ndarray[double,ndim=2] evalpos_a = evalpos.reshape(-1,3)
-        cdef np.ndarray[complex, ndim=3] results = np.empty((evalpos_a.shape[0], fecv_size, 3), dtype=complex)
+        cdef np.ndarray[complex, ndim=3] results = np.empty((evalpos_a.shape[0], basissize, 3), dtype=complex)
         cdef ccart3_t *res
-        res = <ccart3_t *> malloc(fecv_size*sizeof(ccart3_t))
+        res = <ccart3_t *> malloc(basissize*sizeof(ccart3_t))
         cdef cart3_t pos
         cdef Py_ssize_t i, j
         with nogil, wraparound(False), parallel():
@@ -1330,14 +1337,24 @@ cdef class _ScatteringSystemAtOmega:
                 pos.x = evalpos_a[i,0]
                 pos.y = evalpos_a[i,1]
                 pos.z = evalpos_a[i,2]
-                qpms_scatsysw_scattered_field_basis(res, self.ssw, btyp_c, pos)
+                if particle is None:
-                for j in range(fecv_size):
+                    qpms_scatsysw_scattered_field_basis(res, self.ssw, btyp_c, pos)
+                else:
+                    qpms_scatsysw_scattered_field_basis_pi(res, self.ssw, pi, btyp_c, pos)
+                for j in range(basissize):
                     results[i,j,0] = res[j].x
                     results[i,j,1] = res[j].y
                     results[i,j,2] = res[j].z
         free(res)
-        return results.reshape(evalpos.shape[:-1] + (self.fecv_size, 3))
+        return results.reshape(evalpos.shape[:-1] + (basissize, 3))
-    
+
+    def bspec_pi(self, pi):
+        return self.ss_pyref.bspec_pi(pi)
+
+    property bspecs:
+        def __get__(self):
+            return self.ss_pyref.bspecs
+
 
 cdef class ScatteringMatrix:
     '''

qpms/qpms_cdefs.pxd

@@ -712,6 +712,10 @@ cdef extern from "scatsystem.h":
             qpms_bessel_t btyp, cdouble wavenumber, cart3_t where) nogil
     qpms_errno_t qpms_scatsysw_scattered_field_basis(ccart3_t *target, const qpms_scatsys_at_omega_t *ssw,
             qpms_bessel_t btyp, cart3_t where) nogil
+    qpms_errno_t qpms_scatsys_scattered_field_basis_pi(ccart3_t *target, const qpms_scatsys_t *ss,
+            qpms_ss_pi_t pi, qpms_bessel_t btyp, cdouble wavenumber, cart3_t where) nogil
+    qpms_errno_t qpms_scatsysw_scattered_field_basis_pi(ccart3_t *target, const qpms_scatsys_at_omega_t *ssw,
+            qpms_ss_pi_t pi, qpms_bessel_t btyp, cart3_t where) nogil
     qpms_errno_t qpms_scatsyswk_scattered_field_basis(ccart3_t *target, const qpms_scatsys_at_omega_k_t *sswk,
             qpms_bessel_t btyp, cart3_t where) nogil
     double qpms_ss_adjusted_eta(const qpms_scatsys_t *ss, cdouble wavenumber, const double *wavevector) nogil

qpms/scatsystem.c

@@ -2049,6 +2049,37 @@ ccart3_t qpms_scatsysw_scattered_E(const qpms_scatsys_at_omega_t *ssw,
       cvf, where);
 }
 
+qpms_errno_t qpms_scatsys_scattered_field_basis_pi(
+    ccart3_t *target,
+    const qpms_scatsys_t *ss,
+    const qpms_ss_pi_t pi,
+    const qpms_bessel_t btyp,
+    const complex double k,
+    const cart3_t where
+    ) {
+  qpms_ss_ensure_nonperiodic_a(ss, "qpms_scatsyswk_scattered_field_basis()");
+  QPMS_UNTESTED;
+  const qpms_vswf_set_spec_t *bspec = qpms_ss_bspec_pi(ss, pi);
+  csphvec_t *vswfs_sph; //Single particle contributions in spherical coordinates
+  QPMS_CRASHING_MALLOC(vswfs_sph, bspec->n * sizeof(*vswfs_sph));
+  const cart3_t particle_pos = ss->p[pi].pos;
+  const csph_t kr = sph_cscale(k, cart2sph(
+      cart3_substract(where, particle_pos)));
+  QPMS_ENSURE_SUCCESS(qpms_uvswf_fill(vswfs_sph, bspec, kr, btyp));
+  for(size_t i = 0; i < bspec->n; ++i)
+      target[i] = csphvec2ccart_csph(vswfs_sph[i], kr);
+  
+  free(vswfs_sph);
+  return QPMS_SUCCESS;
+}
+
+qpms_errno_t qpms_scatsysw_scattered_field_basis_pi(
+    ccart3_t *target, const qpms_scatsys_at_omega_t *ssw, const qpms_ss_pi_t pi,
+    const qpms_bessel_t btyp, const cart3_t where) {
+  return qpms_scatsys_scattered_field_basis_pi(target, ssw->ss, pi, btyp,
+      ssw->wavenumber, where);
+}
+
 qpms_errno_t qpms_scatsys_scattered_field_basis(
     ccart3_t *target,
     const qpms_scatsys_t *ss,
@@ -2058,20 +2089,9 @@ qpms_errno_t qpms_scatsys_scattered_field_basis(
     ) {
   qpms_ss_ensure_nonperiodic_a(ss, "qpms_scatsyswk_scattered_field_basis()");
   QPMS_UNTESTED;
-  csphvec_t *vswfs_sph; //Single particle contributions in spherical coordinates
+  for (qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi) 
-  QPMS_CRASHING_CALLOC(vswfs_sph, ss->max_bspecn, sizeof(*vswfs_sph));
+    QPMS_ENSURE_SUCCESS(qpms_scatsys_scattered_field_basis_pi(
-
+          target + ss->fecv_pstarts[pi], ss, pi, btyp, k, where));
-  for (qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi) {
-    const qpms_vswf_set_spec_t *bspec = qpms_ss_bspec_pi(ss, pi);
-    const cart3_t particle_pos = ss->p[pi].pos;
-    const csph_t kr = sph_cscale(k, cart2sph(
-          cart3_substract(where, particle_pos)));
-    QPMS_ENSURE_SUCCESS(qpms_uvswf_fill(vswfs_sph, bspec, kr, btyp));
-    for(size_t i = 0; i < bspec->n; ++i)
-      target[ss->fecv_pstarts[pi] + i] = csphvec2ccart_csph(vswfs_sph[i], kr);
-  }
-  
-  free(vswfs_sph);
   return QPMS_SUCCESS;
 }
 

qpms/scatsystem.h

@@ -757,12 +757,77 @@ qpms_errno_t qpms_scatsys_scattered_field_basis(
 		cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
 		);
 
+/// Evaluates a "basis" for electric field of a given particle at a given point.
+/**
+ * This function evaluates all the included VSWFs from a particle in the system, evaluated 
+ * at a given point. 
+ *
+ * Finite systems only.
+ *
+ * \see qpms_scatsys_scattered_field_basis() for the whole-array equivalent.
+ *
+ * \see qpms_scatsysw_scattered_field_basis_pi()
+ */
+qpms_errno_t qpms_scatsys_scattered_field_basis_pi(
+		ccart3_t *target, ///< Target array of length at least `qpms_ss_bspec_pi(ssw->ss, pi)->n`
+		const qpms_scatsys_t *ss,
+		qpms_ss_pi_t pi, ///< Particle index
+		qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS).
+		complex double wavenumber, ///< Wavenumber of the background medium
+		cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
+		);
+
 /// Evaluates a "basis" for electric field at a given point.
 /**
  * This function evaluates all the included VSWFs from the particles in the system, evaluated 
  * at a given point. Taking a linear combination of these with the coefficients \a scattcoeff_full[]
  * would be equivalent to the result of qpms_scatsysw_scattered_E().
  *
+ * Note that this might require a relatively large amount of memory. Depending
+ * on the application, qpms_scatsysw_scattered_field_basis_pi() might be a better
+ * alternative.
+ *
+ * \see qpms_scatsysw_scattered_field_basis_pi() for the single-particle equivalent.
+ *
+ * \see qpms_scatsys_scattered_field_basis()
+ *
+ * \see qpms_scatsyswk_scattered_field_basis() for periodic systems.
+ *
+ */
+qpms_errno_t qpms_scatsysw_scattered_field_basis(
+		ccart3_t *target, ///< Target array of length \a ss->fecv_size
+		const qpms_scatsys_at_omega_t *ssw,
+		qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS).
+		cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
+		);
+
+/// Evaluates a "basis" for electric field of a given particle at a given point.
+/**
+ * This function evaluates all the included VSWFs from a particle in the system, evaluated 
+ * at a given point. 
+ *
+ * Finite systems only.
+ *
+ * \see qpms_scatsysw_scattered_field_basis() for the whole-array equivalent.
+ *
+ * \see qpms_scatsys_scattered_field_basis_pi()
+ */
+qpms_errno_t qpms_scatsysw_scattered_field_basis_pi(
+		ccart3_t *target, ///< Target array of length at least `qpms_ss_bspec_pi(ssw->ss, pi)->n`
+		const qpms_scatsys_at_omega_t *ssw,
+		qpms_ss_pi_t pi, ///< Particle index
+		qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS).
+		cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
+		);
+
+/// Evaluates a "basis" for electric field at a given point.
+/**
+ * This function evaluates all the included VSWFs from the particles in the system, evaluated 
+ * at a given point. Taking a linear combination of these with the coefficients \a scattcoeff_full[]
+ * would be equivalent to the result of qpms_scatsysw_scattered_E().
+ *
+ * \see qpms_scatsysw_scattered_field_basis_pi()
+ *
  * \see qpms_scatsys_scattered_field_basis()
  *
  * \see qpms_scatsyswk_scattered_field_basis() for periodic systems.