diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx index af0b511..fe8d5fd 100644 --- a/qpms/qpms_c.pyx +++ b/qpms/qpms_c.pyx @@ -1571,12 +1571,12 @@ cdef class ScatteringSystem: qpms_scatsys_build_modeproblem_matrix_irrep_packed(&target_view[0][0], self.s, iri, k) return target - def translation_matrix_full(self, double k): + def translation_matrix_full(self, double k, J = QPMS_HANKEL_PLUS): cdef size_t flen = self.s[0].fecv_size cdef np.ndarray[np.complex_t, ndim=2] target = np.empty( (flen,flen),dtype=complex, order='C') cdef cdouble[:,::1] target_view = target - qpms_scatsys_build_translation_matrix_full(&target_view[0][0], self.s, k) + qpms_scatsys_build_translation_matrix_full_e(&target_view[0][0], self.s, k, J) return target def fullvec_psizes(self): diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index c03c7dc..1bcbf69 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -367,6 +367,8 @@ cdef extern from "scatsystem.h": const qpms_scatsys_t *ss, double k) cdouble *qpms_scatsys_build_translation_matrix_full(cdouble *target, const qpms_scatsys_t *ss, double k) + cdouble *qpms_scatsys_build_translation_matrix_full_e(cdouble *target, + const qpms_scatsys_t *ss, double k, qpms_bessel_t J) cdouble *qpms_scatsys_build_modeproblem_matrix_irrep_packed(cdouble *target, const qpms_scatsys_t *ss, qpms_iri_t iri, double k) cdouble *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorderR( diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index 6ed93f9..f79e9e0 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -964,6 +964,18 @@ complex double *qpms_scatsys_build_translation_matrix_full( const qpms_scatsys_t *ss, double k ///< Wave number to use in the translation matrix. ) +{ + return qpms_scatsys_build_translation_matrix_full_e( + target, ss, k, QPMS_HANKEL_PLUS); +} + +complex double *qpms_scatsys_build_translation_matrix_full_e( + /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated. + complex double *target, + const qpms_scatsys_t *ss, + double k, ///< Wave number to use in the translation matrix. + qpms_bessel_t J ///< Bessel function type. + ) { const size_t full_len = ss->fecv_size; if(!target) @@ -982,7 +994,7 @@ complex double *qpms_scatsys_build_translation_matrix_full( QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_lc3p(ss->c, target + fullvec_offsetR*full_len + fullvec_offsetC, bspecR, full_len, bspecC, 1, - k, posR, posC)); + k, posR, posC, J)); } fullvec_offsetC += bspecC->n; } @@ -1025,7 +1037,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_full( QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_lc3p(ss->c, tmp, // tmp is S(piR<-piC) bspecR, bspecC->n, bspecC, 1, - k, posR, posC)); + k, posR, posC, QPMS_HANKEL_PLUS)); cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, bspecR->n /*m*/, bspecC->n /*n*/, bspecR->n /*k*/, &one/*alpha*/, tmmR/*a*/, bspecR->n/*lda*/, @@ -1117,7 +1129,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed( QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_lc3p(ss->c, Sblock, // Sblock is S(piR->piC) bspecR, bspecC->n, bspecC, 1, - k, posR, posC)); + k, posR, posC, QPMS_HANKEL_PLUS)); cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, bspecR->n /*m*/, bspecC->n /*n*/, bspecR->n /*k*/, @@ -1236,7 +1248,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorderR( QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_lc3p(ss->c, Sblock, // Sblock is S(piR->piC) bspecR, bspecC->n, bspecC, 1, - k, posR, posC)); + k, posR, posC, QPMS_HANKEL_PLUS)); cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, bspecR->n /*m*/, bspecC->n /*n*/, bspecR->n /*k*/, @@ -1368,7 +1380,7 @@ static void *qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR_thread QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_lc3p(ss->c, Sblock, // Sblock is S(piR->piC) bspecR, bspecC->n, bspecC, 1, - a->k, posR, posC)); + a->k, posR, posC, QPMS_HANKEL_PLUS)); SERIAL_ZGEMM(CblasRowMajor, CblasNoTrans, CblasNoTrans, bspecR->n /*m*/, bspecC->n /*n*/, bspecR->n /*k*/, diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index 945c305..7f5c5b0 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -234,6 +234,18 @@ complex double *qpms_scatsys_build_translation_matrix_full( const qpms_scatsys_t *ss, double k ///< Wave number to use in the translation matrix. ); + +/// As qpms_scatsys_build_translation_full() but with choice of Bessel function type. +/** Might be useful for evaluation of cross sections and testing. + */ +complex double *qpms_scatsys_build_translation_matrix_full_e( + /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated. + complex double *target, + const qpms_scatsys_t *ss, + double k, ///< Wave number to use in the translation matrix. + qpms_bessel_t J + ); + 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, diff --git a/qpms/translations.c b/qpms/translations.c index 852e45e..73b3186 100644 --- a/qpms/translations.c +++ b/qpms/translations.c @@ -1119,8 +1119,6 @@ qpms_errno_t qpms_trans_calculator_get_trans_array(const qpms_trans_calculator * return retval; } -/// Version with \a k and cartesian particle positions -/// and with automatic \a kdlj = false and \a J = QPMS_HANKEL_PLUS. qpms_errno_t qpms_trans_calculator_get_trans_array_lc3p( const qpms_trans_calculator *c, complex double *target, @@ -1128,7 +1126,7 @@ qpms_errno_t qpms_trans_calculator_get_trans_array_lc3p( const qpms_vswf_set_spec_t *destspec, size_t deststride, /// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm. const qpms_vswf_set_spec_t *srcspec, size_t srcstride, - double k, cart3_t destpos, cart3_t srcpos + double k, cart3_t destpos, cart3_t srcpos, qpms_bessel_t J /// Workspace has to be at least 2 * c->neleme**2 long ) { @@ -1136,7 +1134,7 @@ qpms_errno_t qpms_trans_calculator_get_trans_array_lc3p( kdlj.r *= k; return qpms_trans_calculator_get_trans_array(c, target, destspec, deststride, srcspec, srcstride, kdlj, - false, QPMS_HANKEL_PLUS); + false, J); } #ifdef LATTICESUMS31 diff --git a/qpms/translations.h b/qpms/translations.h index 3af105d..f5974ca 100644 --- a/qpms/translations.h +++ b/qpms/translations.h @@ -155,7 +155,7 @@ qpms_errno_t qpms_trans_calculator_get_trans_array(const qpms_trans_calculator * sph_t kdlj, bool r_ge_d, qpms_bessel_t J); /// Version with \a k and cartesian particle positions -/// and with automatic \a kdlj = false and \a J = QPMS_HANKEL_PLUS. +/// and with automatic \a r_ge_d = `false`. qpms_errno_t qpms_trans_calculator_get_trans_array_lc3p( const qpms_trans_calculator *c, complex double *target, @@ -163,7 +163,8 @@ qpms_errno_t qpms_trans_calculator_get_trans_array_lc3p( const qpms_vswf_set_spec_t *destspec, size_t deststride, /// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm. const qpms_vswf_set_spec_t *srcspec, size_t srcstride, - double k, cart3_t destpos, cart3_t srcpos + double k, cart3_t destpos, cart3_t srcpos, + qpms_bessel_t J /// Workspace has to be at least 2 * c->neleme**2 long );