diff --git a/qpms/qpms_error.h b/qpms/qpms_error.h index 1a72cf2..6d7e76c 100644 --- a/qpms/qpms_error.h +++ b/qpms/qpms_error.h @@ -55,7 +55,7 @@ qpms_dbgmsg_flags qpms_dbgmsg_enable(qpms_dbgmsg_flags types); #define QPMS_UNTESTED {\ static bool already_bitched = false; \ if (!already_bitched) {\ - qpms_warn_at_flf(__FILE__,__LINE__,__func__,"Warning: untested function/feature!")\ + qpms_warn_at_flf(__FILE__,__LINE__,__func__,"Warning: untested function/feature!");\ already_bitched = true;\ }\ } diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index e5107b9..e11b524 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -15,6 +15,7 @@ #include "translations.h" #include "tmatrices.h" #include +#include "kahansum.h" #ifdef QPMS_SCATSYSTEM_USE_OWN_BLAS #include "qpmsblas.h" @@ -1460,23 +1461,34 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR( return target_packed; } -#if 0 -ccart3_t qpms_scatsys_eval_field(const qpms_scatsys_t *ss, +#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, - complex double omega, complex double refindex) { + const complex double k) { QPMS_UNTESTED; - ccart3_t result = {0,0,0}; - ccart3_t result_kahanc = {0,0,0}; + ccart3_t res = {0,0,0}; + ccart3_t res_kc = {0,0,0}; // kahan sum compensation for (qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi) { - const cart3_t particle_pos = ss->p[pi]; - const csph_t kr = sph_cscale(...); - ...; + const qpms_vswf_set_spec_t *bspec = ss->tm[ss->p[pi].tmatrix_id]->spec; + const cart3_t particle_pos = ss->p[pi].pos; + const complex double *particle_cv = cvf + ss->fecv_pstarts[pi]; + + const csph_t kr = sph_cscale(k, cart2sph( + cart3_substract(where, particle_pos))); + const csphvec_t E_sph = qpms_eval_uvswf(bspec, particle_cv, kr, + QPMS_HANKEL_PLUS); + const ccart3_t E_cart = csphvec2ccart_csph(E_sph, kr); + ckahanadd(&(res.x), &(res_kc.x), E_cart.x); + ckahanadd(&(res.y), &(res_kc.y), E_cart.y); + ckahanadd(&(res.z), &(res_kc.z), E_cart.z); } + return res; } +#endif - -ccart3_t qpms_scatsys_eval_field_irrep(const qpms_scatsys_t *ss, +#if 0 +ccart3_t qpms_scatsys_eval_E_irrep(const qpms_scatsys_t *ss, qpms_iri_t iri, const complex double *cvr, cart3_t where) { TODO; } diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index 85f56bf..9cca6f1 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -319,27 +319,31 @@ complex double *qpms_orbit_irrep_basis( /// The index of the irreducible representation of sym. const qpms_iri_t iri); -/** Evaluates scattered fields (corresponding to a given excitation vector) - * at a given point. +/// Evaluates scattered fields (corresponding to a given excitation vector) at a given point. +/** + * By scattered field, one assumes a linear combination of positive-Hankel-type + * spherical waves. * * \return Complex electric field at the point defined by \a where. */ -ccart3_t qpms_scatsys_eval_field(const qpms_scatsys_t *ss, +ccart3_t qpms_scatsys_eval_E(const qpms_scatsys_t *ss, const complex double *coeff_vector, ///< A full-length excitation vector. cart3_t where, ///< Evaluation point. complex double k ///< Wave number. ); +#if 0 /** Evaluates partial scattered fields (corresponding to a given irrep-reduced excitation vector) * at a given point. * * \return Complex electric field at the point defined by \a where. */ -ccart3_t qpms_scatsys_eval_field_irrep(const qpms_scatsys_t *ss, +ccart3_t qpms_scatsys_eval_E_irrep(const qpms_scatsys_t *ss, qpms_iri_t iri, ///< Irreducible representation const complex double *coeff_vector, ///< A reduced excitation vector, corresponding to \a iri. cart3_t where, ///< Evaluation point. complex double k ///< Wave number. ); +#endif #endif //QPMS_SCATSYSTEM_H