diff --git a/vec_wrappers.py b/vec_wrappers.py
index e3963a9b250fe3ed6f7022c613ea62d127ca9a26..f30f1713ee7f5b4805bec60f75e813a3c9771b8f 100644
--- a/vec_wrappers.py
+++ b/vec_wrappers.py
@@ -468,39 +468,52 @@ def gen_vector_wrappers():
 
     for dp_or_no in ["", "_dp"]:
         for dims in [2, 3]:
-            for eqn in [cgh.Laplace(dims), cgh.Helmholtz(dims)]:
-                func_name = "%s%ddformta%s" % (eqn.lh_letter(), dims,  dp_or_no)
-                gen_vector_wrapper(func_name,
-                Template("""
-                        integer ier(nvcount)
-                        % if eqn.lh_letter() == "h":
-                            complex*16 zk
-                        % endif
-                        real*8 rscale(nvcount)
-                        real *8 sources(${dims},*INDIRECT_MANY)
-
-                        % if dp_or_no:
-                            complex *16 dipstr(*INDIRECT_MANY)
-                            %if not (eqn.lh_letter() == "l" and dims == 2):
-                                real *8 dipvec(${dims}, *INDIRECT_MANY)
-                            %endif
-                        % else:
-                            complex *16 charge(*INDIRECT_MANY)
-                        % endif
+            if dims == 3:
+                trunc_or_no = ["", "_trunc"]
+            else:
+                trunc_or_no = [""]
 
-                        integer nsources(*INDIRECT_MANY)
-                        real*8 center(${dims}, nvcount)
-                        integer nterms
-                        complex*16 expn(${eqn.expansion_dims("nterms")},nvcount)
-                        """, strict_undefined=True).render(
-                            dims=dims,
-                            eqn=eqn,
-                            dp_or_no=dp_or_no,
-                            ),
-                        ["ier", "expn"],
-                        output_reductions={"expn": "sum", "ier": "max"},
-                        tmp_init={"ier": "0"},
-                        vec_func_name=func_name + "_imany")
+            for trunc in trunc_or_no:
+                for eqn in [cgh.Laplace(dims), cgh.Helmholtz(dims)]:
+                    func_name = "%s%ddformta%s%s" % (eqn.lh_letter(), dims,
+                            dp_or_no, trunc)
+                    gen_vector_wrapper(func_name,
+                    Template("""
+                            integer ier(nvcount)
+                            % if eqn.lh_letter() == "h":
+                                complex*16 zk
+                            % endif
+                            real*8 rscale(nvcount)
+                            real *8 sources(${dims},*INDIRECT_MANY)
+
+                            % if dp_or_no:
+                                complex *16 dipstr(*INDIRECT_MANY)
+                                %if not (eqn.lh_letter() == "l" and dims == 2):
+                                    real *8 dipvec(${dims}, *INDIRECT_MANY)
+                                %endif
+                            % else:
+                                complex *16 charge(*INDIRECT_MANY)
+                            % endif
+
+                            integer nsources(*INDIRECT_MANY)
+                            real*8 center(${dims}, nvcount)
+                            integer nterms
+                            complex*16 expn(${eqn.expansion_dims("nterms")},nvcount)
+
+                            % if trunc:
+                                real *8 wlege(0:nlege,0:nlege)
+                                integer nlege
+                            % endif
+                            """, strict_undefined=True).render(
+                                dims=dims,
+                                eqn=eqn,
+                                dp_or_no=dp_or_no,
+                                trunc=trunc,
+                                ),
+                            ["ier", "expn"],
+                            output_reductions={"expn": "sum", "ier": "max"},
+                            tmp_init={"ier": "0"},
+                            vec_func_name=func_name + "_imany")
 
     # }}}
 
diff --git a/wrappers.pyf.mako b/wrappers.pyf.mako
index f10dd8e25de32265cf224c755720157c2e347e27..defcfd19732818ea3fcf932d53d7f62b17005b4b 100644
--- a/wrappers.pyf.mako
+++ b/wrappers.pyf.mako
@@ -7,6 +7,8 @@
 
 python module _internal
   interface
+    ! {{{ special functions
+
       subroutine jfuns2d(ier,nterms,z,scale,fjs,ifder,fjder, &
           lwfjs,iscale,ntop)
         ! implicit real *8 (a-h,o-z)
@@ -22,8 +24,6 @@ python module _internal
         complex*16, intent(out) :: ntop
     end subroutine
 
-    ! {{{ special functions
-
     subroutine legeexps(itype,n,x,u,v,whts)
       integer, intent(in) :: itype, n
       real*8, intent(out) :: x(n),whts(n),u(n,n),v(n,n)
@@ -43,6 +43,11 @@ python module _internal
       dimension ts(n),whts(n)
     end subroutine
 
+    subroutine ylgndrini(nmax, rat1, rat2)
+      integer, intent(in) :: nmax
+      real *8, intent(out) :: rat1(0:nmax,0:nmax), rat2(0:nmax,0:nmax)
+    end subroutine
+
     ! }}}
 
     ! {{{ form{mp,ta} entrypoints