diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx index 68ba71e..3a60070 100644 --- a/qpms/qpms_c.pyx +++ b/qpms/qpms_c.pyx @@ -1066,6 +1066,7 @@ cdef class _ScatteringSystemAtOmegaK: results[i,2] = res.z return results.reshape(evalpos.shape) + @boundscheck(False) def scattered_field_basis(self, evalpos, btyp=QPMS_HANKEL_PLUS): # TODO examples """Evaluate scattered field "basis" (periodic system) @@ -1237,7 +1238,6 @@ cdef class _ScatteringSystemAtOmega: @boundscheck(False) def scattered_E(self, scatcoeffvector_full, evalpos, blochvector=None, btyp=QPMS_HANKEL_PLUS, bint alt=False): - # FIXME TODO this obviously does not work for periodic systems """Evaluate electric field for a given excitation coefficient vector Parameters @@ -1288,6 +1288,56 @@ cdef class _ScatteringSystemAtOmega: results[i,2] = res.z return results.reshape(evalpos.shape) + @boundscheck(False) + def scattered_field_basis(self, evalpos, blochvector=None, btyp=QPMS_HANKEL_PLUS): + # TODO examples + """Evaluate scattered field "basis" + + This function enables the evaluation of "scattered" fields + generated by the system for many different excitation + coefficients vectors, without the expensive re-evaluation of the respective + translation operators for each excitation coefficient vector. + + Parameters + ---------- + evalpos: array_like of floats and shape (..., 3) + Evaluation points in cartesian coordinates. + blochvector: array_like or None + Bloch vector, must be supplied (non-None) for periodic systems, else None. + 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)` + "Basis" fields at the positions given in `evalpos`, in cartesian coordinates. + """ + if(btyp != QPMS_HANKEL_PLUS): + 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 + 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 ccart3_t *res + res = malloc(fecv_size*sizeof(ccart3_t)) + cdef cart3_t pos + cdef Py_ssize_t i, j + with nogil, wraparound(False), parallel(): + for i in prange(evalpos_a.shape[0]): + 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) + for j in range(fecv_size): + 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)) + cdef class ScatteringMatrix: ''' diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index 1ca0476..7c3cd79 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -708,6 +708,10 @@ cdef extern from "scatsystem.h": const cdouble *f_excitation_vector_full, cart3_t where) nogil ccart3_t qpms_scatsyswk_scattered_E(const qpms_scatsys_at_omega_k_t *sswk, qpms_bessel_t btyp, const cdouble *f_excitation_vector_full, cart3_t where) nogil + qpms_errno_t qpms_scatsys_scattered_field_basis(ccart3_t *target, const qpms_scatsys_t *ss, + 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_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 diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index 81f1ab5..a3d845b 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -2049,6 +2049,38 @@ ccart3_t qpms_scatsysw_scattered_E(const qpms_scatsys_at_omega_t *ssw, cvf, where); } +qpms_errno_t qpms_scatsys_scattered_field_basis( + ccart3_t *target, + const qpms_scatsys_t *ss, + 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; + csphvec_t *vswfs_sph; //Single particle contributions in spherical coordinates + QPMS_CRASHING_CALLOC(vswfs_sph, ss->max_bspecn, sizeof(*vswfs_sph)); + + 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_errno_t qpms_scatsysw_scattered_field_basis( + ccart3_t *target, const qpms_scatsys_at_omega_t *ssw, + const qpms_bessel_t btyp, const cart3_t where) { + return qpms_scatsys_scattered_field_basis(target, ssw->ss, btyp, + ssw->wavenumber, where); +} #define DIPSPECN 3 // We have three basis vectors // Evaluates the regular electric dipole waves in the origin. The returned diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index bcb3649..31596fa 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -706,7 +706,7 @@ complex double *qpms_scatsys_incident_field_vector_irrep_packed( * * \return Complex electric field at the point defined by \a evalpoint. * - * \see qpms_scatsysw_scattered_E() + * \see qpms_scatsysw_scattered_E(), qpms_scatsys_scattered_field_basis() * * \see qpms_scatsyswk_scattered_E() for periodic systems. * @@ -727,7 +727,7 @@ ccart3_t qpms_scatsys_scattered_E( * * \return Complex electric field at the point defined by \a evalpoint. * - * \see qpms_scatsys_scattered_E() + * \see qpms_scatsys_scattered_E(), qpms_scatsysw_scattered_field_basis() * * \see qpms_scatsyswk_scattered_E() for periodic systems. */ @@ -738,6 +738,43 @@ ccart3_t qpms_scatsysw_scattered_E( 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_scatsys_scattered_E(). + * + * \see qpms_scatsysw_scattered_field_basis() + * + * \see qpms_scatsyswk_scattered_field_basis() for periodic systems. + * + */ +qpms_errno_t qpms_scatsys_scattered_field_basis( + ccart3_t *target, ///< Target array of length \a ss->fecv_size + 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. + 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_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 scattered electric field at a point, given a full vector of scattered field coefficients. /** * This is an alternative implementation of qpms_scatsys_scattered_E(), and should give the same results