diff --git a/notes/Scattering and Shit.lyx b/notes/Scattering and Shit.lyx
index 1e53c39..77c4869 100644
--- a/notes/Scattering and Shit.lyx	
+++ b/notes/Scattering and Shit.lyx	
@@ -683,6 +683,19 @@ In the usual cases we encounter, the part described by the magnetic waves
  is pretty small.
 \end_layout
 
+\begin_layout Standard
+The expression with Bessel function derivatives appearing below in the electric
+ waves can be rewritten using (DLMF 10.51.2)
+\begin_inset Formula 
+\[
+\frac{1}{kr}\frac{\ud\left(kr\,z_{n}^{j}\left(kr\right)\right)}{\ud(kr)}=\frac{\ud z_{n}^{j}\left(kr\right)}{\ud(kr)}+\frac{z_{n}^{j}\left(kr\right)}{kr}=z_{n-1}^{j}\left(kr\right)-n\frac{z_{n}^{j}\left(kr\right)}{kr}.
+\]
+
+\end_inset
+
+
+\end_layout
+
 \begin_layout Subsection
 Taylor
 \end_layout
diff --git a/qpms/bessel.c b/qpms/bessel.c
new file mode 100644
index 0000000..5035a56
--- /dev/null
+++ b/qpms/bessel.c
@@ -0,0 +1,38 @@
+#include "translations.h"
+#include <stdlib.h>
+#include <gsl/gsl_sf_bessel.h>
+#include <assert.h>
+
+int qpms_sph_bessel_fill(qpms_bessel_t typ, int lmax, double x, complex double *result_array) {
+	int retval;
+	double tmparr[lmax+1];
+	switch(typ) {
+		case QPMS_BESSEL_REGULAR:
+			retval = gsl_sf_bessel_jl_steed_array(lmax, x, tmparr);
+			for (int l = 0; l <= lmax; ++l) result_array[l] = tmparr[l];
+			return retval;
+			break;
+		case QPMS_BESSEL_SINGULAR: //FIXME: is this precise enough? Would it be better to do it one-by-one?
+			retval = gsl_sf_bessel_yl_array(lmax,x,tmparr);
+			for (int l = 0; l <= lmax; ++l) result_array[l] = tmparr[l];
+			return retval;
+			break;
+		case QPMS_HANKEL_PLUS:
+		case QPMS_HANKEL_MINUS:
+			retval = gsl_sf_bessel_jl_steed_array(lmax, x, tmparr);
+			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 (typ==QPMS_HANKEL_PLUS)
+				for (int l = 0; l <= lmax; ++l) result_array[l] += I * tmparr[l];
+			else 
+				for (int l = 0; l <= lmax; ++l) result_array[l] +=-I * tmparr[l];
+			return retval;
+			break;
+		default:
+			abort();
+			//return GSL_EDOM;
+	}
+	assert(0);
+}
+
diff --git a/qpms/qpms_types.h b/qpms/qpms_types.h
index bed9fd6..e041b73 100644
--- a/qpms/qpms_types.h
+++ b/qpms/qpms_types.h
@@ -63,6 +63,10 @@ typedef struct {
 	double x, y, z;
 } cart3_t;
 
+typedef struct {
+	complex double x, y, z;
+} ccart3_t;
+
 typedef struct {
 	double x, y;
 } cart2_t;
@@ -80,4 +84,6 @@ typedef struct {
 	double r, phi;
 } pol_t;
 
+
+#define lmcheck(l,m) assert((l) >= 1 && abs(m) <= (l))
 #endif // QPMS_TYPES
diff --git a/qpms/test_vswf.c b/qpms/test_vswf_pitau.c
similarity index 100%
rename from qpms/test_vswf.c
rename to qpms/test_vswf_pitau.c
diff --git a/qpms/translations.c b/qpms/translations.c
index 8bc2422..177e799 100644
--- a/qpms/translations.c
+++ b/qpms/translations.c
@@ -1,4 +1,5 @@
 #include <math.h>
+#include "qpms_types.h"
 #include "gaunt.h"
 #include "translations.h"
 #include "indexing.h" // TODO replace size_t and int with own index types here
@@ -34,7 +35,7 @@ double qpms_legendreD0(int m, int n) {
 	return -2 * qpms_legendre0(m, n);
 }
 
-int qpms_sph_bessel_array(qpms_bessel_t typ, int lmax, double x, complex double *result_array) {
+int qpms_sph_bessel_fill(qpms_bessel_t typ, int lmax, double x, complex double *result_array) {
 	int retval;
 	double tmparr[lmax+1];
 	switch(typ) {
@@ -86,7 +87,7 @@ complex double qpms_trans_single_A_Xu(int m, int n, int mu, int nu, sph_t kdlj,
 	double leg[gsl_sf_legendre_array_n(n+nu)];
 	if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu,costheta,-1,leg)) abort();
 	complex double bes[n+nu+1];
-	if (qpms_sph_bessel_array(J, n+nu, kdlj.r, bes)) abort();
+	if (qpms_sph_bessel_fill(J, n+nu, kdlj.r, bes)) abort();
 	complex double sum = 0;
 	for(int q = 0; q <= qmax; ++q) {
 		int p = n+nu-2*q;
@@ -132,7 +133,7 @@ complex double qpms_trans_single_A_Kristensson(int m, int n, int mu, int nu, sph
 	double leg[gsl_sf_legendre_array_n(n+nu)];
 	if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu,costheta,-1,leg)) abort();
 	complex double bes[n+nu+1];
-	if (qpms_sph_bessel_array(J, n+nu, kdlj.r, bes)) abort();
+	if (qpms_sph_bessel_fill(J, n+nu, kdlj.r, bes)) abort();
 	complex double sum = 0;
 	for(int q = 0; q <= qmax; ++q) {
 		int p = n+nu-2*q;
@@ -177,7 +178,7 @@ complex double qpms_trans_single_A_Taylor(int m, int n, int mu, int nu, sph_t kd
 	double leg[gsl_sf_legendre_array_n(n+nu)];
 	if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu,costheta,-1,leg)) abort();
 	complex double bes[n+nu+1];
-	if (qpms_sph_bessel_array(J, n+nu, kdlj.r, bes)) abort();
+	if (qpms_sph_bessel_fill(J, n+nu, kdlj.r, bes)) abort();
 	complex double sum = 0;
 	for(int q = 0; q <= qmax; ++q) {
 		int p = n+nu-2*q;
@@ -229,7 +230,7 @@ complex double qpms_trans_single_B_Xu(int m, int n, int mu, int nu, sph_t kdlj,
 	double leg[gsl_sf_legendre_array_n(n+nu+1)];
 	if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1,costheta,-1,leg)) abort();
 	complex double bes[n+nu+2];
-	if (qpms_sph_bessel_array(J, n+nu+1, kdlj.r, bes)) abort();
+	if (qpms_sph_bessel_fill(J, n+nu+1, kdlj.r, bes)) abort();
 
 	complex double sum = 0;
 	for (int q = 0; q <= Qmax; ++q) {
@@ -285,7 +286,7 @@ complex double qpms_trans_single_B_Taylor(int m, int n, int mu, int nu, sph_t kd
 	double leg[gsl_sf_legendre_array_n(n+nu+1)];
 	if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1,costheta,-1,leg)) abort();
 	complex double bes[n+nu+2];
-	if (qpms_sph_bessel_array(J, n+nu+1, kdlj.r, bes)) abort();
+	if (qpms_sph_bessel_fill(J, n+nu+1, kdlj.r, bes)) abort();
 
 	complex double sum = 0;
 	for (int q = 0; q <= Qmax; ++q) {
@@ -338,10 +339,6 @@ void qpms_trans_calculator_free(qpms_trans_calculator *c) {
 	free(c);
 }
 
-static inline size_t qpms_mn2y(int m, int n) {
-	return (size_t) n * (n + 1) + m - 1;
-}
-
 static inline size_t qpms_trans_calculator_index_mnmunu(const qpms_trans_calculator *c,
 		int m, int n, int mu, int nu){
 	return c->nelem * qpms_mn2y(m,n) + qpms_mn2y(mu,nu);
@@ -470,9 +467,9 @@ static void qpms_trans_calculator_multipliers_B_Taylor(
 	}
 }
 
-int qpms_trans_calculator_multipliers_A(qpms_normalization_t norm, complex double *dest, int m, int n, int mu, int nu, int qmax) {
+int qpms_trans_calculator_multipliers_A(qpms_normalisation_t norm, complex double *dest, int m, int n, int mu, int nu, int qmax) {
 	switch (norm) {
-		case QPMS_NORMALIZATION_TAYLOR:
+		case QPMS_NORMALISATION_TAYLOR:
 			qpms_trans_calculator_multipliers_A_Taylor(dest,m,n,mu,nu,qmax);
 			return 0;
 			break;
@@ -482,9 +479,9 @@ int qpms_trans_calculator_multipliers_A(qpms_normalization_t norm, complex doubl
 	assert(0);
 }
 
-int qpms_trans_calculator_multipliers_B(qpms_normalization_t norm, complex double *dest, int m, int n, int mu, int nu, int qmax) {
+int qpms_trans_calculator_multipliers_B(qpms_normalisation_t norm, complex double *dest, int m, int n, int mu, int nu, int qmax) {
 	switch (norm) {
-		case QPMS_NORMALIZATION_TAYLOR:
+		case QPMS_NORMALISATION_TAYLOR:
 			qpms_trans_calculator_multipliers_B_Taylor(dest,m,n,mu,nu,qmax);
 			return 0;
 			break;
@@ -495,14 +492,14 @@ int qpms_trans_calculator_multipliers_B(qpms_normalization_t norm, complex doubl
 }
 
 qpms_trans_calculator
-*qpms_trans_calculator_init (int lMax, qpms_normalization_t normalization) {
+*qpms_trans_calculator_init (int lMax, qpms_normalisation_t normalisation) {
 	assert(lMax > 0);
 	qpms_trans_calculator *c = malloc(sizeof(qpms_trans_calculator));
 	c->lMax = lMax;
 	c->nelem = lMax * (lMax+2);
 	c->A_multipliers = malloc((1+SQ(c->nelem)) * sizeof(complex double *));
 	c->B_multipliers = malloc((1+SQ(c->nelem)) * sizeof(complex double *));
-	c->normalization = normalization;
+	c->normalisation = normalisation;
 	size_t *qmaxes = malloc(SQ(c->nelem) * sizeof(size_t));
 	size_t qmaxsum = 0;
 	for(size_t y = 0; y < c->nelem; y++)
@@ -525,7 +522,7 @@ qpms_trans_calculator
 			int m, n, mu, nu;
 			qpms_y2mn_p(y, &m, &n);
 			qpms_y2mn_p(yu, &mu, &nu);
-			qpms_trans_calculator_multipliers_A(normalization,
+			qpms_trans_calculator_multipliers_A(normalisation,
 					c->A_multipliers[i], m, n, mu, nu, qmaxes[i]);
 		}
 
@@ -550,7 +547,7 @@ qpms_trans_calculator
 			int m, n, mu, nu;
 			qpms_y2mn_p(y, &m, &n);
 			qpms_y2mn_p(yu, &mu, &nu);
-			qpms_trans_calculator_multipliers_B(normalization,
+			qpms_trans_calculator_multipliers_B(normalisation,
 					c->B_multipliers[i], m, n, mu, nu, qmaxes[i]);
 		}
 
@@ -586,13 +583,13 @@ complex double qpms_trans_calculator_get_A_buf(const qpms_trans_calculator *c,
 	if (0 == kdlj.r && J != QPMS_BESSEL_REGULAR) 
 		// TODO warn? 
 		return NAN+I*NAN;
-	switch(c->normalization) {
-		case QPMS_NORMALIZATION_TAYLOR:
+	switch(c->normalisation) {
+		case QPMS_NORMALISATION_TAYLOR:
 			{
 				double costheta = cos(kdlj.theta);
 				if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu,
 							costheta,-1,legendre_buf)) abort();
-				if (qpms_sph_bessel_array(J, n+nu+1, kdlj.r, bessel_buf)) abort();
+				if (qpms_sph_bessel_fill(J, n+nu+1, kdlj.r, bessel_buf)) abort();
 				return qpms_trans_calculator_get_A_precalcbuf(c,m,n,mu,nu,
 						kdlj,r_ge_d,J,bessel_buf,legendre_buf);
 			}
@@ -631,13 +628,13 @@ complex double qpms_trans_calculator_get_B_buf(const qpms_trans_calculator *c,
 	if (0 == kdlj.r && J != QPMS_BESSEL_REGULAR) 
 		// TODO warn? 
 		return NAN+I*NAN;
-	switch(c->normalization) {
-		case QPMS_NORMALIZATION_TAYLOR:
+	switch(c->normalisation) {
+		case QPMS_NORMALISATION_TAYLOR:
 			{
 				double costheta = cos(kdlj.theta);
 				if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1,
 							costheta,-1,legendre_buf)) abort();
-				if (qpms_sph_bessel_array(J, n+nu+2, kdlj.r, bessel_buf)) abort();
+				if (qpms_sph_bessel_fill(J, n+nu+2, kdlj.r, bessel_buf)) abort();
 				return qpms_trans_calculator_get_B_precalcbuf(c,m,n,mu,nu,
 						kdlj,r_ge_d,J,bessel_buf,legendre_buf);
 			}
@@ -660,13 +657,13 @@ int qpms_trans_calculator_get_AB_buf_p(const qpms_trans_calculator *c,
 		// TODO warn? different return value?
 		return 0;
 	}
-	switch(c->normalization) {
-		case QPMS_NORMALIZATION_TAYLOR:
+	switch(c->normalisation) {
+		case QPMS_NORMALISATION_TAYLOR:
 			{
 				double costheta = cos(kdlj.theta);
 				if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1,
 							costheta,-1,legendre_buf)) abort();
-				if (qpms_sph_bessel_array(J, n+nu+2, kdlj.r, bessel_buf)) abort();
+				if (qpms_sph_bessel_fill(J, n+nu+2, kdlj.r, bessel_buf)) abort();
 				*Adest = qpms_trans_calculator_get_A_precalcbuf(c,m,n,mu,nu,
 						kdlj,r_ge_d,J,bessel_buf,legendre_buf);
 				*Bdest = qpms_trans_calculator_get_B_precalcbuf(c,m,n,mu,nu,
@@ -698,13 +695,13 @@ int qpms_trans_calculator_get_AB_arrays_buf(const qpms_trans_calculator *c,
 		// TODO warn? different return value?
 		return 0;
 	}
-	switch(c->normalization) {
-		case QPMS_NORMALIZATION_TAYLOR:
+	switch(c->normalisation) {
+		case QPMS_NORMALISATION_TAYLOR:
 			{
 				double costheta = cos(kdlj.theta);
 				if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,2*c->lMax+1,
 							costheta,-1,legendre_buf)) abort();
-				if (qpms_sph_bessel_array(J, 2*c->lMax+2, kdlj.r, bessel_buf)) abort();
+				if (qpms_sph_bessel_fill(J, 2*c->lMax+2, kdlj.r, bessel_buf)) abort();
 				size_t desti = 0, srci = 0;
 				for (int n = 1; n <= c->lMax; ++n) for (int m = -n; m <= n; ++m) {
 					for (int nu = 1; nu <= c->lMax; ++nu) for (int mu = -nu; mu <= nu; ++mu) {
diff --git a/qpms/translations.h b/qpms/translations.h
index f63a288..6e69e85 100644
--- a/qpms/translations.h
+++ b/qpms/translations.h
@@ -9,7 +9,8 @@
 double qpms_legendre0(int m, int n); 
 // Associated Legendre polynomial derivative at zero argument (DLMF 14.5.2)
 double qpms_legendred0(int m, int n); 
-int qpms_sph_bessel_array(qpms_bessel_t typ, int lmax, double x, complex double *result_array);
+// TODO unify types
+int qpms_sph_bessel_fill(qpms_bessel_t typ, int lmax, double x, complex double *result_array);
 
 // TODO replace the xplicit "Taylor" functions with general,
 // taking qpms_bessel_t argument.
@@ -30,10 +31,10 @@ typedef struct qpms_trans_calculator {
         size_t nelem;
         complex double **A_multipliers;
         complex double **B_multipliers;
-        qpms_normalization_t normalization;
+        qpms_normalisation_t normalisation;
 } qpms_trans_calculator;
 
-qpms_trans_calculator *qpms_trans_calculator_init(int lMax, qpms_normalization_t nt);
+qpms_trans_calculator *qpms_trans_calculator_init(int lMax, qpms_normalisation_t nt);
 void qpms_trans_calculator_free(qpms_trans_calculator *);
 
 complex double qpms_trans_calculator_get_A(const qpms_trans_calculator *c,
diff --git a/qpms/vectors.h b/qpms/vectors.h
index a998741..1724f23 100644
--- a/qpms/vectors.h
+++ b/qpms/vectors.h
@@ -6,11 +6,11 @@
 
 //static inline double vectors_h_sq(double x) {return x*x;}
 
-static inline double cart3norm(cart3_t v) {
+static inline double cart3norm(const cart3_t v) {
 	return sqrt(v.x*v.x + v.y*v.y + v.z*v.z);
 }
 
-static inline sph_t cart2sph(cart3_t cart) {
+static inline sph_t cart2sph(const cart3_t cart) {
 	sph_t sph;
 	sph.r = cart3norm(cart);
 	sph.theta = sph.r ? acos(cart.z / sph.r) : M_PI_2;
@@ -18,7 +18,7 @@ static inline sph_t cart2sph(cart3_t cart) {
 	return sph;
 }
 
-static inline cart3_t sph2cart(sph_t sph) {
+static inline cart3_t sph2cart(const sph_t sph) {
 	cart3_t cart;
 	double sin_th = sin(sph.theta);
 	cart.x = sph.r * sin_th * cos(sph.phi);
@@ -27,5 +27,28 @@ static inline cart3_t sph2cart(sph_t sph) {
 	return cart;
 }
 
+// equivalent to sph_loccart2cart in qpms_p.py
+static inline ccart3_t csphvec2cart(const csphvec_t sphvec, const sph_t at) {
+	ccart3_t res = {0, 0, 0};
+	const double st = sin(at.theta);
+	const double ct = cos(at.theta);
+	const double sf = sin(at.phi);
+	const double cf = cos(at.phi);
+	const double rx = st * cf;
+	const double ry = st * sf;
+	const double rz = ct;
+	const double tx = ct * cf;
+	const double ty = ct * sf;
+	const double tz = -st;
+	const double fx = -sf;
+	const double fy = cf;
+	const double fz = 0.;
+	res.x = rx * sphvec.rc + tx * sphvec.thetac + fx * sphvec.phic;
+	res.y = ry * sphvec.rc + ty * sphvec.thetac + fy * sphvec.phic;
+	res.z = rz * sphvec.rc + tz * sphvec.thetac + fz * sphvec.phic;
+	return sphvec;
+}
+
+
 
 #endif //VECTORS_H
diff --git a/qpms/vswf.c b/qpms/vswf.c
index 4837490..f49e700 100644
--- a/qpms/vswf.c
+++ b/qpms/vswf.c
@@ -1,10 +1,12 @@
 #include "vswf.h"
 #include "indexing.h"
+#include "translations.h" // TODO move qpms_sph_bessel_fill elsewhere
 #include <math.h>
 #include <gsl/gsl_math.h>
 #include <gsl/gsl_sf_legendre.h>
 #include <stdlib.h>
 #include <string.h>
+#include "assert_cython_workaround.h"
 
 // Legendre functions also for negative m, see DLMF 14.9.3
 qpms_errno_t qpms_legendre_deriv_y_fill(double *target, double *target_deriv, double x, qpms_l_t lMax,
@@ -156,3 +158,93 @@ void qpms_pitau_free(qpms_pitau_t x) {
 	free(x.pi);
 	free(x.tau);
 }
+
+
+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);
+	csphvec_t N;
+	complex double *bessel = malloc((l+1)*sizeof(complex double));
+	if(qpms_sph_bessel_fill(btyp, l, kdlj.r, bessel)) abort();
+	qpms_pitau_t pt = qpms_pitau_get(kdlj.theta, l, norm);
+	complex double eimf = cexp(m * kdlj.phi * I);
+	qpms_y_t y = qpms_mn2y(m,l);
+
+	N.rc = l*(l+1) * pt.leg[y] * bessel[l] / kdlj.r * eimf;
+	complex double besselfac = bessel[l-1] - l * bessel[l] / kdlj.r;
+	N.thetac = pt.tau[y] * besselfac * eimf;
+	N.phic = pt.pi[y] * besselfac * I * eimf;
+
+	qpms_pitau_free(pt);
+	free(bessel);
+	return N;
+}
+csphvec_t qpms_vswf_single_mg(qpms_m_t m, qpms_l_t l, sph_t kdlj,
+		qpms_bessel_t btyp, qpms_normalisation_t norm) {
+	lmcheck(l,m);
+	csphvec_t M;
+	complex double *bessel = malloc((l+1)*sizeof(complex double));
+	if(qpms_sph_bessel_fill(btyp, l, kdlj.r, bessel)) abort();
+	qpms_pitau_t pt = qpms_pitau_get(kdlj.theta, l, norm);
+	complex double eimf = cexp(m * kdlj.phi * I);
+	qpms_y_t y = qpms_mn2y(m,l);
+
+	M.rc = 0.;
+	M.thetac = pt.pi[y] * bessel[l] * I * eimf;
+	M.phic = -pt.tau[y] * bessel[l] * eimf;
+
+	qpms_pitau_free(pt);
+	free(bessel);
+	return M;
+}
+
+qpms_vswfset_sph_t *qpms_vswfset_make(qpms_l_t lMax, sph_t kdlj, 
+		qpms_bessel_t btyp, qpms_normalisation_t norm) {
+	qpms_vswfset_sph_t *res = malloc(sizeof(qpms_vswfset_sph_t));
+	res->lMax = lMax;
+	qpms_y_t nelem = qpms_lMax2nelem(lMax);
+	res->el = malloc(sizeof(csphvec_t)*nelem);
+	res->mg = malloc(sizeof(csphvec_t)*nelem);
+	if(QPMS_SUCCESS != qpms_vswf_fill(res->mg, res->el, lMax, kdlj, btyp, norm))
+		abort(); // or return NULL? or rather assert?
+	return res;
+}
+
+void qpms_vswfset_sph_pfree(qpms_vswfset_sph_t *w) {
+	assert(NULL != w && NULL != w->el && NULL != w->mg);
+	free(w->el);
+	free(w->mg);
+	free(w);
+}
+
+qpms_errno_t qpms_vswf_fill(csphvec_t *mgtarget, csphvec_t *eltarget,
+		qpms_l_t lMax, sph_t kr,
+		qpms_bessel_t btyp, qpms_normalisation_t norm) {
+	assert(lMax >= 1);
+	complex double *bessel = malloc((lMax+1)*sizeof(complex double));
+	if(qpms_sph_bessel_fill(btyp, lMax, kr.r, bessel)) abort();
+	qpms_pitau_t pt = qpms_pitau_get(kr.theta, lMax, norm);
+	complex double const *pbes = bessel + 1; // starting from l = 1
+	double const *pleg = pt.leg;
+	double const *ppi = pt.pi;
+	double const *ptau = pt.tau;
+	csphvec_t *pmg = mgtarget, *pel = eltarget;
+	for(qpms_l_t l = 1; l <= lMax; ++l) {
+		complex double besfac = *pbes / kr.r;
+		complex double besderfac = *(pbes-1) - l * besfac;
+		for(qpms_m_t m = -l; m <= l; ++m) {
+			complex double eimf = cexp(m * kr.phi * I);
+			pel->rc = l*(l+1) * (*pleg) * besfac * eimf;
+			pel->thetac = *ptau * besderfac * eimf;
+			pel->phic = *ppi * besderfac * I * eimf;
+			pmg->rc = 0.;
+			pmg->thetac = *ppi * (*pbes) * I * eimf;
+			pmg->phic = - *ptau * (*pbes) * eimf;
+			++pleg; ++ppi; ++ptau; ++pel; ++pmg;
+		}
+		++pbes;
+	}
+	free(bessel);
+	qpms_pitau_free(pt);
+	return QPMS_SUCCESS;
+}
diff --git a/qpms/vswf.h b/qpms/vswf.h
index 852b540..d76d8cb 100644
--- a/qpms/vswf.h
+++ b/qpms/vswf.h
@@ -31,6 +31,9 @@ qpms_errno_t qpms_legendre_deriv_y_get(double **result, double **result_deriv, d
 qpms_errno_t qpms_legendre_deriv_y_fill(double *where, double *where_deriv, double x, 
 		qpms_l_t lMax, gsl_sf_legendre_t lnorm, double csphase); 
 
+qpms_errno_t qpms_vswf_fill(csphvec_t *resultM, csphvec_t *resultN, qpms_l_t lMax, sph_t kdrj,
+		qpms_bessel_t btyp, qpms_normalisation_t norm);
+
 qpms_vswfset_sph_t *qpms_vswfset_make(qpms_l_t lMax, sph_t kdlj,
 	qpms_bessel_t btyp, qpms_normalisation_t norm);//NI
 void qpms_vswfset_sph_pfree(qpms_vswfset_sph_t *);//NI