diff --git a/qpms/__init__.py b/qpms/__init__.py
index cc43bed..5f73da3 100644
--- a/qpms/__init__.py
+++ b/qpms/__init__.py
@@ -7,7 +7,7 @@ from sys import platform as __platform
 import warnings as __warnings
 
 try:
-    from .qpms_c import PointGroup, FinitePointGroup, FinitePointGroupElement, Particle, scatsystem_set_nthreads, ScatteringSystem, ScatteringMatrix
+    from .qpms_c import PointGroup, FinitePointGroup, FinitePointGroupElement, Particle, scatsystem_set_nthreads, ScatteringSystem, ScatteringMatrix, pitau
 except ImportError as ex:
     if __platform == "linux" or __platform == "linux2":
         if 'LD_LIBRARY_PATH' not in __os.environ.keys():
diff --git a/qpms/legendre.c b/qpms/legendre.c
index 9d3f19a..a440897 100644
--- a/qpms/legendre.c
+++ b/qpms/legendre.c
@@ -53,50 +53,67 @@ qpms_errno_t qpms_legendre_deriv_y_get(double **target, double **dtarget, double
 
 qpms_pitau_t qpms_pitau_get(double theta, qpms_l_t lMax, const double csphase)
 {
-  QPMS_ENSURE(fabs(csphase) == 1, "The csphase argument must be either 1 or -1, not %g.", csphase);
   qpms_pitau_t res;
   const qpms_y_t nelem = qpms_lMax2nelem(lMax);
   QPMS_CRASHING_MALLOC(res.leg, nelem * sizeof(double));
   QPMS_CRASHING_MALLOC(res.pi, nelem * sizeof(double));
   QPMS_CRASHING_MALLOC(res.tau, nelem * sizeof(double));
+  qpms_pitau_fill(res.leg, res.pi, res.tau, theta, lMax, csphase);
+  return res;
+}
+
+qpms_errno_t qpms_pitau_fill(double *target_leg, double *pi, double *tau, double theta, qpms_l_t lMax, double csphase)
+{
+  QPMS_ENSURE(fabs(csphase) == 1, "The csphase argument must be either 1 or -1, not %g.", csphase);
+  const qpms_y_t nelem = qpms_lMax2nelem(lMax);
+ 
   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));
-    memset(res.leg, 0, nelem * sizeof(double));
+    if(pi) memset(pi, 0, nelem * sizeof(double));
+    if(tau) memset(tau, 0, nelem * sizeof(double));
+    if(target_leg) memset(target_leg, 0, nelem * sizeof(double));
     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)));
+      if(target_leg) target_leg[qpms_mn2y(0, l)] = ((l%2)?ct:1.)*sqrt((2*l+1)/(4*M_PI *l*(l+1)));
       double fl = 0.25 * sqrt((2*l+1)*M_1_PI);
       int lpar = (l%2)?-1:1;
-      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;
+      if(pi) {
+        pi [qpms_mn2y(+1, l)] = -((ct>0) ? -1 : lpar) * fl * csphase;
+        pi [qpms_mn2y(-1, l)] = -((ct>0) ? -1 : lpar) * fl * csphase;
+      }
+      if(tau) {
+        tau[qpms_mn2y(+1, l)] = ((ct>0) ? +1 : lpar) * fl * csphase;
+        tau[qpms_mn2y(-1, l)] = -((ct>0) ? +1 : lpar) * fl * csphase;
+      }
     }
   }
   else { // cos(theta) in (-1,1), use normal calculation
-    double *legder;
+    double *leg, *legder;
+    if (target_leg)
+      leg = target_leg;
+    else
+      QPMS_CRASHING_MALLOC(leg, nelem*sizeof(double));
     QPMS_CRASHING_MALLOC(legder, nelem * sizeof(double));
-    QPMS_ENSURE_SUCCESS(qpms_legendre_deriv_y_fill(res.leg, legder, ct, lMax,
+    QPMS_ENSURE_SUCCESS(qpms_legendre_deriv_y_fill(leg, legder, ct, lMax,
           GSL_SF_LEGENDRE_SPHARM, csphase));
     // Multiply by the "power normalisation" 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;
+        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)];
+        if(pi)  pi [qpms_mn2y(m,l)] = m / st * leg[qpms_mn2y(m,l)];
+        if(tau) tau[qpms_mn2y(m,l)] = - st * legder[qpms_mn2y(m,l)];
       }
     }
     free(legder);
+    if(!target_leg)
+      free(leg);
   }
-  res.lMax = lMax;
-  return res;
+  return QPMS_SUCCESS;
 }
 
 void qpms_pitau_free(qpms_pitau_t x) {
diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx
index cd23411..58941cf 100644
--- a/qpms/qpms_c.pyx
+++ b/qpms/qpms_c.pyx
@@ -488,4 +488,15 @@ cdef class ScatteringMatrix:
         qpms_scatsys_scatter_solve(&f_view[0], &a_view[0], self.lu)
         return f
 
-
+def pitau(double theta, qpms_l_t lMax, double csphase = 1):
+    if(abs(csphase) != 1):
+        raise ValueError("csphase must be 1 or -1, is %g" % csphase)
+    cdef size_t nelem = qpms_lMax2nelem(lMax)
+    cdef np.ndarray[np.float_t, ndim=1] lega = np.empty((nelem,), dtype=float)
+    cdef np.ndarray[np.float_t, ndim=1] pia = np.empty((nelem,), dtype=float)
+    cdef np.ndarray[np.float_t, ndim=1] taua = np.empty((nelem,), dtype=float)
+    cdef double[::1] leg = lega
+    cdef double[::1] pi = pia
+    cdef double[::1] tau = taua
+    qpms_pitau_fill(&leg[0], &pi[0], &tau[0], theta, lMax, csphase)
+    return (lega, pia, taua)
diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd
index 0aee2a3..fc7d7b6 100644
--- a/qpms/qpms_cdefs.pxd
+++ b/qpms/qpms_cdefs.pxd
@@ -65,6 +65,7 @@ cdef extern from "qpms_types.h":
     ctypedef int qpms_lm_t
     ctypedef int qpms_l_t
     ctypedef int qpms_m_t
+    ctypedef size_t qpms_y_t
     struct qpms_quat_t:
         cdouble a
         cdouble b
@@ -161,6 +162,7 @@ cdef extern from "vswf.h":
             const cart3_t evalpoint, const void *args, bint add)
 
 cdef extern from "indexing.h":
+    qpms_y_t qpms_lMax2nelem(qpms_l_t lMax)
     qpms_uvswfi_t qpms_tmn2uvswfi(qpms_vswf_type_t t, qpms_m_t m, qpms_l_t n)
     qpms_errno_t qpms_uvswfi2tmn(qpms_uvswfi_t u, qpms_vswf_type_t* t, qpms_m_t* m, qpms_l_t* n)
     qpms_m_t qpms_uvswfi2m(qpms_uvswfi_t u)
@@ -313,6 +315,17 @@ cdef extern from "translations.h":
                 qpms_bessel_t J
                 );
 
+cdef extern from "qpms_specfunc.h":
+    struct qpms_pitau_t:
+        qpms_l_t lMax
+        double *leg
+        double *pi
+        double *tau
+    qpms_pitau_t qpms_pitau_get(double theta, qpms_l_t lMax, double csphase)
+    qpms_errno_t qpms_pitau_fill(double *leg, double *pi, double *tau, double theta, qpms_l_t lMax, double csphase)
+    void qpms_pitau_free(qpms_pitau_t pitau)
+
+
 cdef extern from "gsl/gsl_interp.h":
     struct gsl_interp_type:
         pass
diff --git a/qpms/qpms_specfunc.h b/qpms/qpms_specfunc.h
index 9db24dd..3c07a91 100644
--- a/qpms/qpms_specfunc.h
+++ b/qpms/qpms_specfunc.h
@@ -68,7 +68,7 @@ typedef struct {
 	double *leg, *pi, *tau;
 } qpms_pitau_t;
 
-/// Returns an array of normalised Legendre and and auxillary \f$\pi_{lm}, \tau_{lm} \f$ functions.
+/// Returns an array of normalised Legendre and auxillary \f$\pi_{lm}, \tau_{lm} \f$ functions.
 /**
  * The normalised Legendre function here is defined as
  * \f[ 
@@ -88,8 +88,18 @@ typedef struct {
  *
  */
 qpms_pitau_t qpms_pitau_get(double theta, qpms_l_t lMax, double csphase);
+
+/// Directly fills (pre-allocated) arrays of normalised Legendre and auxillary \f$\pi_{lm}, \tau_{lm} \f$ functions.
+/**
+ * Arrays must be preallocated for `lMax * (lMax + 2)` elements. `NULL` targets are skipped. 
+ * For details, see qpms_pitau_get().
+ */
+qpms_errno_t qpms_pitau_fill(double *target_leg, double *target_pi, double *target_tau,
+		double theta, qpms_l_t lMax, double csphase);
+
 /// Frees the dynamically allocated arrays from qpms_pitau_t.
 void qpms_pitau_free(qpms_pitau_t);
+
 //void qpms_pitau_pfree(qpms_pitau_t*);//NI
 
 // Associated Legendre polynomial at zero argument (DLMF 14.5.1) DEPRECATED?