New struct qpms_vswf_set_spec_t for specifying VSWF sets.

Former-commit-id: 56b44a931e68b57a2c786e69e859b3e0514b8a08

qpms/indexing.h

@@ -58,15 +58,16 @@ static inline qpms_uvswfi_t qpms_tmn2uvswfi(
 }
 
 /// Conversion from universal VSWF index u to type, order and degree.
-/** Crashes (abort()) the program if the u value is invalid. */
-static inline void qpms_uvswfi2tmn(qpms_uvswfi_t u,
+/** Returns a non-zero value if the u value is invalid. */
+static inline qpms_errno_t qpms_uvswfi2tmn(qpms_uvswfi_t u,
 		qpms_vswf_type_t *t, qpms_m_t *m, qpms_l_t *n) {
 	*t = u & 3;
 	qpms_y_sc_t y_sc = u / 4;
 	// Test validity
-	if (*t == 3) abort(); // VSWF type code invalid
-	if (*t && !y_sc) abort(); // l == 0 for transversal wave
+	if (*t == 3) return QPMS_ERROR; // VSWF type code invalid, TODO WARN
+	if (*t && !y_sc) return QPMS_ERROR; // l == 0 for transversal wave, TODO WARN
 	qpms_y2mn_sc_p(y_sc, m, n);
+	return QPMS_SUCCESS;
 }
 
 #endif //QPMS_INDEXING_H

qpms/vswf.c

@@ -8,6 +8,52 @@
 #include <stdlib.h>
 #include <string.h>
 
+
+qpms_vswf_set_spec_t *qpms_vswf_set_spec_init() {
+  qpms_vswf_set_spec_t *s = calloc(1,sizeof(qpms_vswf_set_spec_t));
+  if (s == NULL) return NULL; // TODO warn
+  // The rest are zeroes automatically because of calloc:
+  s->lMax_L = -1;
+  return s;
+}
+
+#define MAX(x,y) (((x)<(y))?(y):(x))
+
+qpms_errno_t qpms_vswf_set_spec_append(qpms_vswf_set_spec_t *s, const qpms_uvswfi_t u) {
+  qpms_l_t l;
+  qpms_m_t m;
+  qpms_vswf_type_t t;
+  if (qpms_uvswfi2tmn(u, &t, &m, &l)!=QPMS_SUCCESS) return QPMS_ERROR; // TODO WARN
+  if (s->n + 1 > s->capacity) {
+    size_t newcap = (s->capacity > 32) ? 32 : 2*s->capacity;
+    qpms_uvswfi_t *newmem = realloc(s->ilist, newcap * sizeof(qpms_uvswfi_t));
+    if (newmem == NULL) return QPMS_ENOMEM; // TODO WARN
+    s->ilist = newmem;
+  }
+  s->ilist[s->n] = u;
+  ++s->n;
+  switch(t) {
+    case QPMS_VSWF_ELECTRIC:
+      s->lMax_N = MAX(s->lMax_N, l);
+      break;
+    case QPMS_VSWF_MAGNETIC:
+      s->lMax_M = MAX(s->lMax_M, l);
+      break;
+    case QPMS_VSWF_LONGITUDINAL:
+      s->lMax_L = MAX(s->lMax_L, l);
+      break;
+    default:
+      abort();
+  }
+  s->lMax = MAX(s->lMax, l);
+  return QPMS_SUCCESS;
+}
+
+void qpms_vswf_set_spec_free(qpms_vswf_set_spec_t *s) {
+  if(s) free(s->ilist);
+  free(s);
+}
+
 csphvec_t qpms_vswf_single_el(qpms_m_t m, qpms_l_t l, sph_t kdlj,
     qpms_bessel_t btyp, qpms_normalisation_t norm) {
   lmcheck(l,m);

qpms/vswf.h

@@ -1,17 +1,56 @@
+/*! \file vswf.h
+ * \brief Vector spherical wavefunctions.
+ *
+ * N.B. for the Legendre polynomial norm definitions, see
+ * <a href="https://www.gnu.org/software/gsl/doc/html/specfunc.html#associated-legendre-polynomials-and-spherical-harmonics">the corresponding section of GSL docs</a>
+ * or <a href="http://git.savannah.gnu.org/cgit/gsl.git/tree/specfunc/legendre_source.c">gsl/specfunc/legendre_source.c</a>.
+ */
 #ifndef QPMS_VSWF_H
 #define QPMS_VSWF_H
 #include "qpms_types.h"
 #include <gsl/gsl_sf_legendre.h>
 
-// Electric wave N; NI
+/// Specifies a finite set of VSWFs.
+/**
+ * When for example not all the M and N -type waves up to a degree lMax
+ * need to be computed, this will specify the subset. 
+ *
+ * A typical use case would be when only even/odd fields wrt. z-plane
+ * mirror symmetry are considered.
+ */
+typedef struct qpms_vswf_set_spec_t {
+	size_t n; ///< Actual number of VSWF indices included in ilist.
+	qpms_uvswfi_t *ilist; ///< List of wave indices.
+	qpms_l_t lMax; ///< Maximum degree of the waves specified in ilist.
+	qpms_l_t lMax_M, ///< Maximum degree of the magnetic (M-type) waves.
+		 lMax_N, ///< Maximum degree of the electric (N-type) waves.
+		 lMax_L; ///< Maximum degree of the longitudinal (L-type) waves.
+	size_t capacity; ///< Allocated capacity of ilist.
+	qpms_normalisation_t norm; ///< Normalisation convention. To be set manually if needed.
+} qpms_vswf_set_spec_t;
+
+/// Creates a qpms_vswf_set_spec_t structure with an empty list of wave indices.
+qpms_vswf_set_spec_t *qpms_vswf_set_spec_init();
+
+/// Appends a VSWF index to a \ref qpms_vswf_set_spec_t, also updating metadata.
+qpms_errno_t qpms_vswf_set_spec_append(qpms_vswf_set_spec_t *self, qpms_uvswfi_t u);
+
+/// Destroys a \ref qpms_vswf_set_spec_t.
+void qpms_vswf_set_spec_free(qpms_vswf_set_spec_t *);
+
+/// Electric wave N.
 csphvec_t qpms_vswf_single_el(int m, int n, sph_t kdlj,
 		qpms_bessel_t btyp, qpms_normalisation_t norm);
-// Magnetic wave M; NI
+/// Magnetic wave M.
 csphvec_t qpms_vswf_single_mg(int m, int n, sph_t kdlj,
 		qpms_bessel_t btyp, qpms_normalisation_t norm);
 
-// Set of electric and magnetic VSWF in spherical coordinate basis
-typedef struct {
+/// Set of electric and magnetic VSWF values in spherical coordinate basis.
+/* This is supposed to contain all the waves up to $l = lMax$.
+ * for a custom set of waves, use \ref qpms_uvswfset_sph_t instead.
+ */
+
+typedef struct qpms_vswfset_sph_t {
 	//qpms_normalisation_t norm;
 	qpms_l_t lMax;
 	//qpms_y_t nelem;
@@ -20,11 +59,6 @@ typedef struct {
 } qpms_vswfset_sph_t;
 
 
-/*
- * N.B. for the norm definitions, see
- * https://www.gnu.org/software/gsl/manual/html_node/Associated-Legendre-Polynomials-and-Spherical-Harmonics.html
- * ( gsl/specfunc/legendre_source.c and 7.24.2 of gsl docs
- */
 
 qpms_errno_t qpms_legendre_deriv_y_get(double **result, double **result_deriv, double x, qpms_l_t lMax, 
 		gsl_sf_legendre_t lnorm, double csphase); // free() result and result_deriv yourself!