reorganisation of the code

Former-commit-id: 2efe6a07fa27aa2a159dfa83959e9580b3f53809

qpms/bessel.c

@@ -1,9 +1,9 @@
 #include <assert.h>
-#include "translations.h"
+#include "qpms_specfunc.h"
 #include <stdlib.h>
 #include <gsl/gsl_sf_bessel.h>
 
-int qpms_sph_bessel_fill(qpms_bessel_t typ, int lmax, double x, complex double *result_array) {
+qpms_errno_t qpms_sph_bessel_fill(qpms_bessel_t typ, qpms_l_t lmax, double x, complex double *result_array) {
 	int retval;
 	double tmparr[lmax+1];
 	switch(typ) {

qpms/legendre.c

@@ -0,0 +1,158 @@
+#include "qpms_specfunc.h"
+#include "qpms_types.h"
+#include <gsl/gsl_sf_legendre.h>
+#include <gsl/gsl_math.h>
+#include <stdlib.h>
+#include "indexing.h"
+#include <string.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,
+                gsl_sf_legendre_t lnorm, double csphase)
+{
+        size_t n = gsl_sf_legendre_array_n(lMax);
+        	double *legendre_tmp = malloc(n * sizeof(double));
+        double *legendre_deriv_tmp = malloc(n * sizeof(double));
+        int gsl_errno = gsl_sf_legendre_deriv_array_e(
+                        lnorm, (size_t)lMax, x, csphase, legendre_tmp,legendre_deriv_tmp);
+        for (qpms_l_t l = 1; l <= lMax; ++l)
+                        for (qpms_m_t m = 0; m <= l; ++m) {
+                                qpms_y_t y = qpms_mn2y(m,l);
+                                size_t i = gsl_sf_legendre_array_index(l,m);
+                                target[y] = legendre_tmp[i];
+                                target_deriv[y] = legendre_deriv_tmp[i];
+                        }
+        switch(lnorm) {
+                case GSL_SF_LEGENDRE_NONE:
+                        for (qpms_l_t l = 1; l <= lMax; ++l)
+                                for (qpms_m_t m = 1; m <= l; ++m) {
+                                qpms_y_t y = qpms_mn2y(-m,l);
+                                size_t i = gsl_sf_legendre_array_index(l,m);
+                                // viz DLMF 14.9.3, čert ví, jak je to s cs fasí.
+                                double factor = exp(lgamma(l-m+1)-lgamma(l+m+1))*((m%2)?-1:1);
+                                target[y] = factor * legendre_tmp[i];
+                                target_deriv[y] = factor * legendre_deriv_tmp[i];
+                                }
+                        break;
+                case GSL_SF_LEGENDRE_SCHMIDT:
+                case GSL_SF_LEGENDRE_SPHARM:
+                case GSL_SF_LEGENDRE_FULL:
+                        for (qpms_l_t l = 1; l <= lMax; ++l)
+                                for (qpms_m_t m = 1; m <= l; ++m) {
+                                qpms_y_t y = qpms_mn2y(-m,l);
+                                size_t i = gsl_sf_legendre_array_index(l,m);
+                                // viz DLMF 14.9.3, čert ví, jak je to s cs fasí.
+                                double factor = ((m%2)?-1:1); // this is the difference from the unnormalised case
+                                target[y] = factor * legendre_tmp[i];
+                                target_deriv[y] = factor * legendre_deriv_tmp[i];
+                                }
+                        break;
+                default:
+                        abort(); //NI
+                        break;
+        }
+        free(legendre_tmp);
+        free(legendre_deriv_tmp);
+        return QPMS_SUCCESS;
+}
+
+qpms_errno_t qpms_legendre_deriv_y_get(double **target, double **dtarget, double x, qpms_l_t lMax, gsl_sf_legendre_t lnorm,
+                double csphase)
+{
+
+        *target = malloc(sizeof(double)*qpms_lMax2nelem(lMax));
+        *dtarget = malloc(sizeof(double)*qpms_lMax2nelem(lMax));
+        return qpms_legendre_deriv_y_fill(*target, *dtarget, x, lMax, lnorm, csphase);
+}
+
+qpms_pitau_t qpms_pitau_get(double theta, qpms_l_t lMax, qpms_normalisation_t norm)
+{
+        const double csphase = qpms_normalisation_t_csphase(norm);
+        norm = qpms_normalisation_t_normonly(norm);
+        qpms_pitau_t res;
+        qpms_y_t nelem = qpms_lMax2nelem(lMax);
+        res.pi = malloc(nelem * sizeof(double));
+        res.tau = malloc(nelem * sizeof(double));
+        double ct = cos(theta), st = sin(theta);
+        if (1 == fabs(ct)) { // singular case, use DLMF 14.8.2
+                memset(res.pi, 0, nelem*sizeof(double));
+                memset(res.tau, 0, nelem*sizeof(double));
+                res.leg = calloc(nelem, sizeof(double));
+                switch(norm) {
+                        case QPMS_NORMALISATION_XU:
+                                for (qpms_l_t l = 1; l <= lMax; ++l) {
+                                        res.leg[qpms_mn2y(0, l)] = (l%2)?ct:1.;
+                                        double p = l*(l+1)/2;
+                                        const double n = 0.5;
+                                        int lpar = (l%2)?-1:1;
+                                        res.pi [qpms_mn2y(+1, l)] = -((ct>0) ? -1 : lpar) * p * csphase;
+                                        res.pi [qpms_mn2y(-1, l)] = -((ct>0) ? -1 : lpar) * n * csphase;
+                                        res.tau[qpms_mn2y(+1, l)] = ((ct>0) ? +1 : lpar) * p * csphase;
+                                        res.tau[qpms_mn2y(-1, l)] = -((ct>0) ? +1 : lpar) * n * csphase;
+                                }
+                                break;
+                        case QPMS_NORMALISATION_TAYLOR:
+                                for (qpms_l_t l = 1; l <= lMax; ++l) {
+                                        res.leg[qpms_mn2y(0, l)] = ((l%2)?ct:1.)*sqrt((2*l+1)*0.25*M_1_PI);
+                                        int lpar = (l%2)?-1:1;
+                                        double fl = 0.25 * sqrt((2*l+1)*l*(l+1)*M_1_PI);
+                                        res.pi [qpms_mn2y(+1, l)] = -((ct>0) ? -1 : lpar) * fl * csphase;
+                                        res.pi [qpms_mn2y(-1, l)] = -((ct>0) ? -1 : lpar) * fl * csphase;
+                                        res.tau[qpms_mn2y(+1, l)] = ((ct>0) ? +1 : lpar) * fl * csphase;
+                                        res.tau[qpms_mn2y(-1, l)] = -((ct>0) ? +1 : lpar) * fl * csphase;
+                                }
+                                break;
+                        case QPMS_NORMALISATION_POWER:
+                                for (qpms_l_t l = 1; l <= lMax; ++l) {
+                                        res.leg[qpms_mn2y(0, l)] = ((l%2)?ct:1.)*sqrt((2*l+1)/(4*M_PI *l*(l+1)));
+                                        int lpar = (l%2)?-1:1;
+                                        double fl = 0.25 * sqrt((2*l+1)*M_1_PI);
+                                        res.pi [qpms_mn2y(+1, l)] = -((ct>0) ? -1 : lpar) * fl * csphase;
+                                        res.pi [qpms_mn2y(-1, l)] = -((ct>0) ? -1 : lpar) * fl * csphase;
+                                        res.tau[qpms_mn2y(+1, l)] = ((ct>0) ? +1 : lpar) * fl * csphase;
+                                        res.tau[qpms_mn2y(-1, l)] = -((ct>0) ? +1 : lpar) * fl * csphase;
+
+                                }
+                                break;
+                        default:
+                                abort();
+                }
+        }
+        else { // cos(theta) in (-1,1), use normal calculation
+                double *legder = malloc(sizeof(double)*qpms_lMax2nelem(lMax));
+                res.leg = malloc(sizeof(double)*qpms_lMax2nelem(lMax));
+                if (qpms_legendre_deriv_y_fill(res.leg, legder, ct, lMax,
+                                        norm == QPMS_NORMALISATION_XU ? GSL_SF_LEGENDRE_NONE
+                                                : GSL_SF_LEGENDRE_SPHARM, csphase))
+                        abort();
+                if (norm == QPMS_NORMALISATION_POWER)
+                        /* for Xu (=non-normalized) and Taylor (=sph. harm. normalized)
+                         * the correct normalisation is already obtained from gsl_sf_legendre_deriv_array_e().
+                         * However, Kristensson ("power") normalisation differs from Taylor
+                         * by 1/sqrt(l*(l+1)) factor.
+                         */
+                        for (qpms_l_t l = 1; l <= lMax; ++l) {
+                                double prefac = 1./sqrt(l*(l+1));
+                                for (qpms_m_t m = -l; m <= l; ++m) {
+                                        res.leg[qpms_mn2y(m,l)] *= prefac;
+                                        legder[qpms_mn2y(m,l)] *= prefac;
+                                }
+                        }
+                for (qpms_l_t l = 1; l <= lMax; ++l) {
+                        for (qpms_m_t m = -l; m <= l; ++m) {
+                                                res.pi [qpms_mn2y(m,l)] = m / st * res.leg[qpms_mn2y(m,l)];
+                                                res.tau[qpms_mn2y(m,l)] = - st * legder[qpms_mn2y(m,l)];
+                        }
+                }
+                free(legder);
+        }
+        res.lMax = lMax;
+        return res;
+}
+
+void qpms_pitau_free(qpms_pitau_t x) {
+        free(x.leg);
+        free(x.pi);
+        free(x.tau);
+}
+

qpms/qpms_specfunc.h

@@ -0,0 +1,54 @@
+#ifndef QPMS_SPECFUNC_H
+#define QPMS_SPECFUNC_H
+#include "qpms_types.h"
+#include <gsl/gsl_sf_legendre.h>
+
+/******************************************************************************
+ *                    Spherical Bessel functions                              *
+ ******************************************************************************/
+
+// TODO unify types
+qpms_errno_t qpms_sph_bessel_fill(qpms_bessel_t typ, qpms_l_t lmax, double x, complex double *result_array);
+
+
+/******************************************************************************
+ *         Legendre functions and their "angular derivatives"                 *
+ ******************************************************************************/
+
+/*
+ * 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!
+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); 
+
+
+double *qpms_legendre_y_get(double x, qpms_l_t lMax, qpms_normalisation_t norm);//NI
+double *qpms_legendre0d_y_get(qpms_l_t lMax, qpms_normalisation_t norm); //NI
+double *qpms_legendre_plus1d_y_get(qpms_l_t lMax, qpms_normalisation_t norm); //NI
+double *qpms_legendre_minus1d_y_get(qpms_l_t lMax, qpms_normalisation_t norm); //NI
+
+
+
+// array of Legendre and pi, tau auxillary functions (see [1,(37)])
+// This should handle correct evaluation for theta -> 0 and theta -> pi
+typedef struct {
+	//qpms_normalisation_t norm;
+	qpms_l_t lMax;
+	//qpms_y_t nelem;
+	double *leg, *pi, *tau;
+} qpms_pitau_t;
+qpms_pitau_t qpms_pitau_get(double theta, qpms_l_t lMax, qpms_normalisation_t norm);
+void qpms_pitau_free(qpms_pitau_t);//NI
+void qpms_pitau_pfree(qpms_pitau_t*);//NI
+
+// Associated Legendre polynomial at zero argument (DLMF 14.5.1) DEPRECATED?
+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);
+
+#endif // QPMS_SPECFUNC_H

qpms/qpms_types.h

@@ -57,7 +57,9 @@ static inline int qpms_normalisation_t_normonly(qpms_normalisation_t norm) {
 }
 
 
-// TODO move elsewhere
+// TODO move the inlines elsewhere
+#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) {

qpms/translations.h

@@ -5,65 +5,71 @@
 #include <complex.h>
 #include <stdbool.h>
 #include <stddef.h>
-// Associated Legendre polynomial at zero argument (DLMF 14.5.1)
-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); 
-// 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.
-complex double qpms_trans_single_A_Taylor(int m, int n, int mu, int nu, sph_t kdlj,
+complex double qpms_trans_single_A_Taylor(qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, sph_t kdlj,
                 bool r_ge_d, qpms_bessel_t J);
 
-complex double qpms_trans_single_B_Taylor(int m, int n, int mu, int nu, sph_t kdlj,
+complex double qpms_trans_single_B_Taylor(qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, sph_t kdlj,
                 bool r_ge_d, qpms_bessel_t J);
 
-complex double qpms_trans_single_A_Taylor_ext(int m, int n, int mu, int nu, double kdlj_r,
+complex double qpms_trans_single_A_Taylor_ext(qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, double kdlj_r,
 		double kdlj_th, double kdlj_phi, int r_ge_d, int J);
 
-complex double qpms_trans_single_B_Taylor_ext(int m, int n, int mu, int nu, double kdlj_r,
+complex double qpms_trans_single_B_Taylor_ext(qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, double kdlj_r,
 		double kdlj_th, double kdlj_phi, int r_ge_d, int J);
 
 typedef struct qpms_trans_calculator {
-        int lMax;
+        qpms_normalisation_t normalisation;
-        size_t nelem;
+        qpms_l_t lMax;
+        qpms_y_t nelem;
         complex double **A_multipliers;
         complex double **B_multipliers;
-        qpms_normalisation_t normalisation;
+#if 0
+	// Normalised values of the Legendre functions and derivatives
+	// for θ == π/2, i.e. for the 2D case.
+	double *leg0; 
+	double *pi0;
+	double *tau0;
+	// Spherical Bessel function coefficients:
+	// TODO
+#endif
 } qpms_trans_calculator;
 
-qpms_trans_calculator *qpms_trans_calculator_init(int lMax, qpms_normalisation_t nt);
+qpms_trans_calculator *qpms_trans_calculator_init(qpms_l_t 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,
-		int m, int n, int mu, int nu, sph_t kdlj,
+		qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, sph_t kdlj,
 		bool r_ge_d, qpms_bessel_t J);
 complex double qpms_trans_calculator_get_B(const qpms_trans_calculator *c,
-		int m, int n, int mu, int nu, sph_t kdlj,
+		qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, sph_t kdlj,
 		bool r_ge_d, qpms_bessel_t J);
+int qpms_trans_calculator_get_AB_p(const qpms_trans_calculator *c,
+		complex double *Adest, complex double *Bdest,
+		qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, sph_t kdlj,
+		bool r_ge_d, qpms_bessel_t J);
+int qpms_trans_calculator_get_AB_arrays(const qpms_trans_calculator *c,
+		complex double *Adest, complex double *Bdest,
+		size_t deststride, size_t srcstride,
+		sph_t kdlj, bool r_ge_d, qpms_bessel_t J); 
 
+
+// TODO update the types later
 complex double qpms_trans_calculator_get_A_ext(const qpms_trans_calculator *c,
 		int m, int n, int mu, int nu, double kdlj_r,
 		double kdlj_th, double kdlj_phi, int r_ge_d, int J);
+
 complex double qpms_trans_calculator_get_B_ext(const qpms_trans_calculator *c,
 		int m, int n, int mu, int nu, double kdlj_r,
 		double kdlj_th, double kdlj_phi, int r_ge_d, int J);
 
-int qpms_trans_calculator_get_AB_p(const qpms_trans_calculator *c,
-		complex double *Adest, complex double *Bdest,
-		int m, int n, int mu, int nu, sph_t kdlj,
-		bool r_ge_d, qpms_bessel_t J);
 int qpms_trans_calculator_get_AB_p_ext(const qpms_trans_calculator *c,
 		complex double *Adest, complex double *Bdest,
 		int m, int n, int mu, int nu, double kdlj_r,
 		double kdlj_th, double kdlj_phi, int r_ge_d, int J);
 
-int qpms_trans_calculator_get_AB_arrays(const qpms_trans_calculator *c,
-		complex double *Adest, complex double *Bdest,
-		size_t deststride, size_t srcstride,
-		sph_t kdlj, bool r_ge_d, qpms_bessel_t J); 
 int qpms_trans_calculator_get_AB_arrays_ext(const qpms_trans_calculator *c,
 		complex double *Adest, complex double *Bdest,
 		size_t deststride, size_t srcstride,

qpms/vswf.c

@@ -4,162 +4,10 @@
 #include "vswf.h"
 #include "indexing.h"
 #include "translations.h" // TODO move qpms_sph_bessel_fill elsewhere
-#include <gsl/gsl_sf_legendre.h>
+#include "qpms_specfunc.h"
 #include <stdlib.h>
 #include <string.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,
-		gsl_sf_legendre_t lnorm, double csphase)
-{
-	size_t n = gsl_sf_legendre_array_n(lMax);
-	double *legendre_tmp = malloc(n * sizeof(double));
-	double *legendre_deriv_tmp = malloc(n * sizeof(double));
-	int gsl_errno = gsl_sf_legendre_deriv_array_e(
-			lnorm, (size_t)lMax, x, csphase, legendre_tmp,legendre_deriv_tmp);
-	for (qpms_l_t l = 1; l <= lMax; ++l)
-			for (qpms_m_t m = 0; m <= l; ++m) {
-				qpms_y_t y = qpms_mn2y(m,l);
-				size_t i = gsl_sf_legendre_array_index(l,m);
-				target[y] = legendre_tmp[i];
-				target_deriv[y] = legendre_deriv_tmp[i];
-			}
-	switch(lnorm) {
-		case GSL_SF_LEGENDRE_NONE:
-			for (qpms_l_t l = 1; l <= lMax; ++l)
-				for (qpms_m_t m = 1; m <= l; ++m) {
-				qpms_y_t y = qpms_mn2y(-m,l);
-				size_t i = gsl_sf_legendre_array_index(l,m);
-				// viz DLMF 14.9.3, čert ví, jak je to s cs fasí.
-				double factor = exp(lgamma(l-m+1)-lgamma(l+m+1))*((m%2)?-1:1);
-				target[y] = factor * legendre_tmp[i];
-				target_deriv[y] = factor * legendre_deriv_tmp[i];
-				}
-			break;
-		case GSL_SF_LEGENDRE_SCHMIDT:
-		case GSL_SF_LEGENDRE_SPHARM:
-		case GSL_SF_LEGENDRE_FULL:
-			for (qpms_l_t l = 1; l <= lMax; ++l)
-				for (qpms_m_t m = 1; m <= l; ++m) {
-				qpms_y_t y = qpms_mn2y(-m,l);
-				size_t i = gsl_sf_legendre_array_index(l,m);
-				// viz DLMF 14.9.3, čert ví, jak je to s cs fasí.
-				double factor = ((m%2)?-1:1); // this is the difference from the unnormalised case
-				target[y] = factor * legendre_tmp[i];
-				target_deriv[y] = factor * legendre_deriv_tmp[i];
-				}
-			break;
-		default:
-			abort(); //NI
-			break;
-	}
-	free(legendre_tmp);
-	free(legendre_deriv_tmp);
-	return QPMS_SUCCESS;
-}
-
-qpms_errno_t qpms_legendre_deriv_y_get(double **target, double **dtarget, double x, qpms_l_t lMax, gsl_sf_legendre_t lnorm,
-		double csphase)
-{
-
-	*target = malloc(sizeof(double)*qpms_lMax2nelem(lMax));
-	*dtarget = malloc(sizeof(double)*qpms_lMax2nelem(lMax));
-	return qpms_legendre_deriv_y_fill(*target, *dtarget, x, lMax, lnorm, csphase);
-}
-
-
-qpms_pitau_t qpms_pitau_get(double theta, qpms_l_t lMax, qpms_normalisation_t norm) 
-{
-	const double csphase = qpms_normalisation_t_csphase(norm);
-	norm = qpms_normalisation_t_normonly(norm);
-	qpms_pitau_t res;
-	qpms_y_t nelem = qpms_lMax2nelem(lMax);
-	res.pi = malloc(nelem * sizeof(double));
-	res.tau = malloc(nelem * sizeof(double));
-	double ct = cos(theta), st = sin(theta);
-	if (1 == fabs(ct)) { // singular case, use DLMF 14.8.2
-		memset(res.pi, 0, nelem*sizeof(double));
-		memset(res.tau, 0, nelem*sizeof(double));
-		res.leg = calloc(nelem, sizeof(double));
-		switch(norm) {
-			case QPMS_NORMALISATION_XU:
-				for (qpms_l_t l = 1; l <= lMax; ++l) {
-					res.leg[qpms_mn2y(0, l)] = (l%2)?ct:1.;
-					double p = l*(l+1)/2;
-					const double n = 0.5;
-					int lpar = (l%2)?-1:1;
-					res.pi [qpms_mn2y(+1, l)] = -((ct>0) ? -1 : lpar) * p * csphase;
-					res.pi [qpms_mn2y(-1, l)] = -((ct>0) ? -1 : lpar) * n * csphase;
-					res.tau[qpms_mn2y(+1, l)] = ((ct>0) ? +1 : lpar) * p * csphase;
-					res.tau[qpms_mn2y(-1, l)] = -((ct>0) ? +1 : lpar) * n * csphase;
-				}
-				break;
-			case QPMS_NORMALISATION_TAYLOR:
-				for (qpms_l_t l = 1; l <= lMax; ++l) {
-					res.leg[qpms_mn2y(0, l)] = ((l%2)?ct:1.)*sqrt((2*l+1)*0.25*M_1_PI);
-					int lpar = (l%2)?-1:1;
-					double fl = 0.25 * sqrt((2*l+1)*l*(l+1)*M_1_PI);
-					res.pi [qpms_mn2y(+1, l)] = -((ct>0) ? -1 : lpar) * fl * csphase;
-					res.pi [qpms_mn2y(-1, l)] = -((ct>0) ? -1 : lpar) * fl * csphase;
-					res.tau[qpms_mn2y(+1, l)] = ((ct>0) ? +1 : lpar) * fl * csphase;
-					res.tau[qpms_mn2y(-1, l)] = -((ct>0) ? +1 : lpar) * fl * csphase;
-				}
-				break;
-			case QPMS_NORMALISATION_POWER:
-				for (qpms_l_t l = 1; l <= lMax; ++l) {
-					res.leg[qpms_mn2y(0, l)] = ((l%2)?ct:1.)*sqrt((2*l+1)/(4*M_PI *l*(l+1)));
-					int lpar = (l%2)?-1:1;
-					double fl = 0.25 * sqrt((2*l+1)*M_1_PI);
-					res.pi [qpms_mn2y(+1, l)] = -((ct>0) ? -1 : lpar) * fl * csphase;
-					res.pi [qpms_mn2y(-1, l)] = -((ct>0) ? -1 : lpar) * fl * csphase;
-					res.tau[qpms_mn2y(+1, l)] = ((ct>0) ? +1 : lpar) * fl * csphase;
-					res.tau[qpms_mn2y(-1, l)] = -((ct>0) ? +1 : lpar) * fl * csphase;
-	
-				}
-				break;
-			default:
-				abort();
-		}
-	}
-	else { // cos(theta) in (-1,1), use normal calculation
-		double *legder = malloc(sizeof(double)*qpms_lMax2nelem(lMax));
-		res.leg = malloc(sizeof(double)*qpms_lMax2nelem(lMax));
-		if (qpms_legendre_deriv_y_fill(res.leg, legder, ct, lMax, 
-					norm == QPMS_NORMALISATION_XU ? GSL_SF_LEGENDRE_NONE
-						: GSL_SF_LEGENDRE_SPHARM, csphase)) 
-			abort();
-		if (norm == QPMS_NORMALISATION_POWER) 
-			/* for Xu (=non-normalized) and Taylor (=sph. harm. normalized) 
-			 * the correct normalisation is already obtained from gsl_sf_legendre_deriv_array_e().
-			 * However, Kristensson ("power") normalisation differs from Taylor
-			 * by 1/sqrt(l*(l+1)) factor.
-			 */
-			for (qpms_l_t l = 1; l <= lMax; ++l) {
-				double prefac = 1./sqrt(l*(l+1));
-				for (qpms_m_t m = -l; m <= l; ++m) {
-					res.leg[qpms_mn2y(m,l)] *= prefac;
-					legder[qpms_mn2y(m,l)] *= prefac;
-				}
-			}
-		for (qpms_l_t l = 1; l <= lMax; ++l) {
-			for (qpms_m_t m = -l; m <= l; ++m) {
-						res.pi [qpms_mn2y(m,l)] = m / st * res.leg[qpms_mn2y(m,l)];
-						res.tau[qpms_mn2y(m,l)] = - st * legder[qpms_mn2y(m,l)];
-			}
-		}
-		free(legder);
-	}
-	res.lMax = lMax;
-	return res;
-}
-
-void qpms_pitau_free(qpms_pitau_t x) {
-	free(x.leg);
-	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);

qpms/vswf.h

@@ -60,23 +60,4 @@ 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
 
-double *qpms_legendre_y_get(double x, qpms_l_t lMax, qpms_normalisation_t norm);//NI
-double *qpms_legendre0d_y_get(qpms_l_t lMax, qpms_normalisation_t norm); //NI
-double *qpms_legendre_plus1d_y_get(qpms_l_t lMax, qpms_normalisation_t norm); //NI
-double *qpms_legendre_minus1d_y_get(qpms_l_t lMax, qpms_normalisation_t norm); //NI
-
-
-
-// array of Legendre and pi, tau auxillary functions (see [1,(37)])
-// This should handle correct evaluation for theta -> 0 and theta -> pi
-typedef struct {
-	//qpms_normalisation_t norm;
-	qpms_l_t lMax;
-	//qpms_y_t nelem;
-	double *leg, *pi, *tau;
-} qpms_pitau_t;
-qpms_pitau_t qpms_pitau_get(double theta, qpms_l_t lMax, qpms_normalisation_t norm);
-void qpms_pitau_free(qpms_pitau_t);//NI
-void qpms_pitau_pfree(qpms_pitau_t*);//NI
-
 #endif // QPMS_VSWF_H