From 45b82338622fec1d9ee7eeec25993b13659c6403 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Sun, 14 Jul 2019 09:33:30 +0300 Subject: [PATCH 1/3] Tag symmetries.c with non-implemented macros etc. Former-commit-id: 0cf1dc127b43bdded6471501a666ce325e5528f5 --- qpms/symmetries.c | 32 +++++++++++++++++++++++++------- 1 file changed, 25 insertions(+), 7 deletions(-) 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; From 0610f988542ba17fc9a2e6e71f78efcb41313115 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Sun, 14 Jul 2019 09:48:29 +0300 Subject: [PATCH 2/3] Move some legacy code from qpms_types.h to normalisation.h; to be pruned later. Former-commit-id: 0078e465757bc40f574837f196a7ab0c08529db8 --- qpms/groups.h | 1 + qpms/normalisation.h | 73 +++++++++++++++++++++++++++++++++++++++++++ qpms/qpms_types.h | 74 -------------------------------------------- 3 files changed, 74 insertions(+), 74 deletions(-) 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_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) From 005bdaa2be6bbc327fea4b4f2ab3f658d65c8da0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Sun, 14 Jul 2019 10:02:37 +0300 Subject: [PATCH 3/3] Use the new qpms_normalisation_t also in cython parts Former-commit-id: 58eb1c22724d3a94abb88d4c7b9af668a2faf116 --- qpms/qpms_c.pyx | 22 +++++++++------------- qpms/qpms_cdefs.pxd | 30 ++++++++++++++++++------------ 2 files changed, 27 insertions(+), 25 deletions(-) diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx index 81f2a67..f560d52 100644 --- a/qpms/qpms_c.pyx +++ b/qpms/qpms_c.pyx @@ -20,19 +20,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 @@ -721,7 +717,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_cdefs.pxd b/qpms/qpms_cdefs.pxd index 9c61760..5df5512 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -25,18 +25,24 @@ cdef extern from "qpms_types.h": cart2_t cart2 pol_t pol ctypedef enum qpms_normalisation_t: - QPMS_NORMALISATION_XU - QPMS_NORMALISATION_XU_CS - QPMS_NORMALISATION_NONE - QPMS_NORMALISATION_NONE_CS - QPMS_NORMALISATION_KRISTENSSON - QPMS_NORMALISATION_KRISTENSSON_CS - QPMS_NORMALISATION_POWER - QPMS_NORMALISATION_POWER_CS - QPMS_NORMALISATION_TAYLOR - QPMS_NORMALISATION_TAYLOR_CS - QPMS_NORMALISATION_SPHARM - QPMS_NORMALISATION_SPHARM_CS + QPMS_NORMALISATION_CONVENTION_KRISTENSSON + QPMS_NORMALISATION_CONVENTION_KRISTENSSON_REAL + QPMS_NORMALISATION_CONVENTION_SCUFF + QPMS_NORMALISATION_CSPHASE + QPMS_NORMALISATION_INVERSE + QPMS_NORMALISATION_L_I + QPMS_NORMALISATION_L_MINUS + QPMS_NORMALISATION_M_I + QPMS_NORMALISATION_M_MINUS + QPMS_NORMALISATION_NORM_BITS + QPMS_NORMALISATION_NORM_BITSTART + QPMS_NORMALISATION_NORM_NONE + QPMS_NORMALISATION_NORM_POWER + QPMS_NORMALISATION_NORM_SPHARM + QPMS_NORMALISATION_N_I + QPMS_NORMALISATION_N_MINUS + QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE + QPMS_NORMALISATION_SPHARM_REAL QPMS_NORMALISATION_UNDEF ctypedef enum qpms_bessel_t: QPMS_BESSEL_REGULAR