diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx index 09032a6..db78208 100644 --- a/qpms/qpms_c.pyx +++ b/qpms/qpms_c.pyx @@ -1463,6 +1463,15 @@ cdef class ScatteringSystem: qpms_scatsys_build_modeproblem_matrix_full(&target_view[0][0], self.s, k) return target + def translation_matrix_full(self, double k): + 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) + return target + + def tlm2uvswfi(t, l, m): ''' TODO doc diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index a9691ed..80b60e9 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -287,6 +287,8 @@ cdef extern from "scatsystem.h": const cdouble *orig_packed, const qpms_scatsys_t *ss, qpms_iri_t iri, bint add) cdouble *qpms_scatsys_build_modeproblem_matrix_full(cdouble *target, 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_modeproblem_matrix_irrep_packed(cdouble *target, const qpms_scatsys_t *ss, qpms_iri_t iri, double k) diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index 0a38109..8721f3f 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -1032,6 +1032,43 @@ complex double *qpms_scatsys_irrep_unpack_vector(complex double *target_full, return target_full; } +complex double *qpms_scatsys_build_translation_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, + double k ///< Wave number to use in the translation matrix. + ) +{ + const size_t full_len = ss->fecv_size; + if(!target) + QPMS_CRASHING_MALLOC(target, SQ(full_len) * sizeof(complex double)); + memset(target, 0, SQ(full_len) * sizeof(complex double)); //unnecessary? + { // Non-diagonal part; M[piR, piC] = T[piR] S(piR<-piC) + size_t fullvec_offsetR = 0; + for(qpms_ss_pi_t piR = 0; piR < ss->p_count; ++piR) { + const qpms_vswf_set_spec_t *bspecR = ss->tm[ss->p[piR].tmatrix_id]->spec; + const cart3_t posR = ss->p[piR].pos; + size_t fullvec_offsetC = 0; + for(qpms_ss_pi_t piC = 0; piC < ss->p_count; ++piC) { + if(piC == piR) continue; // The diagonal will be dealt with later. + const qpms_vswf_set_spec_t *bspecC = ss->tm[ss->p[piC].tmatrix_id]->spec; + const cart3_t posC = ss->p[piC].pos; + 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)); + fullvec_offsetC += bspecC->n; + } + assert(fullvec_offsetC = full_len); + fullvec_offsetR += bspecR->n; + } + assert(fullvec_offsetR == full_len); + } + + return target; +} + + 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/scatsystem.h b/qpms/scatsystem.h index 31a0e37..be503f3 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -450,6 +450,12 @@ complex double *qpms_scatsys_irrep_unpack_vector(complex double *target_full, const complex double *orig_packed, const qpms_scatsys_t *ss, qpms_iri_t iri, bool add); +complex double *qpms_scatsys_build_translation_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, + double k ///< Wave number to use in the translation matrix. + ); 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/tests/bspectransl.c b/tests/bspectransl.c new file mode 100644 index 0000000..0210b95 --- /dev/null +++ b/tests/bspectransl.c @@ -0,0 +1,42 @@ +// c99 -g -I.. bspectransl.c ../qpms/translations.c ../qpms/legendre.c ../qpms/vswf.c ../qpms/gaunt.c ../qpms/error.c -lm -lgsl -lblas +#include +#include +#include +#include + +int main() { + cart3_t pos1={0,1,2}, pos2={3,5,2}; + qpms_vswf_set_spec_t + *bspec1 = qpms_vswf_set_spec_from_lMax(1, QPMS_NORMALISATION_POWER), + *bspec2 = qpms_vswf_set_spec_from_lMax(2, QPMS_NORMALISATION_POWER); + + qpms_trans_calculator *c = qpms_trans_calculator_init(3, QPMS_NORMALISATION_POWER); + + complex double s_2_1[bspec2->n][bspec1->n]; + complex double s_1_2[bspec1->n][bspec2->n]; + const double k = 1.8; + + qpms_trans_calculator_get_trans_array_lc3p(c, s_2_1[0], bspec2, bspec1->n, + bspec1, 1, k, pos2, pos1); + + qpms_trans_calculator_get_trans_array_lc3p(c, s_1_2[0], bspec1, bspec2->n, + bspec2, 1, k, pos1, pos2); + + for(size_t R = 0; R < bspec2->n; ++R) { + for(size_t C = 0; C < bspec1->n; ++C) + printf("%.3lg+%.3lgj\t", creal(s_2_1[R][C]), cimag(s_2_1[C][R])); + putchar('\n'); + } + putchar('\n'); + + for(size_t R = 0; R < bspec1->n; ++R) { + for(size_t C = 0; C < bspec2->n; ++C) + printf("%.3lg+%.3lgj\t", creal(s_1_2[R][C]), cimag(s_1_2[R][C])); + putchar('\n'); + } + + qpms_trans_calculator_free(c); + qpms_vswf_set_spec_free(bspec1); + qpms_vswf_set_spec_free(bspec2); + return 0; +}