diff --git a/qpms/bessel.c b/qpms/bessel.c
index 1d11805..4b3edb1 100644
--- a/qpms/bessel.c
+++ b/qpms/bessel.c
@@ -59,6 +59,11 @@ qpms_errno_t qpms_sph_bessel_realx_fill(qpms_bessel_t typ, qpms_l_t lmax, double
       return retval;
       break;
     case QPMS_BESSEL_SINGULAR: //FIXME: is this precise enough? Would it be better to do it one-by-one?
+      if(QPMS_UNLIKELY(x == 0)) {
+        for (int l = 0; l <= lmax; ++l)
+          result_array[l] = NAN;
+        return QPMS_ESING; // GSL would have returned GSL_EDOM without setting NANs.
+      }
       retval = gsl_sf_bessel_yl_array(lmax,x,tmparr);
       for (int l = 0; l <= lmax; ++l) result_array[l] = tmparr[l];
       return retval;
@@ -69,6 +74,11 @@ qpms_errno_t qpms_sph_bessel_realx_fill(qpms_bessel_t typ, qpms_l_t lmax, double
       for (int l = 0; l <= lmax; ++l) result_array[l] = tmparr[l];
       if(retval) return retval;
       retval = gsl_sf_bessel_yl_array(lmax, x, tmparr);
+      if(QPMS_UNLIKELY(x == 0)) {
+        for (int l = 0; l <= lmax; ++l)
+          result_array[l] += I * NAN;
+        return QPMS_ESING; // GSL would have returned GSL_EDOM without setting NANs.
+      }
       if (typ==QPMS_HANKEL_PLUS)
         for (int l = 0; l <= lmax; ++l) result_array[l] += I * tmparr[l];
       else 
@@ -76,10 +86,9 @@ qpms_errno_t qpms_sph_bessel_realx_fill(qpms_bessel_t typ, qpms_l_t lmax, double
       return retval;
       break;
     default:
-      abort();
-      //return GSL_EDOM;
+      QPMS_INVALID_ENUM(typ);
   }
-  assert(0);
+  QPMS_WTF;
 }
 
 // TODO DOC
@@ -119,7 +128,7 @@ qpms_errno_t qpms_sph_bessel_fill(qpms_bessel_t typ, qpms_l_t lmax, complex doub
         }
         break;
       default:
-        QPMS_WTF;
+        QPMS_INVALID_ENUM(typ);
     }
     // TODO check for underflows? (nz != 0)
     if (ierr == 0 || ierr == 3) {
diff --git a/qpms/qpms_error.h b/qpms/qpms_error.h
index e40e105..14b68b9 100644
--- a/qpms/qpms_error.h
+++ b/qpms/qpms_error.h
@@ -83,6 +83,8 @@ qpms_dbgmsg_flags qpms_dbgmsg_enable(qpms_dbgmsg_flags types);
 
 #define QPMS_WTF qpms_pr_error_at_flf(__FILE__,__LINE__,__func__,"Unexpected error.")
 
+#define QPMS_INVALID_ENUM(x) qpms_pr_error_at_flf(__FILE__,__LINE__,__func__,"Invalid enum value (" #x " == %d)", (int) (x))
+
 #define QPMS_UNTESTED {\
 	static _Bool already_bitched = 0; \
 	if (QPMS_UNLIKELY(!already_bitched)) {\
@@ -93,9 +95,9 @@ qpms_dbgmsg_flags qpms_dbgmsg_enable(qpms_dbgmsg_flags types);
 
 #define QPMS_PR_ERROR(msg, ...) qpms_pr_error_at_flf(__FILE__,__LINE__,__func__,msg, ##__VA_ARGS__)
 
-// TODO replace this macro with an inline function?
+/// Raises an error if \a x is not zero.
 #define QPMS_ENSURE_SUCCESS(x) { \
-	int errorcode = (x); \
+	int errorcode = (x); /* evaluate x only once */ \
 	if(QPMS_UNLIKELY(errorcode)) \
 		qpms_pr_error_at_flf(__FILE__,__LINE__,__func__,"Unexpected error (%d)", errorcode); \
 }
@@ -107,6 +109,17 @@ qpms_dbgmsg_flags qpms_dbgmsg_enable(qpms_dbgmsg_flags types);
 		qpms_pr_error_at_flf(__FILE__,__LINE__,__func__,msg, ##__VA_ARGS__); \
 }
 
+/// Raises an error if \a x is not zero or one of the values listed in arguments.
+#define QPMS_ENSURE_SUCCESS_OR(x, ...) { \
+	int errorcode = (x); /* evaluate x only once */ \
+	static const int allowed_errorcodes[] = {0, ##__VA_ARGS__};\
+	static const int n_allowed = sizeof(allowed_errorcodes) / sizeof(int); \
+	int i; \
+	for(i = 0; i < n_allowed; ++i) \
+		if (errorcode == allowed_errorcodes[i]) break; \
+	if (QPMS_UNLIKELY(i >= n_allowed))  \
+		qpms_pr_error_at_flf(__FILE__,__LINE__,__func__,"Unexpected error (%d)", errorcode); \
+}
 
 #define QPMS_ENSURE(x, msg, ...) {if(QPMS_UNLIKELY(!(x))) qpms_pr_error_at_flf(__FILE__,__LINE__,__func__,msg, ##__VA_ARGS__); }
 
diff --git a/qpms/qpms_types.h b/qpms/qpms_types.h
index 8b169b2..1d8f7dd 100644
--- a/qpms/qpms_types.h
+++ b/qpms/qpms_types.h
@@ -78,13 +78,14 @@ typedef enum {
 typedef unsigned long long qpms_uvswfi_t; 
 
 /// Error codes / return values for certain numerical functions.
-/** These are de facto a subset of the GSL error codes. */
+/** Those with values between -2 and 32 are subset of the GSL error codes. */
 typedef enum {
 	QPMS_SUCCESS = 0, ///< Success.
-	QPMS_ERROR = 1, ///< Unspecified error.
+	QPMS_ERROR = -1, ///< Unspecified error.
+	QPMS_FAILURE = QPMS_ERROR,
 	QPMS_ENOMEM = 8, ///< Out of memory.
+	QPMS_ESING = 21, ///< Apparent singularity detected.
 	QPMS_NAN_ENCOUNTERED = 1024 ///< NaN value encountered in data processing.
-
 } qpms_errno_t;
 
 
diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c
index 17a5e47..bc6d60b 100644
--- a/qpms/scatsystem.c
+++ b/qpms/scatsystem.c
@@ -2039,7 +2039,7 @@ ccart3_t qpms_scatsys_scattered_E__alt(const qpms_scatsys_t *ss,
       qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC,  0, 1),
       qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, +1, 1),
     },
-    .lMax=1, .lMax_M=1, .lMax_N=1, .lMax_L=-1,
+    .lMax=1, .lMax_M=0, .lMax_N=1, .lMax_L=-1,
     .capacity=0,
     .norm = ss->c->normalisation,
   };
diff --git a/qpms/vswf.c b/qpms/vswf.c
index 19a87d0..f827364 100644
--- a/qpms/vswf.c
+++ b/qpms/vswf.c
@@ -231,7 +231,8 @@ qpms_errno_t qpms_vswf_fill_csph(csphvec_t *const longtarget,
     csph_t kr, qpms_bessel_t btyp, const qpms_normalisation_t norm) {
   assert(lMax >= 1);
   complex double *bessel = malloc((lMax+1)*sizeof(complex double));
-  QPMS_ENSURE_SUCCESS(qpms_sph_bessel_fill(btyp, lMax, kr.r, bessel));
+  qpms_errno_t bessel_retval = qpms_sph_bessel_fill(btyp, lMax, kr.r, bessel);
+  QPMS_ENSURE_SUCCESS_OR(bessel_retval, QPMS_ESING);
   qpms_pitau_t pt = qpms_pitau_get(kr.theta, lMax, qpms_normalisation_t_csphase(norm));
   complex double const *pbes = bessel + 1; // starting from l = 1
   double const *pleg = pt.leg;
@@ -282,7 +283,7 @@ qpms_errno_t qpms_vswf_fill_csph(csphvec_t *const longtarget,
   }
   free(bessel);
   qpms_pitau_free(pt);
-  return QPMS_SUCCESS;
+  return bessel_retval;
 }
 
 qpms_errno_t qpms_vswf_fill(csphvec_t *const longtarget, 
diff --git a/qpms/vswf.h b/qpms/vswf.h
index 33dae0e..d2d01ac 100644
--- a/qpms/vswf.h
+++ b/qpms/vswf.h
@@ -177,6 +177,11 @@ csphvec_t qpms_vswf_L00(
  * 
  * Does not evaluate the zeroth-order wave \f$ \mathbf{L}_0^0 \f$. 
  * If you need that, use qpms_vswf_L00().
+ *
+ * If \a kdrj.r == 0 and \a btyp != QPMS_BESSEL_REGULAR, returns QPMS_ESING and fills
+ * targets with NaNs.
+ *
+ * \return 0 (QPMS_SUCCESS) on success, otherwise.
  */
 qpms_errno_t qpms_vswf_fill(
 	csphvec_t *resultL, //< Target array for longitudinal VSWFs.
diff --git a/tests/test_vswf_transop_consistency.py b/tests/test_vswf_transop_consistency.py
index d6ea7ba..2674109 100755
--- a/tests/test_vswf_transop_consistency.py
+++ b/tests/test_vswf_transop_consistency.py
@@ -6,6 +6,7 @@ waves at origin
 """
 
 from qpms import Particle, CTMatrix, lorentz_drude, EpsMuGenerator, TMatrixGenerator, BaseSpec, FinitePointGroup, ScatteringSystem, TMatrixInterpolator, EpsMu, dbgmsg_enable, dbgmsg_disable, dbgmsg_active, BesselType,eV, hbar
+from qpms.qpms_c import set_gsl_pythonic_error_handling
 from qpms.symmetries import point_group_info
 import numpy as np
 eh = eV/hbar
@@ -15,6 +16,8 @@ dbgmsg_enable(2)
 part_radius = 80e-9
 p = 1580e-9
 
+set_gsl_pythonic_error_handling()
+
 sym = FinitePointGroup(point_group_info['D4h'])
 bspec1 = BaseSpec(lMax=3)
 medium=EpsMuGenerator(EpsMu(1.52**2))
@@ -36,7 +39,7 @@ for i in range(ss.fecv_size):
     E = ssw.scattered_E(fvc, points)
     E_alt = ssw.scattered_E(fvc, points,alt=True)
     diff = abs(E-E_alt)
-    reldiffavg = np.average(diff/(abs(E)+abs(E_alt)))
+    reldiffavg = np.nanmean(diff/(abs(E)+abs(E_alt)))
     fail = reldiffavg > 1e-3
     fails += fail