diff --git a/qpms/indexing.h b/qpms/indexing.h index 89398b2..5e97c5e 100644 --- a/qpms/indexing.h +++ b/qpms/indexing.h @@ -72,6 +72,18 @@ static inline qpms_errno_t qpms_uvswfi2tmn(qpms_uvswfi_t u, return QPMS_SUCCESS; } +/// Conversion from universal VSWF index u to type and the traditional layout index. +/** Does not for allow longitudinal waves. */ +static inline qpms_errno_t qpms_uvswfi2ty(qpms_uvswfi_t u, + qpms_vswf_type_t *t, qpms_y_t *y) { + *t = u & 3; + *y = u / 4 - 1; + if (*t == 0 || *t == 3) return QPMS_ERROR; + if (*y < 0) return QPMS_ERROR; + return QPMS_SUCCESS; +} + + /// Extract degree \a m from an universal VSWF index \a u. static inline qpms_m_t qpms_uvswfi2m(qpms_uvswfi_t u) { qpms_vswf_type_t t; qpms_m_t m; qpms_l_t n; diff --git a/qpms/translations.c b/qpms/translations.c index ee0f0bf..04b2da5 100644 --- a/qpms/translations.c +++ b/qpms/translations.c @@ -11,6 +11,7 @@ #include "kahansum.h" #include //abort() #include +#include "qpms_error.h" /* @@ -1153,7 +1154,59 @@ int qpms_trans_calculator_get_AB_arrays(const qpms_trans_calculator *c, bes, leg); } +// Convenience functions using VSWF base specs +qpms_errno_t qpms_trans_calculator_get_trans_array(const qpms_trans_calculator *c, + complex double *target, + /// Must be destspec->lMax <= c-> lMax && destspec->norm == c->norm. + 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) +{ + assert(c->normalisation == destspec->norm && c->normalisation == srcspec->norm); + assert(c->lMax >= destspec->lMax && c->lMax >= srcspec->lMax); + assert(!destspec->lMax_L && !srcspec->lMax_L); + complex double A[c->nelem][c->nelem]; + complex double B[c->nelem][c->nelem]; + qpms_errno_t retval = qpms_trans_calculator_get_AB_arrays(c, + A[0], B[0], c->nelem, 1, + kdlj, r_ge_d, J); + for (size_t desti = 0; desti < destspec->n; ++desti) { + qpms_y_t desty; qpms_vswf_type_t destt; + if(QPMS_SUCCESS != qpms_uvswfi2ty(destspec->ilist[desti], &destt, &desty)) + qpms_pr_error_at_flf(__FILE__,__LINE__,__func__, + "Invalid u. vswf index %llx.", destspec->ilist[desti]); + for (size_t srci = 0; srci < srcspec->n; ++srci){ + qpms_y_t srcy; qpms_vswf_type_t srct; + if(QPMS_SUCCESS != qpms_uvswfi2ty(srcspec->ilist[srci], &srct, &srcy)) + qpms_pr_error_at_flf(__FILE__,__LINE__,__func__, + "Invalid u. vswf index %llx.", srcspec->ilist[srci]); + target[srci * srcstride + desti * deststride] + = (srct == destt) ? A[desty][srcy] : B[desty][srcy]; + } + } + return QPMS_SUCCESS; +} +/// 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, + /// Must be destspec->lMax <= c-> lMax && destspec->norm == c->norm. + 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 + /// Workspace has to be at least 2 * c->neleme**2 long + ) +{ + sph_t kdlj = cart2sph(cart3_substract(destpos, srcpos)); + kdlj.r *= k; + return qpms_trans_calculator_get_trans_array(c, target, + destspec, deststride, srcspec, srcstride, kdlj, + false, QPMS_HANKEL_PLUS); +} #ifdef LATTICESUMS31 int qpms_trans_calculator_get_AB_arrays_e31z_both_points_and_shift(const qpms_trans_calculator *c, @@ -1270,7 +1323,7 @@ int qpms_trans_calculator_get_AB_arrays_e31z_both_points_and_shift(const qpms_tr } return 0; } -#endif LATTICESUMS_31 +#endif // LATTICESUMS_31 #ifdef LATTICESUMS32 diff --git a/qpms/translations.h b/qpms/translations.h index a80fde7..2368548 100644 --- a/qpms/translations.h +++ b/qpms/translations.h @@ -105,7 +105,6 @@ int qpms_trans_calculator_get_AB_arrays(const qpms_trans_calculator *c, size_t deststride, size_t srcstride, sph_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, @@ -126,6 +125,28 @@ int qpms_trans_calculator_get_AB_arrays_ext(const qpms_trans_calculator *c, double kdlj_r, double kdlj_theta, double kdlj_phi, int r_ge_d, int J); +// Convenience functions using VSWF base specs +qpms_errno_t qpms_trans_calculator_get_trans_array(const qpms_trans_calculator *c, + complex double *target, + /// Must be destspec->lMax <= c-> lMax && destspec->norm == c->norm. + 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); + +/// 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, + /// Must be destspec->lMax <= c-> lMax && destspec->norm == c->norm. + 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 + /// Workspace has to be at least 2 * c->neleme**2 long + ); + #ifdef LATTICESUMS_OLD // Short-range parts of the translation coefficients int qpms_trans_calculator_get_shortrange_AB_p(const qpms_trans_calculator *c, diff --git a/qpms/vswf.c b/qpms/vswf.c index 10af7ff..66cfb3b 100644 --- a/qpms/vswf.c +++ b/qpms/vswf.c @@ -82,6 +82,8 @@ qpms_vswf_set_spec_t *qpms_vswf_set_spec_from_lMax(qpms_l_t lMax, c->ilist[i++] = qpms_tmn2uvswfi(it ? QPMS_VSWF_MAGNETIC : QPMS_VSWF_ELECTRIC, m, n); c->norm = norm; + c->lMax = c->lMax_M = c->lMax_N = lMax; + c->lMax_L = 0; return c; }