QPMS
/
qpms
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Complexify wavenumber in (finite) translation operators. 

Former-commit-id: f84e1588cdda916d8feda6d807c12bca69512e5f
This commit is contained in:
Marek Nečada 2019-11-06 18:13:50 +02:00
parent cc9dbb6cc5
commit 28f4e7f3d2
7 changed files with 86 additions and 76 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
qpms/cytranslations.pyx
View File

@ -329,14 +329,14 @@ cdef class trans_calculator:
# have to be cdef public in order that __init__ can set these attributes # have to be cdef public in order that __init__ can set these attributes
object get_A, get_B, get_AB 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): if (lMax <= 0):
raise ValueError('lMax has to be greater than 0.') raise ValueError('lMax has to be greater than 0.')
self.c = qpms_trans_calculator_init(lMax, normalization) self.c = qpms_trans_calculator_init(lMax, normalization)
if self.c is NULL: if self.c is NULL:
raise MemoryError 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: if self.c is NULL:
raise MemoryError() raise MemoryError()
self.get_A_data[0].c = self.c self.get_A_data[0].c = self.c
@ -401,6 +401,7 @@ cdef class trans_calculator:
def nelem(self): def nelem(self):
return self.c[0].nelem 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, def get_AB_arrays(self, r, theta, phi, r_ge_d, int J,
destaxis=None, srcaxis=None, expand=True): 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( cdef np.ndarray[np.complex_t, ndim=2] target = np.empty(
(destn, srcn), dtype=complex, order='C') (destn, srcn), dtype=complex, order='C')
cdef cdouble[:,::1] target_view = target cdef cdouble[:,::1] target_view = target
cdef sph_t kdlj_sph cdef csph_t kdlj_sph
kdlj_sph.r = kdlj[0] kdlj_sph.r = kdlj[0]
kdlj_sph.theta = kdlj[1] kdlj_sph.theta = kdlj[1].real
kdlj_sph.phi = kdlj[2] kdlj_sph.phi = kdlj[2].real
qpms_trans_calculator_get_trans_array(self.c, &target_view[0][0], qpms_trans_calculator_get_trans_array(self.c, &target_view[0][0],
destspec.rawpointer(), srcn, srcspec.rawpointer(), 1, destspec.rawpointer(), srcn, srcspec.rawpointer(), 1,
kdlj_sph, False, J) kdlj_sph, False, J)
return target return target
def get_trans_array_bspec_c3pos(self, BaseSpec destspec, BaseSpec srcspec, 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) destpos = np.array(destpos)
srcpos = np.array(srcpos) srcpos = np.array(srcpos)
if destpos.shape != (3,) or srcpos.shape != (3,): if destpos.shape != (3,) or srcpos.shape != (3,):

30
qpms/qpms_cdefs.pxd
View File

@ -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 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) void qpms_trans_calculator_free(qpms_trans_calculator* c)
cdouble qpms_trans_calculator_get_A_ext(const 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 int r_ge_d, int J) nogil
cdouble qpms_trans_calculator_get_B_ext(const qpms_trans_calculator* c, 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 r_ge_d, int J) nogil
int qpms_trans_calculator_get_AB_p_ext(const qpms_trans_calculator* c, int qpms_trans_calculator_get_AB_p_ext(const qpms_trans_calculator* c,
cdouble *Adest, cdouble *Bdest, 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 r_ge_d, int J) nogil
int qpms_trans_calculator_get_AB_arrays_ext(const qpms_trans_calculator *c, int qpms_trans_calculator_get_AB_arrays_ext(const qpms_trans_calculator *c,
cdouble *Adest, cdouble *Bdest, cdouble *Adest, cdouble *Bdest,
size_t deststride, size_t srcstride, 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 r_ge_d, int J) nogil
int qpms_cython_trans_calculator_get_AB_arrays_loop(qpms_trans_calculator *c, int qpms_cython_trans_calculator_get_AB_arrays_loop(qpms_trans_calculator *c,
int J, int resnd, int J, int resnd,
@ -341,14 +341,14 @@ cdef extern from "translations.h":
cdouble *target, cdouble *target,
const qpms_vswf_set_spec_t *destspec, size_t deststride, const qpms_vswf_set_spec_t *destspec, size_t deststride,
const qpms_vswf_set_spec_t *srcspec, size_t srcstride, 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( int qpms_trans_calculator_get_trans_array_lc3p(
const qpms_trans_calculator *c, const qpms_trans_calculator *c,
cdouble *target, cdouble *target,
const qpms_vswf_set_spec_t *destspec, size_t deststride, const qpms_vswf_set_spec_t *destspec, size_t deststride,
const qpms_vswf_set_spec_t *srcspec, size_t srcstride, 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 qpms_bessel_t J
); );
@ -520,19 +520,19 @@ cdef extern from "scatsystem.h":
cdouble *qpms_scatsys_irrep_unpack_vector(cdouble *target_full, cdouble *qpms_scatsys_irrep_unpack_vector(cdouble *target_full,
const cdouble *orig_packed, const qpms_scatsys_t *ss, qpms_iri_t iri, bint add) const cdouble *orig_packed, const qpms_scatsys_t *ss, qpms_iri_t iri, bint add)
cdouble *qpms_scatsys_build_modeproblem_matrix_full(cdouble *target, 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, 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, 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, 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, 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 *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 *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, cdouble *qpms_scatsys_incident_field_vector_full(cdouble *target_full,
const qpms_scatsys_t *ss, qpms_incfield_t field_at_point, const qpms_scatsys_t *ss, qpms_incfield_t field_at_point,
const void *args, bint add) const void *args, bint add)
@ -546,9 +546,9 @@ cdef extern from "scatsystem.h":
int *ipiv int *ipiv
void qpms_ss_LU_free(qpms_ss_LU lu) void qpms_ss_LU_free(qpms_ss_LU lu)
qpms_ss_LU qpms_scatsys_build_modeproblem_matrix_full_LU(cdouble *target, 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, 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, qpms_ss_LU qpms_scatsys_modeproblem_matrix_full_factorise(cdouble *modeproblem_matrix_full,
int *target_piv, const qpms_scatsys_t *ss) int *target_piv, const qpms_scatsys_t *ss)
qpms_ss_LU qpms_scatsys_modeproblem_matrix_irrep_packed_factorise(cdouble *modeproblem_matrix_full, qpms_ss_LU qpms_scatsys_modeproblem_matrix_irrep_packed_factorise(cdouble *modeproblem_matrix_full,

22
qpms/scatsystem.c
View File

@ -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. /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
complex double *target, complex double *target,
const qpms_scatsys_t *ss, 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( 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. /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
complex double *target, complex double *target,
const qpms_scatsys_t *ss, 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. 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. /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
complex double *target, complex double *target,
const qpms_scatsys_t *ss, 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; 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. /// Target memory with capacity for ss->saecv_sizes[iri]**2 elements. If NULL, new will be allocated.
complex double *target_packed, complex double *target_packed,
const qpms_scatsys_t *ss, qpms_iri_t iri, 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]; 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. /// Target memory with capacity for ss->saecv_sizes[iri]**2 elements. If NULL, new will be allocated.
complex double *target_packed, complex double *target_packed,
const qpms_scatsys_t *ss, qpms_iri_t iri, 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]; 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; pthread_mutex_t *opistartR_mutex;
qpms_iri_t iri; qpms_iri_t iri;
complex double *target_packed; complex double *target_packed;
double k; complex double k;
}; };
static void *qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR_thread(void *arg) 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; pthread_mutex_t *opistartR_mutex;
qpms_iri_t iri; qpms_iri_t iri;
complex double *target_packed; complex double *target_packed;
double k; complex double k;
qpms_bessel_t J; qpms_bessel_t J;
}; };
@ -1566,7 +1566,7 @@ complex double *qpms_scatsys_build_translation_matrix_e_irrep_packed(
complex double *target_packed, complex double *target_packed,
const qpms_scatsys_t *ss, const qpms_scatsys_t *ss,
qpms_iri_t iri, 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 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. /// Target memory with capacity for ss->saecv_sizes[iri]**2 elements. If NULL, new will be allocated.
complex double *target_packed, complex double *target_packed,
const qpms_scatsys_t *ss, qpms_iri_t iri, 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]; 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( qpms_ss_LU qpms_scatsys_build_modeproblem_matrix_full_LU(
complex double *target, int *target_piv, 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); target = qpms_scatsys_build_modeproblem_matrix_full(target, ss, k);
return qpms_scatsys_modeproblem_matrix_full_factorise(target, target_piv, ss); return qpms_scatsys_modeproblem_matrix_full_factorise(target, target_piv, ss);
} }
qpms_ss_LU qpms_scatsys_build_modeproblem_matrix_irrep_packed_LU( qpms_ss_LU qpms_scatsys_build_modeproblem_matrix_irrep_packed_LU(
complex double *target, int *target_piv, 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); 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); return qpms_scatsys_modeproblem_matrix_irrep_packed_factorise(target, target_piv, ss, iri);
} }

18
qpms/scatsystem.h
View File

@ -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. /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
complex double *target, complex double *target,
const qpms_scatsys_t *ss, 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. /// 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. /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
complex double *target, complex double *target,
const qpms_scatsys_t *ss, 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 qpms_bessel_t J
); );
@ -261,7 +261,7 @@ complex double *qpms_scatsys_build_translation_matrix_e_irrep_packed(
complex double *target, complex double *target,
const qpms_scatsys_t *ss, const qpms_scatsys_t *ss,
qpms_iri_t iri, 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 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. /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
complex double *target, complex double *target,
const qpms_scatsys_t *ss, 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. /// Creates the mode problem matrix \f$ (I - TS) \f$ directly in the irrep-packed form.
complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed( 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. /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
complex double *target, complex double *target,
const qpms_scatsys_t *ss, qpms_iri_t iri, 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(). /// Alternative implementation of qpms_scatsys_build_modeproblem_matrix_irrep_packed().
complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorderR( 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. /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
complex double *target, complex double *target,
const qpms_scatsys_t *ss, qpms_iri_t iri, 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(). /// Alternative implementation of qpms_scatsys_build_modeproblem_matrix_irrep_packed().
complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorder_parallelR( 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. /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
complex double *target, complex double *target,
const qpms_scatsys_t *ss, qpms_iri_t iri, 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. /// 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). 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). int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated).
const qpms_scatsys_t *ss, 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. /// 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). 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). int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated).
const qpms_scatsys_t *ss, qpms_iri_t iri, 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. /// Computes LU factorisation of a pre-calculated mode/scattering problem matrix, replacing its contents.

46
qpms/translations.c
View File

@ -132,7 +132,7 @@ static inline double qpms_trans_normfac(qpms_normalisation_t norm,
} }
complex double qpms_trans_single_A(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) { bool r_ge_d, qpms_bessel_t J) {
TROPS_ONLY_EIMF_IMPLEMENTED(norm); TROPS_ONLY_EIMF_IMPLEMENTED(norm);
if(r_ge_d) J = QPMS_BESSEL_REGULAR; 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, 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) { bool r_ge_d, qpms_bessel_t J) {
TROPS_ONLY_EIMF_IMPLEMENTED(norm); TROPS_ONLY_EIMF_IMPLEMENTED(norm);
if(r_ge_d) J = QPMS_BESSEL_REGULAR; 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, 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, bool r_ge_d, qpms_bessel_t J,
const complex double *bessel_buf, const double *legendre_buf) { const complex double *bessel_buf, const double *legendre_buf) {
TROPS_ONLY_EIMF_IMPLEMENTED(c->normalisation); 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, 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, bool r_ge_d, qpms_bessel_t J,
complex double *bessel_buf, double *legendre_buf) { complex double *bessel_buf, double *legendre_buf) {
// This functions gets preallocated memory for bessel and legendre functions, but computes them itself // 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, 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, bool r_ge_d, qpms_bessel_t J,
const complex double *bessel_buf, const double *legendre_buf) { const complex double *bessel_buf, const double *legendre_buf) {
TROPS_ONLY_EIMF_IMPLEMENTED(c->normalisation); 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, 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, bool r_ge_d, qpms_bessel_t J,
complex double *bessel_buf, double *legendre_buf) { complex double *bessel_buf, double *legendre_buf) {
// This functions gets preallocated memory for bessel and legendre functions, but computes them itself // 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, int qpms_trans_calculator_get_AB_buf_p(const qpms_trans_calculator *c,
complex double *Adest, complex double *Bdest, 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, bool r_ge_d, qpms_bessel_t J,
complex double *bessel_buf, double *legendre_buf) { complex double *bessel_buf, double *legendre_buf) {
if (r_ge_d) J = QPMS_BESSEL_REGULAR; 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, int qpms_trans_calculator_get_AB_arrays_buf(const qpms_trans_calculator *c,
complex double *Adest, complex double *Bdest, complex double *Adest, complex double *Bdest,
size_t deststride, size_t srcstride, 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) { complex double *bessel_buf, double *legendre_buf) {
if (r_ge_d) J = QPMS_BESSEL_REGULAR; if (r_ge_d) J = QPMS_BESSEL_REGULAR;
if (0 == kdlj.r && 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, 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) { bool r_ge_d, qpms_bessel_t J) {
double leg[gsl_sf_legendre_array_n(n+nu)]; 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. 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, 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) { bool r_ge_d, qpms_bessel_t J) {
double leg[gsl_sf_legendre_array_n(n+nu+1)]; 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. 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, int qpms_trans_calculator_get_AB_p(const qpms_trans_calculator *c,
complex double *Adest, complex double *Bdest, 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) { bool r_ge_d, qpms_bessel_t J) {
double leg[gsl_sf_legendre_array_n(2*c->lMax+1)]; 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. 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, int qpms_trans_calculator_get_AB_arrays(const qpms_trans_calculator *c,
complex double *Adest, complex double *Bdest, complex double *Adest, complex double *Bdest,
size_t deststride, size_t srcstride, 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)]; 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. 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, 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, const qpms_vswf_set_spec_t *destspec, size_t deststride,
/// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm. /// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm.
const qpms_vswf_set_spec_t *srcspec, size_t srcstride, 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); TROPS_ONLY_AB_SYMMETRIC_NORMS_IMPLEMENTED(c->normalisation);
assert(c->normalisation == destspec->norm && c->normalisation == srcspec->norm); 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, const qpms_vswf_set_spec_t *destspec, size_t deststride,
/// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm. /// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm.
const qpms_vswf_set_spec_t *srcspec, size_t srcstride, 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 /// 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; kdlj.r *= k;
return qpms_trans_calculator_get_trans_array(c, target, return qpms_trans_calculator_get_trans_array(c, target,
destspec, deststride, srcspec, srcstride, kdlj, 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, complex double qpms_trans_calculator_get_A_ext(const qpms_trans_calculator *c,
int m, int n, int mu, int nu, 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) { 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); 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, complex double qpms_trans_calculator_get_B_ext(const qpms_trans_calculator *c,
int m, int n, int mu, int nu, 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) { 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); 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, int qpms_trans_calculator_get_AB_p_ext(const qpms_trans_calculator *c,
complex double *Adest, complex double *Bdest, complex double *Adest, complex double *Bdest,
int m, int n, int mu, int nu, 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) { 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); 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, int qpms_trans_calculator_get_AB_arrays_ext(const qpms_trans_calculator *c,
complex double *Adest, complex double *Bdest, complex double *Adest, complex double *Bdest,
size_t deststride, size_t srcstride, 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) { 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, return qpms_trans_calculator_get_AB_arrays(c,Adest,Bdest,deststride,srcstride,
kdlj, r_ge_d, J); kdlj, r_ge_d, J);
} }

24
qpms/translations.h
View File

@ -37,10 +37,10 @@
#endif #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); 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); bool r_ge_d, qpms_bessel_t J);
/// Structure holding the constant factors in normalisation operators. /// 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 *); void qpms_trans_calculator_free(qpms_trans_calculator *);
complex double qpms_trans_calculator_get_A(const qpms_trans_calculator *c, 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); bool r_ge_d, qpms_bessel_t J);
complex double qpms_trans_calculator_get_B(const qpms_trans_calculator *c, 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); bool r_ge_d, qpms_bessel_t J);
int qpms_trans_calculator_get_AB_p(const qpms_trans_calculator *c, int qpms_trans_calculator_get_AB_p(const qpms_trans_calculator *c,
complex double *Adest, complex double *Bdest, 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); bool r_ge_d, qpms_bessel_t J);
int qpms_trans_calculator_get_AB_arrays(const qpms_trans_calculator *c, int qpms_trans_calculator_get_AB_arrays(const qpms_trans_calculator *c,
complex double *Adest, complex double *Bdest, complex double *Adest, complex double *Bdest,
size_t deststride, size_t srcstride, 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 // TODO update the types later
complex double qpms_trans_calculator_get_A_ext(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, int m, int n, int mu, int nu, complex double kdlj_r,
double kdlj_th, double kdlj_phi, int r_ge_d, int J); 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, 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); 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, int qpms_trans_calculator_get_AB_p_ext(const qpms_trans_calculator *c,
complex double *Adest, complex double *Bdest, 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); 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, int qpms_trans_calculator_get_AB_arrays_ext(const qpms_trans_calculator *c,
complex double *Adest, complex double *Bdest, complex double *Adest, complex double *Bdest,
size_t deststride, size_t srcstride, 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); int r_ge_d, int J);
// Convenience functions using VSWF base specs // 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, const qpms_vswf_set_spec_t *destspec, size_t deststride,
/// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm. /// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm.
const qpms_vswf_set_spec_t *srcspec, size_t srcstride, 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 /// Version with \a k and cartesian particle positions
/// and with automatic \a r_ge_d = `false`. /// 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, const qpms_vswf_set_spec_t *destspec, size_t deststride,
/// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm. /// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm.
const qpms_vswf_set_spec_t *srcspec, size_t srcstride, 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 qpms_bessel_t J
/// Workspace has to be at least 2 * c->neleme**2 long /// Workspace has to be at least 2 * c->neleme**2 long
); );

9
qpms/vectors.h
View File

@ -333,6 +333,15 @@ static inline csph_t ccart2csph(const ccart3_t cart) {
return sph; 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 /// Coordinate transform of csph_t to ccart3_t
static inline ccart3_t csph2ccart(const csph_t sph) { static inline ccart3_t csph2ccart(const csph_t sph) {
ccart3_t cart; ccart3_t cart;