Scatsystem scattered E field methods

Former-commit-id: 6521916a81a7c8022bfaefb7102ae657eef99516

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:

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:

qpms/scatsystem.c

@@ -1986,11 +1986,11 @@ complex double *qpms_scatsysw_apply_Tmatrices_full(
   return target_full;
 }
 
-
+ccart3_t qpms_scatsys_scattered_E(const qpms_scatsys_t *ss, 
-
+    const complex double k,
-ccart3_t qpms_scatsys_eval_E(const qpms_scatsys_t *ss, 
+    const complex double *cvf, 
-    const complex double *cvf, const cart3_t where,
+    const cart3_t where
-    const complex double k) {
+    ) {
   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;
 }

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
+ * This function evaluates the field \f$ \vect E (\vect r ) \f$, with any given wavenumber of the
- * spherical waves.
+ * 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,
+ccart3_t qpms_scatsys_scattered_E(
-		const complex double *coeff_vector, ///< A full-length excitation vector (outgoing wave coefficients).
+		const qpms_scatsys_t *ss,
-		cart3_t where, ///< Evaluation point.
+		complex double wavenumber, ///< Wavenumber of the background medium.
-		complex double k ///< Wave number.
+		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.