From 37b864effcc164f295f54f706ee142075db45d40 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Mon, 24 Jun 2019 17:19:25 +0300 Subject: [PATCH] uvswf evaluation with complex kr; ss field ewald done, untested Former-commit-id: 63048ddf2a9260742510a4fa246a59b8360381ce --- qpms/qpms_error.h | 2 +- qpms/scatsystem.c | 2 -- qpms/vswf.c | 46 +++++++++++++++++++++++++++++++--------------- qpms/vswf.h | 33 ++++++++++++++++++++++++++++----- 4 files changed, 60 insertions(+), 23 deletions(-) diff --git a/qpms/qpms_error.h b/qpms/qpms_error.h index 6d7e76c..154cf94 100644 --- a/qpms/qpms_error.h +++ b/qpms/qpms_error.h @@ -66,7 +66,7 @@ qpms_dbgmsg_flags qpms_dbgmsg_enable(qpms_dbgmsg_flags types); #define QPMS_ENSURE(x, msg, ...) {if(!(x)) qpms_pr_error_at_flf(__FILE__,__LINE__,__func__,msg, ##__VA_ARGS__); } -#define QPMS_ASSERT(x) {if(!(x)) qpms_pr_error_at_flf(__FILE__,__LINE__,__func__, "Unexpected error. This is most certainly a bug.") +#define QPMS_ASSERT(x) {if(!(x)) qpms_pr_error_at_flf(__FILE__,__LINE__,__func__, "Unexpected error. This is most certainly a bug.");} #define QPMS_NOT_IMPLEMENTED(msg, ...) qpms_pr_error_at_flf(__FILE__,__LINE__,__func__, \ "Not implemented:" msg, ##__VA_ARGS__) diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index e11b524..65a166d 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -1461,7 +1461,6 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR( return target_packed; } -#if 0 // must allow for complex kr in qpms_eval_uvswf ccart3_t qpms_scatsys_eval_E(const qpms_scatsys_t *ss, const complex double *cvf, const cart3_t where, const complex double k) { @@ -1485,7 +1484,6 @@ ccart3_t qpms_scatsys_eval_E(const qpms_scatsys_t *ss, } return res; } -#endif #if 0 ccart3_t qpms_scatsys_eval_E_irrep(const qpms_scatsys_t *ss, diff --git a/qpms/vswf.c b/qpms/vswf.c index 574fde5..86afc0a 100644 --- a/qpms/vswf.c +++ b/qpms/vswf.c @@ -7,6 +7,7 @@ #include "qpms_specfunc.h" #include #include +#include "qpms_error.h" qpms_vswf_set_spec_t *qpms_vswf_set_spec_init() { @@ -150,18 +151,19 @@ void qpms_vswfset_sph_pfree(qpms_vswfset_sph_t *w) { free(w); } -csphvec_t qpms_vswf_L00(sph_t kdrj, qpms_bessel_t btyp, +csphvec_t qpms_vswf_L00(csph_t kr, qpms_bessel_t btyp, qpms_normalisation_t norm) { QPMS_UNTESTED; // CHECKME Is it OK to ignore norm?? (Is L_0^0 the same for all conventions?) complex double bessel0; QPMS_ENSURE_SUCCESS(qpms_sph_bessel_fill(btyp, 0, kr.r, &bessel0)); - return {0.25 * M_2_SQRTPI * bessel0, 0, 0}; + csphvec_t result = {0.25 * M_2_SQRTPI * bessel0, 0, 0}; + return result; } -qpms_errno_t qpms_vswf_fill(csphvec_t *const longtarget, +qpms_errno_t qpms_vswf_fill_csph(csphvec_t *const longtarget, csphvec_t * const mgtarget, csphvec_t * const eltarget, qpms_l_t lMax, - sph_t kr, qpms_bessel_t btyp, qpms_normalisation_t norm) { + csph_t kr, qpms_bessel_t btyp, qpms_normalisation_t norm) { assert(lMax >= 1); complex double *bessel = malloc((lMax+1)*sizeof(complex double)); if(qpms_sph_bessel_fill(btyp, lMax, kr.r, bessel)) abort(); @@ -214,6 +216,14 @@ qpms_errno_t qpms_vswf_fill(csphvec_t *const longtarget, return QPMS_SUCCESS; } +qpms_errno_t qpms_vswf_fill(csphvec_t *const longtarget, + csphvec_t * const mgtarget, csphvec_t * const eltarget, qpms_l_t lMax, + sph_t kr, qpms_bessel_t btyp, qpms_normalisation_t norm) { + csph_t krc = {kr.r, kr.theta, kr.phi}; + return qpms_vswf_fill_csph(longtarget, mgtarget, eltarget, lMax, + krc, btyp, norm); +} + // consistency check: this should give the same results as the above function (up to rounding errors) qpms_errno_t qpms_vswf_fill_alternative(csphvec_t *const longtarget, csphvec_t * const mgtarget, csphvec_t * const eltarget, qpms_l_t lMax, sph_t kr, @@ -385,7 +395,7 @@ qpms_errno_t qpms_planewave2vswf_fill_cart(cart3_t wavedir_cart /*allow complex longcoeff, mgcoeff, elcoeff, lMax, norm); } -csphvec_t qpms_eval_vswf(sph_t kr, +csphvec_t qpms_eval_vswf_csph(csph_t kr, complex double * const lc, complex double *const mc, complex double *const ec, qpms_l_t lMax, qpms_bessel_t btyp, qpms_normalisation_t norm) { @@ -398,7 +408,7 @@ csphvec_t qpms_eval_vswf(sph_t kr, if(lc) lset = malloc(nelem * sizeof(csphvec_t)); if(mc) mset = malloc(nelem * sizeof(csphvec_t)); if(ec) eset = malloc(nelem * sizeof(csphvec_t)); - qpms_vswf_fill(lset, mset, eset, lMax, kr, btyp, norm); + qpms_vswf_fill_csph(lset, mset, eset, lMax, kr, btyp, norm); if(lc) for(qpms_y_t y = 0; y < nelem; ++y) csphvec_kahanadd(&lsum, &lcomp, csphvec_scale(lc[y], lset[y])); if(mc) for(qpms_y_t y = 0; y < nelem; ++y) @@ -414,18 +424,24 @@ csphvec_t qpms_eval_vswf(sph_t kr, return esum; } +csphvec_t qpms_eval_vswf(sph_t kr, + complex double * const lc, complex double *const mc, complex double *const ec, + qpms_l_t lMax, qpms_bessel_t btyp, qpms_normalisation_t norm) { + csph_t krc = {kr.r, kr.theta, kr.phi}; + return qpms_eval_vswf_csph(krc, lc, mc, ec, lMax, btyp, norm); +} -csphvec_t qpms_eval_uvswf(const qpms_vswf_set_spec_t *setspec, - const complex double *coeffs, const sph_t kr, +csphvec_t qpms_eval_uvswf(const qpms_vswf_set_spec_t *bspec, + const complex double *coeffs, const csph_t kr, const qpms_bessel_t btyp) { QPMS_UNTESTED; - const qpms_l_t lMax = b->lMax; + const qpms_l_t lMax = bspec->lMax; complex double *cM = NULL, *cN = NULL, *cL = NULL, cL00 = 0; - if (b->lMax_L > 0) + if (bspec->lMax_L > 0) QPMS_CRASHING_CALLOC(cL, lMax, sizeof(complex double)); - if (b->lMax_M > 0) + if (bspec->lMax_M > 0) QPMS_CRASHING_CALLOC(cM, lMax, sizeof(complex double)); - if (b->lMax_N > 0) + if (bspec->lMax_N > 0) QPMS_CRASHING_CALLOC(cN, lMax, sizeof(complex double)); for (size_t i = 0; i < bspec->n; ++i) { if (bspec->ilist[i] == 0) // L00, needs special care @@ -452,10 +468,10 @@ csphvec_t qpms_eval_uvswf(const qpms_vswf_set_spec_t *setspec, } } } - csphvec_t result = qpms_eval_vswf(kr, cL, cM, cN, lMax, btyp, norm); + csphvec_t result = qpms_eval_vswf_csph(kr, cL, cM, cN, lMax, btyp, bspec->norm); + free(cM); free(cN); free(cL); if(cL00) result = csphvec_add(result, - scphvec_scale(cL00, qpms_vswf_L00(kr, btyp, norm))); + csphvec_scale(cL00, qpms_vswf_L00(kr, btyp, bspec->norm))); return result; } - diff --git a/qpms/vswf.h b/qpms/vswf.h index 776bc32..ca5da8c 100644 --- a/qpms/vswf.h +++ b/qpms/vswf.h @@ -42,20 +42,22 @@ static inline ssize_t qpms_vswf_set_spec_find_uvswfi(const qpms_vswf_set_spec_t return -1; } -/// NOT IMPLEMENTED Evaluates a set of VSWF basis functions at a given point. +/// Evaluates a set of VSWF basis functions at a given point. /** The list of basis wave indices is specified in \a setspec; * \a setspec->norm must be set as well. */ qpms_errno_t qpms_uvswf_fill( - csphvec_t *const target, + csphvec_t *const target, ///< Target array of size at least setspec->n. const qpms_vswf_set_spec_t *setspec, - sph_t evaluation_point, qpms_bessel_t btyp); + csph_t kr, ///< Evaluation point. + qpms_bessel_t btyp); /// Evaluates field specified by SVWF coefficients at a given point. /** SVWF coefficients in \a coeffs must be ordered according to \a setspec->ilist. */ csphvec_t qpms_eval_uvswf(const qpms_vswf_set_spec_t *setspec, - const complex double *coeffs, sph_t evaluation_point, + const complex double *coeffs, ///< SVWF coefficient vector of size setspec->n. + csph_t kr, ///< Evaluation point. qpms_bessel_t btyp); /// Electric wave N. @@ -86,7 +88,7 @@ qpms_errno_t qpms_legendre_deriv_y_fill(double *where, double *where_deriv, doub /// Evaluate the zeroth-degree longitudinal VSWF \f$ \mathbf{L}_0^0 \f$. csphvec_t qpms_vswf_L00( - sph_t kdrj, //< VSWF evaluation point. + csph_t kdrj, //< VSWF evaluation point. qpms_bessel_t btyp, qpms_normalisation_t norm); @@ -106,10 +108,28 @@ qpms_errno_t qpms_vswf_fill( qpms_l_t lMax, //< Maximum multipole degree to be calculated. sph_t kdrj, //< VSWF evaluation point. qpms_bessel_t btyp, qpms_normalisation_t norm); + // Should give the same results: for consistency checks qpms_errno_t qpms_vswf_fill_alternative(csphvec_t *resultL, csphvec_t *resultM, csphvec_t *resultN, qpms_l_t lMax, sph_t kdrj, qpms_bessel_t btyp, qpms_normalisation_t norm); +/// Evaluate VSWFs at a given point from \a l = 1 up to a given degree \a lMax (complex \a kr version). +/** + * The target arrays \a resultL, \a resultM, \a resultN have to be large enough to contain + * \a lMax * (\a lMax + 2) elements. If NULL is passed instead, the corresponding SVWF type + * is not evaluated. + * + * Does not evaluate the zeroth-order wave \f$ \mathbf{L}_0^0 \f$. + * If you need that, use qpms_vswf_L00(). + */ +qpms_errno_t qpms_vswf_fill_csph( + csphvec_t *resultL, //< Target array for longitudinal VSWFs. + csphvec_t *resultM, //< Target array for magnetic VSWFs. + csphvec_t *resultN, //< Target array for electric VSWFs. + qpms_l_t lMax, //< Maximum multipole degree to be calculated. + csph_t kdrj, //< VSWF evaluation point. + qpms_bessel_t btyp, qpms_normalisation_t norm); + qpms_errno_t qpms_vecspharm_fill(csphvec_t *const a1target, csphvec_t *const a2target, csphvec_t *const a3target, qpms_l_t lMax, sph_t dir, qpms_normalisation_t norm); qpms_errno_t qpms_vecspharm_dual_fill(csphvec_t *const a1target, csphvec_t *const a2target, csphvec_t *const a3target, @@ -127,6 +147,9 @@ csphvec_t qpms_eval_vswf(sph_t where, complex double *longcoeffs, complex double *mgcoeffs, complex double *elcoeffs, qpms_l_t lMax, qpms_bessel_t btyp, qpms_normalisation_t norm); +csphvec_t qpms_eval_vswf_csph(csph_t where, + complex double *longcoeffs, complex double *mgcoeffs, complex double *elcoeffs, + qpms_l_t lMax, qpms_bessel_t btyp, qpms_normalisation_t norm); qpms_vswfset_sph_t *qpms_vswfset_make(qpms_l_t lMax, sph_t kdlj, qpms_bessel_t btyp, qpms_normalisation_t norm);//NI