From 35ffc61fafbd5667880556eeaccaf0c15e1911b3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Sat, 21 Mar 2020 08:11:16 +0200 Subject: [PATCH] Scatsystem scattered E field methods Former-commit-id: 6521916a81a7c8022bfaefb7102ae657eef99516 --- qpms/qpms_c.pyx | 38 ++++++++++++++++++++++++++++++++++++++ qpms/qpms_cdefs.pxd | 4 ++++ qpms/scatsystem.c | 19 +++++++++++++------ qpms/scatsystem.h | 41 +++++++++++++++++++++++++++++++---------- 4 files changed, 86 insertions(+), 16 deletions(-) diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx index 3adb2c9..26a0ad4 100644 --- a/qpms/qpms_c.pyx +++ b/qpms/qpms_c.pyx @@ -710,6 +710,25 @@ cdef class ScatteringSystem: return retdict + def scattered_E(self, cdouble wavenumber, scatcoeffvector_full, evalpos): + evalpos = np.array(evalpos, copy=False) + cdef np.ndarray evalpos_a = evalpos + cdef np.ndarray[dtype=complex] scv = np.array(scatcoeffvector_full, copy=False) + cdef cdouble[::1] scv_view = scv + if evalpos.shape[-1] != 3: + raise ValueError("Last dimension of evalpos has to be 3") + cdef np.ndarray[complex] results = np.empty(evalpos.shape, dtype=complex) + cdef ccart3_t res + cdef cart3_t pos + for i in np.ndindex(evalpos.shape[:-1]): + pos.x = evalpos_a[i,0] + pos.y = evalpos_a[i,1] + pos.z = evalpos_a[i,2] + res = qpms_scatsys_scattered_E(self.s, wavenumber, &scv_view[0], pos) + results[i,0] = res.x + results[i,1] = res.y + results[i,2] = res.z + return results def empty_lattice_modes_xy(EpsMu epsmu, reciprocal_basis, wavevector, double maxomega): '''Empty (2D, xy-plane) lattice mode (diffraction order) frequencies of a non-dispersive medium. @@ -844,6 +863,25 @@ 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): + evalpos = np.array(evalpos, copy=False) + cdef np.ndarray evalpos_a = evalpos + cdef np.ndarray[dtype=complex] scv = np.array(scatcoeffvector_full, copy=False) + cdef cdouble[::1] scv_view = scv + if evalpos.shape[-1] != 3: + raise ValueError("Last dimension of evalpos has to be 3") + cdef np.ndarray[complex] results = np.empty(evalpos.shape, dtype=complex) + cdef ccart3_t res + cdef cart3_t pos + for i in np.ndindex(evalpos.shape[:-1]): + pos.x = evalpos_a[i,0] + pos.y = evalpos_a[i,1] + pos.z = evalpos_a[i,2] + res = qpms_scatsysw_scattered_E(self.ssw, &scv_view[0], pos) + results[i,0] = res.x + results[i,1] = res.y + results[i,2] = res.z + return results cdef class ScatteringMatrix: diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index d5c515d..61ea266 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -625,6 +625,10 @@ 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, + const cdouble *f_excitation_vector_full, cart3_t where) + ccart3_t qpms_scatsysw_scattered_E(const qpms_scatsys_at_omega_t *ssw, + const cdouble *f_excitation_vector_full, cart3_t where) cdef extern from "ewald.h": struct qpms_csf_result: diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index 907e779..12cf040 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -1986,11 +1986,11 @@ complex double *qpms_scatsysw_apply_Tmatrices_full( return target_full; } - - -ccart3_t qpms_scatsys_eval_E(const qpms_scatsys_t *ss, - const complex double *cvf, const cart3_t where, - const complex double k) { +ccart3_t qpms_scatsys_scattered_E(const qpms_scatsys_t *ss, + const complex double k, + const complex double *cvf, + const cart3_t where + ) { QPMS_UNTESTED; ccart3_t res = {0,0,0}; ccart3_t res_kc = {0,0,0}; // kahan sum compensation @@ -2012,8 +2012,15 @@ ccart3_t qpms_scatsys_eval_E(const qpms_scatsys_t *ss, return res; } +ccart3_t qpms_scatsysw_scattered_E(const qpms_scatsys_at_omega_t *ssw, + const complex double *cvf, const cart3_t where) { + return qpms_scatsys_scattered_E(ssw->ss, ssw->wavenumber, + cvf, where); +} + + #if 0 -ccart3_t qpms_scatsys_eval_E_irrep(const qpms_scatsys_t *ss, +ccart3_t qpms_scatsys_scattered_E_irrep(const qpms_scatsys_t *ss, qpms_iri_t iri, const complex double *cvr, cart3_t where) { TODO; } diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index a83b1ad..1dbf952 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -240,7 +240,7 @@ typedef struct qpms_scatsys_at_omega_t { qpms_tmatrix_t **tm; complex double omega; ///< Angular frequency qpms_epsmu_t medium; ///< Background medium optical properties at the given frequency - complex double wavenumber; ///< Background medium wavenumber + complex double wavenumber; ///< Background medium wave number } qpms_scatsys_at_omega_t; @@ -692,19 +692,40 @@ complex double *qpms_scatsys_incident_field_vector_irrep_packed( ); #endif -/// Evaluates scattered fields (corresponding to a given excitation vector) at a given point. +/// Evaluates scattered electric field at a point, given a full vector of scattered field coefficients. /** - * By scattered field, one assumes a linear combination of positive-Hankel-type - * spherical waves. + * This function evaluates the field \f$ \vect E (\vect r ) \f$, with any given wavenumber of the + * background medium and any given vector of the excitation coefficients \f$ \wckcout \f$. * - * \return Complex electric field at the point defined by \a where. + * \return Complex electric field at the point defined by \a evalpoint. + * + * \see qpms_scatsysw_scattered_E() + * + * Not yet implemented for periodic systems. */ -ccart3_t qpms_scatsys_eval_E(const qpms_scatsys_t *ss, - const complex double *coeff_vector, ///< A full-length excitation vector (outgoing wave coefficients). - cart3_t where, ///< Evaluation point. - complex double k ///< Wave number. +ccart3_t qpms_scatsys_scattered_E( + const qpms_scatsys_t *ss, + 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. ); +/// Evaluates scattered electric field at a point, given a full vector of scattered field coefficients. +/** + * This function evaluates the field \f$ \vect E (\vect r ) \f$, with any + * given vector of the excitation coefficients \f$ \wckcout \f$. + * + * \return Complex electric field at the point defined by \a evalpoint. + * + * \see qpms_scatsys_scattered_E() + * + * Not yet implemented for periodic systems. + */ +ccart3_t qpms_scatsysw_scattered_E( + const qpms_scatsys_at_omega_t *ssw, + 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. + ); #if 0 /** Evaluates partial scattered fields (corresponding to a given irrep-reduced excitation vector) @@ -712,7 +733,7 @@ ccart3_t qpms_scatsys_eval_E(const qpms_scatsys_t *ss, * * \return Complex electric field at the point defined by \a where. */ -ccart3_t qpms_scatsys_eval_E_irrep(const qpms_scatsys_t *ss, +ccart3_t qpms_scatsys_scattered_E_irrep(const qpms_scatsys_t *ss, qpms_iri_t iri, ///< Irreducible representation const complex double *coeff_vector, ///< A reduced excitation vector, corresponding to \a iri. cart3_t where, ///< Evaluation point.