Merge branch 'newnorm' of necada:~/repo/qpms into newnorm

Former-commit-id: 19b6714173e1e205b0d9c5fb6b853b09f892bae0

qpms/groups.h

@@ -25,6 +25,7 @@
 #define QPMS_GROUPS_H
 
 #include "qpms_types.h"
+#include <assert.h>
 
 /// To be used only in qpms_finite_group_t
 struct qpms_finite_group_irrep_t {

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 <math.h>
+#include <assert.h>
+// 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

qpms/qpms_c.pyx

@@ -25,19 +25,15 @@ class VSWFType(enum.IntEnum):
     L = QPMS_VSWF_LONGITUDINAL
 
 class VSWFNorm(enum.IntEnum):
-    #XU = QPMS_NORMALISATION_XU
+    # TODO try to make this an enum.IntFlag if supported
-    #XU_CS = QPMS_NORMALISATION_XU_CS
+    # TODO add the other flags from qpms_normalisation_t as well
-    NONE = QPMS_NORMALISATION_NONE
+    UNNORM = QPMS_NORMALISATION_NORM_NONE
-    NONE_CS = QPMS_NORMALISATION_NONE_CS
+    UNNORM_CS = QPMS_NORMALISATION_NORM_NONE | QPMS_NORMALISATION_CSPHASE
-    POWER = QPMS_NORMALISATION_POWER
+    POWERNORM = QPMS_NORMALISATION_NORM_POWER
-    POWER_CS = QPMS_NORMALISATION_POWER_CS
+    POWERNORM_CS = QPMS_NORMALISATION_NORM_POWER | QPMS_NORMALISATION_CSPHASE
-    SPHARM = QPMS_NORMALISATION_SPHARM
+    SPHARMNORM = QPMS_NORMALISATION_NORM_SPHARM
-    SPHARM_CS = QPMS_NORMALISATION_SPHARM_CS
+    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

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 <math.h>
-#include <assert.h>
-// 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)

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)) {
+  switch (s->norm & QPMS_NORMALISATION_NORM_BITS) {
-    case QPMS_NORMALISATION_POWER:
+    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));
+    QPMS_CRASHING_MALLOC(target, n * n * sizeof(complex double));
-  if (target == NULL) abort();
   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;