Apply T-matrices on incident fields vector.

Former-commit-id: c1cfd5a2c24651f55325d73430387a4e87276d48

qpms/qpms_c.pyx

@@ -1698,6 +1698,19 @@ cdef class ScatteringSystem:
                 self.s, qpms_incfield_planewave, <void *>&p, 0)
         return target_np
 
+    def apply_Tmatrices_full(self, a):
+        if len(a) != self.fecv_size: 
+            raise ValueError("Length of a full vector has to be %d, not %d" 
+                    % (self.fecv_size, len(a)))
+        a = np.array(a, dtype=complex, copy=False, order='C')
+        cdef cdouble[::1] a_view = a;
+        cdef np.ndarray[np.complex_t, ndim=1] target_np = np.empty(
+                (self.fecv_size,), dtype=complex, order='C')
+        cdef cdouble[::1] target_view = target_np
+        qpms_scatsys_apply_Tmatrices_full(&target_view[0], &a_view[0], self.s)
+        return target_np
+
+
 
 def tlm2uvswfi(t, l, m):
     ''' TODO doc

qpms/qpms_cdefs.pxd

@@ -396,6 +396,8 @@ cdef extern from "scatsystem.h":
     cdouble *qpms_scatsys_incident_field_vector_full(cdouble *target_full,
             const qpms_scatsys_t *ss, qpms_incfield_t field_at_point, 
             const void *args, bint add)
+    cdouble *qpms_scatsys_apply_Tmatrices_full(cdouble *target_full, const cdouble *inc_full,
+            const qpms_scatsys_t *ss)
 
 
 

qpms/scatsystem.c

@@ -1688,6 +1688,25 @@ complex double *qpms_scatsys_incident_field_vector_irrep_packed(
 #endif
 
 
+complex double *qpms_scatsys_apply_Tmatrices_full(
+		complex double *target_full, const complex double *inc_full, 
+		const qpms_scatsys_t *ss) {
+  QPMS_UNTESTED;
+  if (!target_full) QPMS_CRASHING_CALLOC(target_full, ss->fecv_size,
+      sizeof(complex double));
+  for(qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi) {
+    complex double *ptarget = target_full + ss->fecv_pstarts[pi];
+    const complex double *psrc = inc_full + ss->fecv_pstarts[pi];
+    const qpms_vswf_set_spec_t *bspec = qpms_ss_bspec_pi(ss, pi);
+    // TODO check whether T-matrix is non-virtual after virtual t-matrices are implemented.
+    const qpms_tmatrix_t *T = ss->tm[ss->p[pi].tmatrix_id];
+    qpms_apply_tmatrix(ptarget, psrc, T);
+  }
+  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) {

qpms/scatsystem.h

@@ -265,13 +265,14 @@ complex double *qpms_scatsys_build_translation_matrix_e_irrep_packed(
 		qpms_bessel_t J
 		);
 
+/// Creates the full \f$ (-I + TS) \f$ matrix of the scattering system.
 complex double *qpms_scatsys_build_modeproblem_matrix_full(
 		/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
 		complex double *target,
 		const qpms_scatsys_t *ss,
 		/*COMPLEXIFY*/double k ///< Wave number to use in the translation matrix.
 		);
-/// Creates the mode problem matrix directly in the irrep-packed form.
+/// Creates the mode problem matrix \f$ (-I + TS) \f$ directly in the irrep-packed form.
 complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed(
 		/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
 		complex double *target,
@@ -370,6 +371,14 @@ complex double *qpms_scatsys_incident_field_vector_full(
 		bool add ///< If true, add to target_full; rewrite target_full if false.
 		);
 
+/// Applies T-matrices onto an incident field vector in the full basis.
+complex double *qpms_scatsys_apply_Tmatrices_full(
+		complex double *target_full, /// Target vector array. If NULL, a new one is allocated.
+		const complex double *inc_full, /// Incident field coefficient vector. Must not be NULL.
+		const qpms_scatsys_t *ss
+		);
+
+
 #if 0
 /// Creates a (partial) incident field vector in the symmetry-adapted basis, given a function that evaluates the field expansions at points.
 /** TODO detailed doc! */

qpms/tmatrices.c

@@ -704,3 +704,16 @@ qpms_errno_t qpms_tmatrix_spherical_mu0_fill(
   return qpms_tmatrix_spherical_fill(t, a, k_i, k_e, 1, 1);
 }
 
+complex double *qpms_apply_tmatrix(
+    complex double *f, const complex double *a, const qpms_tmatrix_t *T) 
+{
+  const size_t n = T->spec->n;
+  if(!f)
+    QPMS_CRASHING_CALLOC(f, n, sizeof(complex double));
+  const complex double one = 1;
+  const complex double zero = 0;
+  cblas_zgemv(CblasRowMajor, CblasNoTrans, n, n, &one, T->m, n, a, 1, &zero, NULL, 1);
+  return f;
+}
+
+

qpms/tmatrices.h

@@ -228,6 +228,14 @@ qpms_tmatrix_t *qpms_tmatrix_symmetrise_finite_group_inplace(
 		const qpms_finite_group_t *pointgroup
 		);
 
+/// Application of T-matrix on a vector of incident field coefficients, \f$ f = Ta \f$.
+complex double *qpms_apply_tmatrix(
+		complex double *f_target, ///< Scattered field coefficient array of size T->spec->n; if NULL, a new one is allocated.
+		const complex double *a, ///< Incident field coefficient array of size T->spec->n.
+		const qpms_tmatrix_t *T
+		);
+
+
 /* Fuck this, include the whole <gsl/gsl_spline.h>
 typedef struct gsl_spline gsl_spline; // Forward declaration for the interpolator struct
 typedef struct gsl_interp_type gsl_interp_type;