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

uvswf evaluation with complex kr; ss field ewald done, untested 

Former-commit-id: 63048ddf2a9260742510a4fa246a59b8360381ce
This commit is contained in:
Marek Nečada 2019-06-24 17:19:25 +03:00
parent a3f5f98736
commit 37b864effc
4 changed files with 60 additions and 23 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
qpms/qpms_error.h
View File

@ -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__)

2
qpms/scatsystem.c
View File

@ -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,

46
qpms/vswf.c
View File

@ -7,6 +7,7 @@
#include "qpms_specfunc.h"
#include <stdlib.h>
#include <string.h>
#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;
}

33
qpms/vswf.h
View File

@ -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