diff --git a/qpms/groups.h b/qpms/groups.h index 34f2c44..6b05332 100644 --- a/qpms/groups.h +++ b/qpms/groups.h @@ -25,6 +25,7 @@ #define QPMS_GROUPS_H #include "qpms_types.h" +#include /// To be used only in qpms_finite_group_t struct qpms_finite_group_irrep_t { diff --git a/qpms/normalisation.h b/qpms/normalisation.h index c7ee22b..5ab7099 100644 --- a/qpms/normalisation.h +++ b/qpms/normalisation.h @@ -207,4 +207,77 @@ static inline complex double qpms_spharm_azimuthal_part_derivative_div_m(qpms_no QPMS_WTF; } } + +#if 0 // legacy code moved from qpms_types.h. TODO cleanup +/// Returns the normalisation convention code without the Condon-Shortley phase. +static inline int qpms_normalisation_t_normonly(qpms_normalisation_t norm) { + return norm & (~QPMS_NORMALISATION_T_CSBIT); +} + +/* Normalisation of the spherical waves is now scattered in at least three different files: + * here, we have the norm in terms of radiated power of outgoing wave. + * In file legendre.c, function qpms_pitau_get determines the norm used in the vswf.c + * spherical vector wave norms. The "dual" waves in vswf.c use the ..._abssquare function below. + * In file translations.c, the normalisations are again set by hand using the normfac and lognormfac + * functions. + */ +#include +#include +// relative to QPMS_NORMALISATION_KRISTENSSON_CS, i.e. +// P_l^m[normtype] = P_l^m[Kristensson] +static inline double qpms_normalisation_t_factor(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) { + int csphase = qpms_normalisation_t_csphase(norm); + norm = qpms_normalisation_t_normonly(norm); + double factor; + switch (norm) { + case QPMS_NORMALISATION_KRISTENSSON: + factor = 1.; + break; + case QPMS_NORMALISATION_TAYLOR: + factor = sqrt(l*(l+1)); + break; + case QPMS_NORMALISATION_NONE: + factor = sqrt(l*(l+1) * 4 * M_PI / (2*l+1) * exp(lgamma(l+m+1)-lgamma(l-m+1))); + break; +#ifdef USE_XU_ANTINORMALISATION // broken probably in legendre.c + case QPMS_NORMALISATION_XU: + factor = sqrt(4 * M_PI) / (2*l+1) * exp(lgamma(l+m+1)-lgamma(l-m+1)); + break; +#endif + default: + assert(0); + } + factor *= (m%2)?(-csphase):1; + return factor; +} + + +// TODO move elsewhere +static inline double qpms_normalisation_t_factor_abssquare(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) { + norm = qpms_normalisation_t_normonly(norm); + switch (norm) { + case QPMS_NORMALISATION_KRISTENSSON: + return 1.; + break; + case QPMS_NORMALISATION_TAYLOR: + return l*(l+1); + break; + case QPMS_NORMALISATION_NONE: + return l*(l+1) * 4 * M_PI / (2*l+1) * exp(lgamma(l+m+1)-lgamma(l-m+1)); + break; +#ifdef USE_XU_ANTINORMALISATION // broken probably in legendre.c + case QPMS_NORMALISATION_XU: + { + double fac = sqrt(4 * M_PI) / (2*l+1) * exp(lgamma(l+m+1)-lgamma(l-m+1)); + return fac * fac; + } + break; +#endif + default: + assert(0); + return NAN; + } +} +#endif + #endif //NORMALISATION_H diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx index 3dc266a..f158954 100644 --- a/qpms/qpms_c.pyx +++ b/qpms/qpms_c.pyx @@ -25,19 +25,15 @@ class VSWFType(enum.IntEnum): L = QPMS_VSWF_LONGITUDINAL class VSWFNorm(enum.IntEnum): - #XU = QPMS_NORMALISATION_XU - #XU_CS = QPMS_NORMALISATION_XU_CS - NONE = QPMS_NORMALISATION_NONE - NONE_CS = QPMS_NORMALISATION_NONE_CS - POWER = QPMS_NORMALISATION_POWER - POWER_CS = QPMS_NORMALISATION_POWER_CS - SPHARM = QPMS_NORMALISATION_SPHARM - SPHARM_CS = QPMS_NORMALISATION_SPHARM_CS + # TODO try to make this an enum.IntFlag if supported + # TODO add the other flags from qpms_normalisation_t as well + UNNORM = QPMS_NORMALISATION_NORM_NONE + UNNORM_CS = QPMS_NORMALISATION_NORM_NONE | QPMS_NORMALISATION_CSPHASE + POWERNORM = QPMS_NORMALISATION_NORM_POWER + POWERNORM_CS = QPMS_NORMALISATION_NORM_POWER | QPMS_NORMALISATION_CSPHASE + SPHARMNORM = QPMS_NORMALISATION_NORM_SPHARM + SPHARMNORM_CS = QPMS_NORMALISATION_NORM_SPHARM | QPMS_NORMALISATION_CSPHASE UNDEF = QPMS_NORMALISATION_UNDEF - KRISTENSSON = QPMS_NORMALISATION_KRISTENSSON - KRISTENSSON_CS = QPMS_NORMALISATION_KRISTENSSON_CS - TAYLOR = QPMS_NORMALISATION_TAYLOR - TAYLOR_CS = QPMS_NORMALISATION_TAYLOR_CS try: class DebugFlags(enum.IntFlag): # Should be IntFlag if python version >= 3.6 @@ -726,7 +722,7 @@ cdef class BaseSpec: if 'norm' in kwargs.keys(): self.s.norm = kwargs['norm'] else: - self.s.norm = QPMS_NORMALISATION_POWER_CS + self.s.norm = QPMS_NORMALISATION_NORM_POWER | QPMS_NORMALISATION_CSPHASE # set the other metadata cdef qpms_l_t l self.s.lMax_L = -1 diff --git a/qpms/qpms_types.h b/qpms/qpms_types.h index 0bf4263..e72ef0f 100644 --- a/qpms/qpms_types.h +++ b/qpms/qpms_types.h @@ -166,80 +166,6 @@ static inline int qpms_normalisation_t_csphase(qpms_normalisation_t norm) { return (norm & QPMS_NORMALISATION_CSPHASE)? -1 : 1; } -#if 0 -/// Returns the normalisation convention code without the Condon-Shortley phase. -static inline int qpms_normalisation_t_normonly(qpms_normalisation_t norm) { - return norm & (~QPMS_NORMALISATION_T_CSBIT); -} - -// TODO move the inlines elsewhere -/* Normalisation of the spherical waves is now scattered in at least three different files: - * here, we have the norm in terms of radiated power of outgoing wave. - * In file legendre.c, function qpms_pitau_get determines the norm used in the vswf.c - * spherical vector wave norms. The "dual" waves in vswf.c use the ..._abssquare function below. - * In file translations.c, the normalisations are again set by hand using the normfac and lognormfac - * functions. - */ -#include -#include -// relative to QPMS_NORMALISATION_KRISTENSSON_CS, i.e. -// P_l^m[normtype] = P_l^m[Kristensson] -static inline double qpms_normalisation_t_factor(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) { - int csphase = qpms_normalisation_t_csphase(norm); - norm = qpms_normalisation_t_normonly(norm); - double factor; - switch (norm) { - case QPMS_NORMALISATION_KRISTENSSON: - factor = 1.; - break; - case QPMS_NORMALISATION_TAYLOR: - factor = sqrt(l*(l+1)); - break; - case QPMS_NORMALISATION_NONE: - factor = sqrt(l*(l+1) * 4 * M_PI / (2*l+1) * exp(lgamma(l+m+1)-lgamma(l-m+1))); - break; -#ifdef USE_XU_ANTINORMALISATION // broken probably in legendre.c - case QPMS_NORMALISATION_XU: - factor = sqrt(4 * M_PI) / (2*l+1) * exp(lgamma(l+m+1)-lgamma(l-m+1)); - break; -#endif - default: - assert(0); - } - factor *= (m%2)?(-csphase):1; - return factor; -} - - -// TODO move elsewhere -static inline double qpms_normalisation_t_factor_abssquare(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) { - norm = qpms_normalisation_t_normonly(norm); - switch (norm) { - case QPMS_NORMALISATION_KRISTENSSON: - return 1.; - break; - case QPMS_NORMALISATION_TAYLOR: - return l*(l+1); - break; - case QPMS_NORMALISATION_NONE: - return l*(l+1) * 4 * M_PI / (2*l+1) * exp(lgamma(l+m+1)-lgamma(l-m+1)); - break; -#ifdef USE_XU_ANTINORMALISATION // broken probably in legendre.c - case QPMS_NORMALISATION_XU: - { - double fac = sqrt(4 * M_PI) / (2*l+1) * exp(lgamma(l+m+1)-lgamma(l-m+1)); - return fac * fac; - } - break; -#endif - default: - assert(0); - return NAN; - } -} -#endif - - /// Bessel function kinds. typedef enum { QPMS_BESSEL_REGULAR = 1, ///< regular (spherical) Bessel function \a j (Bessel function of the first kind) diff --git a/qpms/symmetries.c b/qpms/symmetries.c index 3ee2826..342f80b 100644 --- a/qpms/symmetries.c +++ b/qpms/symmetries.c @@ -2,38 +2,46 @@ #include "tiny_inlines.h" #include "indexing.h" #include "wigner.h" +#include "qpms_error.h" // TODO at some point, maybe support also other norms. // (perhaps just use qpms_normalisation_t_factor() at the right places) static inline void check_norm_compat(const qpms_vswf_set_spec_t *s) { - switch (qpms_normalisation_t_normonly(s->norm)) { - case QPMS_NORMALISATION_POWER: + switch (s->norm & QPMS_NORMALISATION_NORM_BITS) { + case QPMS_NORMALISATION_NORM_POWER: break; - case QPMS_NORMALISATION_SPHARM: + case QPMS_NORMALISATION_NORM_SPHARM: break; default: - abort(); // Only SPHARM and POWER norms are supported right now. + QPMS_NOT_IMPLEMENTED("At the moment, only spherical harmonics of spherical harmonics or power normalisations implemented."); } } + +static inline void ONLY_EIMF_IMPLEMENTED(const qpms_normalisation_t norm) +{ + if (norm & QPMS_NORMALISATION_SPHARM_REAL) + QPMS_NOT_IMPLEMENTED("Support for real spherical harmonics not implemented yet."); +} + + // Used in the functions below to ensure memory allocation and checks for bspec validity static inline complex double *ensure_alloc(complex double *target, const qpms_vswf_set_spec_t *bspec) { check_norm_compat(bspec); const size_t n = bspec->n; if (target == NULL) - target = malloc(n * n * sizeof(complex double)); - if (target == NULL) abort(); + QPMS_CRASHING_MALLOC(target, n * n * sizeof(complex double)); return target; } - complex double *qpms_zflip_uvswi_dense( complex double *target, const qpms_vswf_set_spec_t *bspec) { + check_norm_compat(bspec); target = ensure_alloc(target, bspec); const size_t n = bspec->n; @@ -69,6 +77,8 @@ complex double *qpms_yflip_uvswi_dense( complex double *target, const qpms_vswf_set_spec_t *bspec) { + check_norm_compat(bspec); + ONLY_EIMF_IMPLEMENTED(bspec->norm); target = ensure_alloc(target, bspec); const size_t n = bspec->n; @@ -104,6 +114,8 @@ complex double *qpms_xflip_uvswi_dense( complex double *target, const qpms_vswf_set_spec_t *bspec) { + check_norm_compat(bspec); + ONLY_EIMF_IMPLEMENTED(bspec->norm); target = ensure_alloc(target, bspec); const size_t n = bspec->n; @@ -142,6 +154,9 @@ complex double *qpms_zrot_uvswi_dense( double phi ///< Rotation angle ) { + QPMS_UNTESTED; // not sure about the C.-S. phase. Don't forget documenting it as well. + check_norm_compat(bspec); + ONLY_EIMF_IMPLEMENTED(bspec->norm); target = ensure_alloc(target, bspec); const size_t n = bspec->n; @@ -182,6 +197,9 @@ complex double *qpms_irot3_uvswfi_dense( const qpms_vswf_set_spec_t *bspec, const qpms_irot3_t t) { + QPMS_UNTESTED; // not sure about the C.-S. phase. Don't forget documenting it as well. + check_norm_compat(bspec); + ONLY_EIMF_IMPLEMENTED(bspec->norm); target = ensure_alloc(target, bspec); const size_t n = bspec->n;