From e05a79483d741ae1fea51e5b3f0f48305bb681cb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Mon, 22 Jul 2019 16:56:03 +0300 Subject: [PATCH] Apply T-matrices on incident fields vector. Former-commit-id: c1cfd5a2c24651f55325d73430387a4e87276d48 --- qpms/qpms_c.pyx | 13 +++++++++++++ qpms/qpms_cdefs.pxd | 2 ++ qpms/scatsystem.c | 19 +++++++++++++++++++ qpms/scatsystem.h | 11 ++++++++++- qpms/tmatrices.c | 13 +++++++++++++ qpms/tmatrices.h | 8 ++++++++ 6 files changed, 65 insertions(+), 1 deletion(-) diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx index b607488..3fc1250 100644 --- a/qpms/qpms_c.pyx +++ b/qpms/qpms_c.pyx @@ -1698,6 +1698,19 @@ cdef class ScatteringSystem: self.s, qpms_incfield_planewave, &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 diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index 64f8b0f..28a758b 100644 --- a/qpms/qpms_cdefs.pxd +++ b/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) diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index 1c7158e..02e04da 100644 --- a/qpms/scatsystem.c +++ b/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) { diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index ba8646b..a3a88f1 100644 --- a/qpms/scatsystem.h +++ b/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! */ diff --git a/qpms/tmatrices.c b/qpms/tmatrices.c index f5bbbb7..7a5210f 100644 --- a/qpms/tmatrices.c +++ b/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; +} + + diff --git a/qpms/tmatrices.h b/qpms/tmatrices.h index b4aa876..fb2b42b 100644 --- a/qpms/tmatrices.h +++ b/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 typedef struct gsl_spline gsl_spline; // Forward declaration for the interpolator struct typedef struct gsl_interp_type gsl_interp_type;