From 28f4e7f3d2b695ae8632f9c9cf4794e9e1996510 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Wed, 6 Nov 2019 18:13:50 +0200 Subject: [PATCH] Complexify wavenumber in (finite) translation operators. Former-commit-id: f84e1588cdda916d8feda6d807c12bca69512e5f --- qpms/cytranslations.pyx | 13 ++++++------ qpms/qpms_cdefs.pxd | 30 +++++++++++++-------------- qpms/scatsystem.c | 22 ++++++++++---------- qpms/scatsystem.h | 18 ++++++++-------- qpms/translations.c | 46 ++++++++++++++++++++--------------------- qpms/translations.h | 24 ++++++++++----------- qpms/vectors.h | 9 ++++++++ 7 files changed, 86 insertions(+), 76 deletions(-) diff --git a/qpms/cytranslations.pyx b/qpms/cytranslations.pyx index 801346a..3baa551 100644 --- a/qpms/cytranslations.pyx +++ b/qpms/cytranslations.pyx @@ -329,14 +329,14 @@ cdef class trans_calculator: # have to be cdef public in order that __init__ can set these attributes object get_A, get_B, get_AB - def __cinit__(self, int lMax, int normalization = 1): + def __cinit__(self, int lMax, int normalization = QPMS_NORMALISATION_DEFAULT): if (lMax <= 0): raise ValueError('lMax has to be greater than 0.') self.c = qpms_trans_calculator_init(lMax, normalization) if self.c is NULL: raise MemoryError - def __init__(self, int lMax, int normalization = 1): + def __init__(self, int lMax, int normalization = QPMS_NORMALISATION_DEFAULT): if self.c is NULL: raise MemoryError() self.get_A_data[0].c = self.c @@ -401,6 +401,7 @@ cdef class trans_calculator: def nelem(self): return self.c[0].nelem + # THIS FUNCTION MIGHT BE OBSOLETE; NOT SURE WHETHER IT'S WORKING ANY MORE def get_AB_arrays(self, r, theta, phi, r_ge_d, int J, destaxis=None, srcaxis=None, expand=True): """ @@ -520,17 +521,17 @@ cdef class trans_calculator: cdef np.ndarray[np.complex_t, ndim=2] target = np.empty( (destn, srcn), dtype=complex, order='C') cdef cdouble[:,::1] target_view = target - cdef sph_t kdlj_sph + cdef csph_t kdlj_sph kdlj_sph.r = kdlj[0] - kdlj_sph.theta = kdlj[1] - kdlj_sph.phi = kdlj[2] + kdlj_sph.theta = kdlj[1].real + kdlj_sph.phi = kdlj[2].real qpms_trans_calculator_get_trans_array(self.c, &target_view[0][0], destspec.rawpointer(), srcn, srcspec.rawpointer(), 1, kdlj_sph, False, J) return target def get_trans_array_bspec_c3pos(self, BaseSpec destspec, BaseSpec srcspec, - double k, destpos, srcpos, qpms_bessel_t J = QPMS_HANKEL_PLUS): + cdouble k, destpos, srcpos, qpms_bessel_t J = QPMS_HANKEL_PLUS): destpos = np.array(destpos) srcpos = np.array(srcpos) if destpos.shape != (3,) or srcpos.shape != (3,): diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index cf573fa..10ce555 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -313,19 +313,19 @@ cdef extern from "translations.h": qpms_trans_calculator* qpms_trans_calculator_init(int lMax, int nt) # should be qpms_normalization_t void qpms_trans_calculator_free(qpms_trans_calculator* c) cdouble qpms_trans_calculator_get_A_ext(const qpms_trans_calculator* c, - int m, int n, int mu, int nu, double kdlj_r, double kdlj_th, double kdlj_phi, + int m, int n, int mu, int nu, cdouble kdlj_r, double kdlj_th, double kdlj_phi, int r_ge_d, int J) nogil cdouble qpms_trans_calculator_get_B_ext(const qpms_trans_calculator* c, - int m, int n, int mu, int nu, double kdlj_r, double kdlj_th, double kdlj_phi, + int m, int n, int mu, int nu, cdouble kdlj_r, double kdlj_th, double kdlj_phi, int r_ge_d, int J) nogil int qpms_trans_calculator_get_AB_p_ext(const qpms_trans_calculator* c, cdouble *Adest, cdouble *Bdest, - int m, int n, int mu, int nu, double kdlj_r, double kdlj_th, double kdlj_phi, + int m, int n, int mu, int nu, cdouble kdlj_r, double kdlj_th, double kdlj_phi, int r_ge_d, int J) nogil int qpms_trans_calculator_get_AB_arrays_ext(const qpms_trans_calculator *c, cdouble *Adest, cdouble *Bdest, size_t deststride, size_t srcstride, - double kdlj_r, double kdlj_theta, double kdlj_phi, + cdouble kdlj_r, double kdlj_theta, double kdlj_phi, int r_ge_d, int J) nogil int qpms_cython_trans_calculator_get_AB_arrays_loop(qpms_trans_calculator *c, int J, int resnd, @@ -341,14 +341,14 @@ cdef extern from "translations.h": cdouble *target, const qpms_vswf_set_spec_t *destspec, size_t deststride, const qpms_vswf_set_spec_t *srcspec, size_t srcstride, - sph_t kdlj, bint r_ge_d, qpms_bessel_t J); + csph_t kdlj, bint r_ge_d, qpms_bessel_t J); int qpms_trans_calculator_get_trans_array_lc3p( const qpms_trans_calculator *c, cdouble *target, const qpms_vswf_set_spec_t *destspec, size_t deststride, const qpms_vswf_set_spec_t *srcspec, size_t srcstride, - double k, cart3_t destpos, cart3_t srcpos, + cdouble k, cart3_t destpos, cart3_t srcpos, qpms_bessel_t J ); @@ -520,19 +520,19 @@ cdef extern from "scatsystem.h": cdouble *qpms_scatsys_irrep_unpack_vector(cdouble *target_full, 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) + const qpms_scatsys_t *ss, cdouble k) cdouble *qpms_scatsys_build_translation_matrix_full(cdouble *target, - const qpms_scatsys_t *ss, double k) + const qpms_scatsys_t *ss, cdouble k) cdouble *qpms_scatsys_build_translation_matrix_e_full(cdouble *target, - const qpms_scatsys_t *ss, double k, qpms_bessel_t J) + const qpms_scatsys_t *ss, cdouble 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) + const qpms_scatsys_t *ss, qpms_iri_t iri, cdouble k) cdouble *qpms_scatsys_build_translation_matrix_e_irrep_packed(cdouble *target, - const qpms_scatsys_t *ss, qpms_iri_t iri, double k, qpms_bessel_t J) + const qpms_scatsys_t *ss, qpms_iri_t iri, cdouble k, qpms_bessel_t J) cdouble *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorderR( - cdouble *target, const qpms_scatsys_t *ss, qpms_iri_t iri, double k) + cdouble *target, const qpms_scatsys_t *ss, qpms_iri_t iri, cdouble k) cdouble *qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR( - cdouble *target, const qpms_scatsys_t *ss, qpms_iri_t iri, double k) nogil + cdouble *target, const qpms_scatsys_t *ss, qpms_iri_t iri, cdouble k) nogil 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) @@ -546,9 +546,9 @@ cdef extern from "scatsystem.h": int *ipiv void qpms_ss_LU_free(qpms_ss_LU lu) qpms_ss_LU qpms_scatsys_build_modeproblem_matrix_full_LU(cdouble *target, - int *target_piv, const qpms_scatsys_t *ss, double k) + int *target_piv, const qpms_scatsys_t *ss, cdouble k) qpms_ss_LU qpms_scatsys_build_modeproblem_matrix_irrep_packed_LU(cdouble *target, - int *target_piv, const qpms_scatsys_t *ss, qpms_iri_t iri, double k) + int *target_piv, const qpms_scatsys_t *ss, qpms_iri_t iri, cdouble k) qpms_ss_LU qpms_scatsys_modeproblem_matrix_full_factorise(cdouble *modeproblem_matrix_full, int *target_piv, const qpms_scatsys_t *ss) qpms_ss_LU qpms_scatsys_modeproblem_matrix_irrep_packed_factorise(cdouble *modeproblem_matrix_full, diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index de60f77..4c5c8aa 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -967,7 +967,7 @@ 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 k ///< Wave number to use in the translation matrix. ) { return qpms_scatsys_build_translation_matrix_e_full( @@ -978,7 +978,7 @@ complex double *qpms_scatsys_build_translation_matrix_e_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 k, ///< Wave number to use in the translation matrix. qpms_bessel_t J ///< Bessel function type. ) { @@ -1018,7 +1018,7 @@ 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, - double k ///< Wave number to use in the translation matrix. + complex double k ///< Wave number to use in the translation matrix. ) { const size_t full_len = ss->fecv_size; @@ -1067,7 +1067,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed( /// Target memory with capacity for ss->saecv_sizes[iri]**2 elements. If NULL, new will be allocated. complex double *target_packed, const qpms_scatsys_t *ss, qpms_iri_t iri, - double k ///< Wave number to use in the translation matrix. + complex double k ///< Wave number to use in the translation matrix. ) { const size_t packedlen = ss->saecv_sizes[iri]; @@ -1177,7 +1177,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorderR( /// Target memory with capacity for ss->saecv_sizes[iri]**2 elements. If NULL, new will be allocated. complex double *target_packed, const qpms_scatsys_t *ss, qpms_iri_t iri, - double k ///< Wave number to use in the translation matrix. + complex double k ///< Wave number to use in the translation matrix. ) { const size_t packedlen = ss->saecv_sizes[iri]; @@ -1299,7 +1299,7 @@ struct qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR_thread_arg{ pthread_mutex_t *opistartR_mutex; qpms_iri_t iri; complex double *target_packed; - double k; + complex double k; }; static void *qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR_thread(void *arg) @@ -1436,7 +1436,7 @@ struct qpms_scatsys_build_translation_matrix_e_irrep_packed_parallelR_thread_arg pthread_mutex_t *opistartR_mutex; qpms_iri_t iri; complex double *target_packed; - double k; + complex double k; qpms_bessel_t J; }; @@ -1566,7 +1566,7 @@ complex double *qpms_scatsys_build_translation_matrix_e_irrep_packed( complex double *target_packed, const qpms_scatsys_t *ss, qpms_iri_t iri, - double k, ///< Wave number to use in the translation matrix. + complex double k, ///< Wave number to use in the translation matrix. qpms_bessel_t J ) { @@ -1619,7 +1619,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR( /// Target memory with capacity for ss->saecv_sizes[iri]**2 elements. If NULL, new will be allocated. complex double *target_packed, const qpms_scatsys_t *ss, qpms_iri_t iri, - double k ///< Wave number to use in the translation matrix. + complex double k ///< Wave number to use in the translation matrix. ) { const size_t packedlen = ss->saecv_sizes[iri]; @@ -1781,14 +1781,14 @@ qpms_ss_LU qpms_scatsys_modeproblem_matrix_irrep_packed_factorise(complex double qpms_ss_LU qpms_scatsys_build_modeproblem_matrix_full_LU( complex double *target, int *target_piv, - const qpms_scatsys_t *ss, double k){ + const qpms_scatsys_t *ss, complex double k){ target = qpms_scatsys_build_modeproblem_matrix_full(target, ss, k); return qpms_scatsys_modeproblem_matrix_full_factorise(target, target_piv, ss); } qpms_ss_LU qpms_scatsys_build_modeproblem_matrix_irrep_packed_LU( complex double *target, int *target_piv, - const qpms_scatsys_t *ss, qpms_iri_t iri, double k){ + const qpms_scatsys_t *ss, qpms_iri_t iri, complex double k){ target = qpms_scatsys_build_modeproblem_matrix_irrep_packed(target, ss, iri, k); return qpms_scatsys_modeproblem_matrix_irrep_packed_factorise(target, target_piv, ss, iri); } diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index d55002a..c9a9a91 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -237,7 +237,7 @@ 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 k ///< Wave number to use in the translation matrix. ); /// As qpms_scatsys_build_translation_full() but with choice of Bessel function type. @@ -247,7 +247,7 @@ complex double *qpms_scatsys_build_translation_matrix_e_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 k, ///< Wave number to use in the translation matrix. qpms_bessel_t J ); @@ -261,7 +261,7 @@ complex double *qpms_scatsys_build_translation_matrix_e_irrep_packed( complex double *target, const qpms_scatsys_t *ss, qpms_iri_t iri, - double k, ///< Wave number to use in the translation matrix. + complex double k, ///< Wave number to use in the translation matrix. qpms_bessel_t J ); @@ -270,28 +270,28 @@ 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. + complex double k ///< Wave number to use in the translation matrix. ); /// 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, const qpms_scatsys_t *ss, qpms_iri_t iri, - /*COMPLEXIFY*/double k ///< Wave number to use in the translation matrix. + complex double k ///< Wave number to use in the translation matrix. ); /// Alternative implementation of qpms_scatsys_build_modeproblem_matrix_irrep_packed(). complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorderR( /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated. complex double *target, const qpms_scatsys_t *ss, qpms_iri_t iri, - /*COMPLEXIFY*/double k ///< Wave number to use in the translation matrix. + complex double k ///< Wave number to use in the translation matrix. ); /// Alternative implementation of qpms_scatsys_build_modeproblem_matrix_irrep_packed(). complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorder_parallelR( /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated. complex double *target, const qpms_scatsys_t *ss, qpms_iri_t iri, - /*COMPLEXIFY*/double k ///< Wave number to use in the translation matrix. + complex double k ///< Wave number to use in the translation matrix. ); /// LU factorisation (LAPACKE_zgetrf) result holder. @@ -311,7 +311,7 @@ qpms_ss_LU qpms_scatsys_build_modeproblem_matrix_full_LU( complex double *target, ///< Pre-allocated target array. Optional (if NULL, new one is allocated). int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated). const qpms_scatsys_t *ss, - /*COMPLEXIFY*/ double k ///< Wave number to use in the translation matrix. + complex double k ///< Wave number to use in the translation matrix. ); /// Builds an irrep-packed LU-factorised mode/scattering problem matrix from scratch. @@ -319,7 +319,7 @@ qpms_ss_LU qpms_scatsys_build_modeproblem_matrix_irrep_packed_LU( complex double *target, ///< Pre-allocated target array. Optional (if NULL, new one is allocated). int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated). const qpms_scatsys_t *ss, qpms_iri_t iri, - /*COMPLEXIFY*/ double k ///< Wave number to use in the translation matrix. + complex double k ///< Wave number to use in the translation matrix. ); /// Computes LU factorisation of a pre-calculated mode/scattering problem matrix, replacing its contents. diff --git a/qpms/translations.c b/qpms/translations.c index bdc97ce..65b0de3 100644 --- a/qpms/translations.c +++ b/qpms/translations.c @@ -132,7 +132,7 @@ static inline double qpms_trans_normfac(qpms_normalisation_t norm, } complex double qpms_trans_single_A(qpms_normalisation_t norm, - int m, int n, int mu, int nu, sph_t kdlj, + int m, int n, int mu, int nu, csph_t kdlj, bool r_ge_d, qpms_bessel_t J) { TROPS_ONLY_EIMF_IMPLEMENTED(norm); if(r_ge_d) J = QPMS_BESSEL_REGULAR; @@ -201,7 +201,7 @@ complex double qpms_trans_single_A(qpms_normalisation_t norm, complex double qpms_trans_single_B(qpms_normalisation_t norm, - int m, int n, int mu, int nu, sph_t kdlj, + int m, int n, int mu, int nu, csph_t kdlj, bool r_ge_d, qpms_bessel_t J) { TROPS_ONLY_EIMF_IMPLEMENTED(norm); if(r_ge_d) J = QPMS_BESSEL_REGULAR; @@ -509,7 +509,7 @@ qpms_trans_calculator } static inline complex double qpms_trans_calculator_get_A_precalcbuf(const qpms_trans_calculator *c, - int m, int n, int mu, int nu, sph_t kdlj, + int m, int n, int mu, int nu, csph_t kdlj, bool r_ge_d, qpms_bessel_t J, const complex double *bessel_buf, const double *legendre_buf) { TROPS_ONLY_EIMF_IMPLEMENTED(c->normalisation); @@ -530,7 +530,7 @@ static inline complex double qpms_trans_calculator_get_A_precalcbuf(const qpms_t } complex double qpms_trans_calculator_get_A_buf(const qpms_trans_calculator *c, - int m, int n, int mu, int nu, sph_t kdlj, + int m, int n, int mu, int nu, csph_t kdlj, bool r_ge_d, qpms_bessel_t J, complex double *bessel_buf, double *legendre_buf) { // This functions gets preallocated memory for bessel and legendre functions, but computes them itself @@ -549,7 +549,7 @@ complex double qpms_trans_calculator_get_A_buf(const qpms_trans_calculator *c, } static inline complex double qpms_trans_calculator_get_B_precalcbuf(const qpms_trans_calculator *c, - int m, int n, int mu, int nu, sph_t kdlj, + int m, int n, int mu, int nu, csph_t kdlj, bool r_ge_d, qpms_bessel_t J, const complex double *bessel_buf, const double *legendre_buf) { TROPS_ONLY_EIMF_IMPLEMENTED(c->normalisation); @@ -570,7 +570,7 @@ static inline complex double qpms_trans_calculator_get_B_precalcbuf(const qpms_t } complex double qpms_trans_calculator_get_B_buf(const qpms_trans_calculator *c, - int m, int n, int mu, int nu, sph_t kdlj, + int m, int n, int mu, int nu, csph_t kdlj, bool r_ge_d, qpms_bessel_t J, complex double *bessel_buf, double *legendre_buf) { // This functions gets preallocated memory for bessel and legendre functions, but computes them itself @@ -589,7 +589,7 @@ complex double qpms_trans_calculator_get_B_buf(const qpms_trans_calculator *c, int qpms_trans_calculator_get_AB_buf_p(const qpms_trans_calculator *c, complex double *Adest, complex double *Bdest, - int m, int n, int mu, int nu, sph_t kdlj, + int m, int n, int mu, int nu, csph_t kdlj, bool r_ge_d, qpms_bessel_t J, complex double *bessel_buf, double *legendre_buf) { if (r_ge_d) J = QPMS_BESSEL_REGULAR; @@ -614,7 +614,7 @@ int qpms_trans_calculator_get_AB_buf_p(const qpms_trans_calculator *c, int qpms_trans_calculator_get_AB_arrays_buf(const qpms_trans_calculator *c, complex double *Adest, complex double *Bdest, size_t deststride, size_t srcstride, - sph_t kdlj, bool r_ge_d, qpms_bessel_t J, + csph_t kdlj, bool r_ge_d, qpms_bessel_t J, complex double *bessel_buf, double *legendre_buf) { if (r_ge_d) J = QPMS_BESSEL_REGULAR; if (0 == kdlj.r && J != QPMS_BESSEL_REGULAR) { @@ -654,7 +654,7 @@ int qpms_trans_calculator_get_AB_arrays_buf(const qpms_trans_calculator *c, } complex double qpms_trans_calculator_get_A(const qpms_trans_calculator *c, - int m, int n, int mu, int nu, sph_t kdlj, + int m, int n, int mu, int nu, csph_t kdlj, bool r_ge_d, qpms_bessel_t J) { double leg[gsl_sf_legendre_array_n(n+nu)]; complex double bes[n+nu+1]; // maximum order is 2n for A coeffs, plus the zeroth. @@ -663,7 +663,7 @@ complex double qpms_trans_calculator_get_A(const qpms_trans_calculator *c, } complex double qpms_trans_calculator_get_B(const qpms_trans_calculator *c, - int m, int n, int mu, int nu, sph_t kdlj, + int m, int n, int mu, int nu, csph_t kdlj, bool r_ge_d, qpms_bessel_t J) { double leg[gsl_sf_legendre_array_n(n+nu+1)]; complex double bes[n+nu+2]; // maximum order is 2n+1 for B coeffs, plus the zeroth. @@ -673,7 +673,7 @@ complex double qpms_trans_calculator_get_B(const qpms_trans_calculator *c, int qpms_trans_calculator_get_AB_p(const qpms_trans_calculator *c, complex double *Adest, complex double *Bdest, - int m, int n, int mu, int nu, sph_t kdlj, + int m, int n, int mu, int nu, csph_t kdlj, bool r_ge_d, qpms_bessel_t J) { double leg[gsl_sf_legendre_array_n(2*c->lMax+1)]; complex double bes[2*c->lMax+2]; // maximum order is 2n+1 for B coeffs, plus the zeroth. @@ -684,7 +684,7 @@ int qpms_trans_calculator_get_AB_p(const qpms_trans_calculator *c, int qpms_trans_calculator_get_AB_arrays(const qpms_trans_calculator *c, complex double *Adest, complex double *Bdest, size_t deststride, size_t srcstride, - sph_t kdlj, bool r_ge_d, qpms_bessel_t J) { + csph_t kdlj, bool r_ge_d, qpms_bessel_t J) { double leg[gsl_sf_legendre_array_n(c->lMax+c->lMax+1)]; complex double bes[2*c->lMax+2]; // maximum order is 2n+1 for B coeffs, plus the zeroth. return qpms_trans_calculator_get_AB_arrays_buf(c, @@ -700,7 +700,7 @@ qpms_errno_t qpms_trans_calculator_get_trans_array(const qpms_trans_calculator * 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, - sph_t kdlj, bool r_ge_d, qpms_bessel_t J) + csph_t kdlj, bool r_ge_d, qpms_bessel_t J) { TROPS_ONLY_AB_SYMMETRIC_NORMS_IMPLEMENTED(c->normalisation); assert(c->normalisation == destspec->norm && c->normalisation == srcspec->norm); @@ -807,11 +807,11 @@ 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, qpms_bessel_t J + complex double k, cart3_t destpos, cart3_t srcpos, qpms_bessel_t J /// Workspace has to be at least 2 * c->neleme**2 long ) { - sph_t kdlj = cart2sph(cart3_substract(destpos, srcpos)); + csph_t kdlj = cart2csph(cart3_substract(destpos, srcpos)); kdlj.r *= k; return qpms_trans_calculator_get_trans_array(c, target, destspec, deststride, srcspec, srcstride, kdlj, @@ -1098,35 +1098,35 @@ int qpms_trans_calculator_get_AB_arrays_e32(const qpms_trans_calculator *c, complex double qpms_trans_calculator_get_A_ext(const qpms_trans_calculator *c, int m, int n, int mu, int nu, - double kdlj_r, double kdlj_theta, double kdlj_phi, + complex double kdlj_r, double kdlj_theta, double kdlj_phi, int r_ge_d, int J) { - sph_t kdlj = {kdlj_r, kdlj_theta, kdlj_phi}; + csph_t kdlj = {kdlj_r, kdlj_theta, kdlj_phi}; return qpms_trans_calculator_get_A(c,m,n,mu,nu,kdlj,r_ge_d,J); } complex double qpms_trans_calculator_get_B_ext(const qpms_trans_calculator *c, int m, int n, int mu, int nu, - double kdlj_r, double kdlj_theta, double kdlj_phi, + complex double kdlj_r, double kdlj_theta, double kdlj_phi, int r_ge_d, int J) { - sph_t kdlj = {kdlj_r, kdlj_theta, kdlj_phi}; + csph_t kdlj = {kdlj_r, kdlj_theta, kdlj_phi}; return qpms_trans_calculator_get_B(c,m,n,mu,nu,kdlj,r_ge_d,J); } int qpms_trans_calculator_get_AB_p_ext(const qpms_trans_calculator *c, complex double *Adest, complex double *Bdest, int m, int n, int mu, int nu, - double kdlj_r, double kdlj_theta, double kdlj_phi, + complex double kdlj_r, double kdlj_theta, double kdlj_phi, int r_ge_d, int J) { - sph_t kdlj = {kdlj_r, kdlj_theta, kdlj_phi}; + csph_t kdlj = {kdlj_r, kdlj_theta, kdlj_phi}; return qpms_trans_calculator_get_AB_p(c,Adest,Bdest,m,n,mu,nu,kdlj,r_ge_d,J); } int qpms_trans_calculator_get_AB_arrays_ext(const qpms_trans_calculator *c, complex double *Adest, complex double *Bdest, size_t deststride, size_t srcstride, - double kdlj_r, double kdlj_theta, double kdlj_phi, + complex double kdlj_r, double kdlj_theta, double kdlj_phi, int r_ge_d, int J) { - sph_t kdlj = {kdlj_r, kdlj_theta, kdlj_phi}; + csph_t kdlj = {kdlj_r, kdlj_theta, kdlj_phi}; return qpms_trans_calculator_get_AB_arrays(c,Adest,Bdest,deststride,srcstride, kdlj, r_ge_d, J); } diff --git a/qpms/translations.h b/qpms/translations.h index c5ed6e1..7d1785a 100644 --- a/qpms/translations.h +++ b/qpms/translations.h @@ -37,10 +37,10 @@ #endif -complex double qpms_trans_single_A(qpms_normalisation_t norm, qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, sph_t kdlj, +complex double qpms_trans_single_A(qpms_normalisation_t norm, qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, csph_t kdlj, bool r_ge_d, qpms_bessel_t J); -complex double qpms_trans_single_B(qpms_normalisation_t norm, qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, sph_t kdlj, +complex double qpms_trans_single_B(qpms_normalisation_t norm, qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, csph_t kdlj, bool r_ge_d, qpms_bessel_t J); /// Structure holding the constant factors in normalisation operators. @@ -91,38 +91,38 @@ qpms_trans_calculator *qpms_trans_calculator_init(qpms_l_t lMax, ///< Truncation void qpms_trans_calculator_free(qpms_trans_calculator *); complex double qpms_trans_calculator_get_A(const qpms_trans_calculator *c, - qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, sph_t kdlj, + qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, csph_t kdlj, bool r_ge_d, qpms_bessel_t J); complex double qpms_trans_calculator_get_B(const qpms_trans_calculator *c, - qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, sph_t kdlj, + qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, csph_t kdlj, bool r_ge_d, qpms_bessel_t J); int qpms_trans_calculator_get_AB_p(const qpms_trans_calculator *c, complex double *Adest, complex double *Bdest, - qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, sph_t kdlj, + qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, csph_t kdlj, bool r_ge_d, qpms_bessel_t J); int qpms_trans_calculator_get_AB_arrays(const qpms_trans_calculator *c, complex double *Adest, complex double *Bdest, size_t deststride, size_t srcstride, - sph_t kdlj, bool r_ge_d, qpms_bessel_t J); + csph_t kdlj, bool r_ge_d, qpms_bessel_t J); // TODO update the types later complex double qpms_trans_calculator_get_A_ext(const qpms_trans_calculator *c, - int m, int n, int mu, int nu, double kdlj_r, + int m, int n, int mu, int nu, complex double kdlj_r, double kdlj_th, double kdlj_phi, int r_ge_d, int J); complex double qpms_trans_calculator_get_B_ext(const qpms_trans_calculator *c, - int m, int n, int mu, int nu, double kdlj_r, + int m, int n, int mu, int nu, complex double kdlj_r, double kdlj_th, double kdlj_phi, int r_ge_d, int J); int qpms_trans_calculator_get_AB_p_ext(const qpms_trans_calculator *c, complex double *Adest, complex double *Bdest, - int m, int n, int mu, int nu, double kdlj_r, + int m, int n, int mu, int nu, complex double kdlj_r, double kdlj_th, double kdlj_phi, int r_ge_d, int J); int qpms_trans_calculator_get_AB_arrays_ext(const qpms_trans_calculator *c, complex double *Adest, complex double *Bdest, size_t deststride, size_t srcstride, - double kdlj_r, double kdlj_theta, double kdlj_phi, + complex double kdlj_r, double kdlj_theta, double kdlj_phi, int r_ge_d, int J); // Convenience functions using VSWF base specs @@ -132,7 +132,7 @@ qpms_errno_t qpms_trans_calculator_get_trans_array(const qpms_trans_calculator * 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, - sph_t kdlj, bool r_ge_d, qpms_bessel_t J); + csph_t kdlj, bool r_ge_d, qpms_bessel_t J); /// Version with \a k and cartesian particle positions /// and with automatic \a r_ge_d = `false`. @@ -143,7 +143,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, + complex double k, cart3_t destpos, cart3_t srcpos, qpms_bessel_t J /// Workspace has to be at least 2 * c->neleme**2 long ); diff --git a/qpms/vectors.h b/qpms/vectors.h index e1dd31a..cf4b01c 100644 --- a/qpms/vectors.h +++ b/qpms/vectors.h @@ -333,6 +333,15 @@ static inline csph_t ccart2csph(const ccart3_t cart) { return sph; } +/// Real 3D cartesian to spherical (complex r) coordinates conversion. See @ref coord_conversions. +static inline csph_t cart2csph(const cart3_t cart) { + csph_t sph; + sph.r = cart3norm(cart); + sph.theta = sph.r ? acos(cart.z / sph.r) : M_PI_2; + sph.phi = atan2(cart.y, cart.x); + return sph; +} + /// Coordinate transform of csph_t to ccart3_t static inline ccart3_t csph2ccart(const csph_t sph) { ccart3_t cart;