diff --git a/eigensystems.lisp b/eigensystems.lisp
index fc3e5685a136513e3208c0a9fe84504ed7644f15..51d9bd20d4fcc26226b19ab16bb3bbd318109ff4 100644
--- a/eigensystems.lisp
+++ b/eigensystems.lisp
@@ -1,12 +1,16 @@
 ;; Eigenvectors and eigenvalues
 ;; Liam Healy, Sun May 21 2006 - 19:52
-;; Time-stamp: <2009-01-08 21:49:15EST eigensystems.lisp>
+;; Time-stamp: <2009-01-11 10:35:06EST eigensystems.lisp>
 ;; $Id$
 
 (in-package :gsl)
 
 ;;; /usr/include/gsl/gsl_eigen.h
 
+;;; Should symmetric matrices form a subclass of matrices, so that
+;;; both eigenvalues and eigenvalues-nonsymm could be methods of the
+;;; same function?
+
 ;;;;****************************************************************************
 ;;;; Workspace 
 ;;;;****************************************************************************
@@ -18,6 +22,23 @@
   n-by-n real symmetric matrices.  The size of the workspace
   is O(2n).")
 
+;; V 1.9
+(defmobject eigen-nonsymm
+    "gsl_eigen_nonsymm" ((n sizet))
+    "non-symmetric eigenvalue workspace"	; FDL
+    "Make a workspace for computing eigenvalues of
+  n-by-n real non-symmetric matrices.  The size of the workspace
+  is O(2n).")
+
+(cffi:defcstruct gsl-nonsymm-ws
+  (size sizet)		     ; size of matrices
+  (diag :pointer)	     ; diagonal matrix elements from balancing
+  (tau :pointer)	     ; Householder coefficients
+  (Z :pointer)		     ; pointer to Z matrix
+  (balancep :int)	     ; perform balancing transformation?
+  (n-evals sizet)	     ; number of eigenvalues found
+  (francis-ws :pointer))
+
 (defmobject eigen-symmv
     "gsl_eigen_symmv" ((n sizet))
     "symmetric eigensystem workspace"	; FDL
@@ -25,6 +46,14 @@
   eigenvectors of n-by-n real symmetric matrices.  The size of
   the workspace is O(4n).")
 
+;; V 1.9
+(defmobject eigen-nonsymmv
+    "gsl_eigen_nonsymmv" ((n sizet))
+    "non-symmetric eigenvalue workspace"	; FDL
+    "Make a workspace for computing for computing eigenvalues and
+    eigenvectors of n-by-n real nonsymmetric matrices. The size of the
+    workspace is O(5n).")
+
 (defmobject eigen-herm
     "gsl_eigen_herm" ((n sizet))
     "Hermitian eigenvalue workspace"	; FDL
@@ -68,6 +97,81 @@
   referenced.  The eigenvalues are stored in the vector eigenvalues and
   are unordered.")
 
+(defmfun set-parameters-nonsymmetric
+    (ws &optional compute-shur-form balance)
+  "gsl_eigen_nonsymm_params"
+  (((if compute-shur-form 1 0) :int)
+   ((if balance 1 0) :int) ((mpointer ws) :pointer))
+  :gsl-version (1 9)
+  :c-return :void
+  :export nil
+  :index eigenvalues-nonsymm)
+
+(defmfun eigenvalues-nonsymm
+    (A
+     &optional
+     (eigenvalues
+      (make-marray '(complex double-float) :dimensions (dim0 A)))
+     (ws (make-eigen-nonsymm (dim0 A)))
+     compute-shur-form balance shur-vectors
+     &aux
+     (sv
+      (if (eql shur-vectors t)
+	  (make-marray 'double-float :dimensions (dimensions A))
+	  shur-vectors)))
+  ("gsl_eigen_nonsymm" "gsl_eigen_nonsymm_Z")
+  ((((mpointer A) :pointer)
+    ((mpointer eigenvalues) :pointer) ((mpointer ws) :pointer))
+   (((mpointer A) :pointer)
+    ((mpointer eigenvalues) :pointer) ((mpointer ws) :pointer)
+    ((mpointer sv) :pointer)))
+  :before
+  ((set-parameters-nonsymmetric ws compute-shur-form balance))
+  :gsl-version (1 9)
+  :switch (shur-vectors)
+  :inputs (A)
+  :outputs (A eigenvalues)
+  :return
+  (eigenvalues
+   (cffi:foreign-slot-value (mpointer ws) 'gsl-nonsymm-ws 'n-evals))
+  :documentation			; FDL
+  "Compute the eigenvalues of the real nonsymmetric matrix A and
+  stores them in the vector 'eigenvalues. If T is desired, it is
+  stored in the upper portion of A on output. Otherwise, on output,
+  the diagonal of A will contain the 1-by-1 real eigenvalues and
+  2-by-2 complex conjugate eigenvalue systems, and the rest of A is
+  destroyed. In rare cases, this function may fail to find all
+  eigenvalues. If this happens, a warning is signalled and the number
+  of converged eigenvalues is returned as a second value. The
+  converged eigenvalues are stored in the beginning of eval.
+
+  If compute-shur-form is true, the full Schur form T will be computed.
+  If it is set to nil, T will not be computed (this is
+  the default setting). Computing the full Schur form requires
+  approximately 1.5-2 times the number of flops.
+
+  If balance is true, a balancing transformation is applied to the
+  matrix prior to computing eigenvalues. This transformation is
+  designed to make the rows and columns of the matrix have comparable
+  norms, and can result in more accurate eigenvalues for matrices
+  whose entries vary widely in magnitude. See Balancing for more
+  information. Note that the balancing transformation does not
+  preserve the orthogonality of the Schur vectors, so if you wish to
+  compute the Schur vectors with you will obtain the Schur vectors of
+  the balanced matrix instead of the original matrix. The relationship
+  will be
+
+          T = Q^t D^(-1) A D Q
+
+  where Q is the matrix of Schur vectors for the balanced matrix, and D
+  is the balancing transformation. Then this function will compute
+  a matrix Z which satisfies
+
+          T = Z^(-1) A Z
+
+  with Z = D Q. Note that Z will not be orthogonal. For this reason,
+  balancing is not performed by default.")
+
 (defmfun eigenvalues-eigenvectors
     ((A matrix)
      &optional
@@ -98,6 +202,75 @@
   The eigenvectors are guaranteed to be mutually orthogonal and
   normalised to unit magnitude.")
 
+(defmfun eigenvalues-eigenvectors-nonsymm
+    (A
+     &optional
+     (eigenvalues
+      (make-marray '(complex double-float) :dimensions (dim0 A)))
+     (eigenvectors
+      (make-marray  '(complex double-float) :dimensions (dimensions A)))
+     (ws (make-eigen-nonsymmv (dim0 A)))
+     compute-shur-form balance shur-vectors
+     &aux
+     (sv
+      (if (eql shur-vectors t)
+	  (make-marray 'double-float :dimensions (dimensions A))
+	  shur-vectors)))
+  ("gsl_eigen_nonsymmv" "gsl_eigen_nonsymmv_Z")
+  ((((mpointer A) :pointer)
+    ((mpointer eigenvalues) :pointer) ((mpointer eigenvectors) :pointer)
+    ((mpointer ws) :pointer))
+   (((mpointer A) :pointer)
+    ((mpointer eigenvalues) :pointer) ((mpointer eigenvectors) :pointer)
+    ((mpointer ws) :pointer) ((mpointer sv) :pointer)))
+  :before		      ; this applies for evec-eval too, right?
+  ((set-parameters-nonsymmetric ws compute-shur-form balance))
+  :gsl-version (1 9)
+  :switch (shur-vectors)
+  :inputs (A)
+  :outputs (A eigenvalues eigenvectors)
+  :return
+  (eigenvalues
+   (cffi:foreign-slot-value (mpointer ws) 'gsl-nonsymm-ws 'n-evals))
+  :documentation			; FDL
+  "Compute eigenvalues and right eigenvectors of the n-by-n real
+  nonsymmetric matrix A. It first calls #'eigenvalues-nonsymm to
+  compute the eigenvalues, Schur form T, and Schur vectors. Then it
+  finds eigenvectors of T and backtransforms them using the Schur
+  vectors. The Schur vectors are destroyed in the process, but can be
+  saved by specifying binding shur-vectors to a vector of length n, or
+  t to have it automatically made.  The computed eigenvectors are
+  normalized to have unit magnitude. On output, the upper portion of A
+  contains the Schur form T.  If #'eigenvalues-nonsymm fails, no
+  eigenvectors are computed, and an error code is returned.
+
+  If compute-shur-form is true, the full Schur form T will be computed.
+  If it is set to nil, T will not be computed (this is
+  the default setting). Computing the full Schur form requires
+  approximately 1.5-2 times the number of flops.
+
+  If balance is true, a balancing transformation is applied to the
+  matrix prior to computing eigenvalues. This transformation is
+  designed to make the rows and columns of the matrix have comparable
+  norms, and can result in more accurate eigenvalues for matrices
+  whose entries vary widely in magnitude. See Balancing for more
+  information. Note that the balancing transformation does not
+  preserve the orthogonality of the Schur vectors, so if you wish to
+  compute the Schur vectors with you will obtain the Schur vectors of
+  the balanced matrix instead of the original matrix. The relationship
+  will be
+
+          T = Q^t D^(-1) A D Q
+
+  where Q is the matrix of Schur vectors for the balanced matrix, and D
+  is the balancing transformation. Then this function will compute
+  a matrix Z which satisfies
+
+          T = Z^(-1) A Z
+
+  with Z = D Q. Note that Z will not be orthogonal. For this reason,
+  balancing is not performed by default.")
+
 ;;;;****************************************************************************
 ;;;; Sorting Eigenvalues and Eigenvectors
 ;;;;****************************************************************************
diff --git a/init/defmfun-single.lisp b/init/defmfun-single.lisp
index b730b23d00d13e2a1cfe80ccc9e383b1c55bb7c6..7aeff7d9bbb35c0d87083c3027ac208189d53c79 100644
--- a/init/defmfun-single.lisp
+++ b/init/defmfun-single.lisp
@@ -1,6 +1,6 @@
 ;; Helpers that define a single GSL function interface
 ;; Liam Healy 2009-01-07 22:02:20EST defmfun-single.lisp
-;; Time-stamp: <2009-01-10 20:15:42EST defmfun-single.lisp>
+;; Time-stamp: <2009-01-11 10:07:06EST defmfun-single.lisp>
 ;; $Id: $
 
 (in-package :gsl)
@@ -91,7 +91,9 @@
 			 (let ((auxstart (position '&aux arglist)))
 			   ;; &aux bindings are checked
 			   (when auxstart
-			     (mapcan 'rest (subseq arglist (1+ auxstart)))))))))))
+			     (apply
+			      'append
+			      (mapcar 'rest (subseq arglist (1+ auxstart))))))))))))
 	   ,@(when documentation (list documentation))
 	   #-native
 	   ,(funcall body-maker name arglist gsl-name c-arguments key-args)
diff --git a/init/defmfun.lisp b/init/defmfun.lisp
index e3b4502805c7f9ba214d2c8fe5e4496280b2f0f7..11901874755748e910208e7f2ee6f82064c957a8 100644
--- a/init/defmfun.lisp
+++ b/init/defmfun.lisp
@@ -1,6 +1,6 @@
 ;; Macro for defining GSL functions.
 ;; Liam Healy 2008-04-16 20:49:50EDT defmfun.lisp
-;; Time-stamp: <2009-01-10 20:29:54EST defmfun.lisp>
+;; Time-stamp: <2009-01-10 21:47:36EST defmfun.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -182,18 +182,24 @@
    selects one of two GSL functions."
   (let ((optpos (position-if (lambda (s) (member s *defmfun-optk*)) arglist)))
     (if optpos
-	(let ((mandatory-arglist (subseq arglist 0 optpos))
-	      (optional-arglist (subseq arglist (1+ optpos))))
-	  `(if ,(first optional-arglist)
-	       ,(body-no-optional-arg 
-		 name
-		 (append mandatory-arglist optional-arglist)
-		 (second gsl-name)
-		 (second c-arguments)
-		 key-args)
-	       ,(body-no-optional-arg
-		 name
-		 mandatory-arglist
-		 (first gsl-name)
-		 (first c-arguments)
-		 key-args))))))
+	(with-defmfun-key-args key-args
+	  (let ((mandatory-arglist (subseq arglist 0 optpos))
+		(optional-arglist (subseq arglist (1+ optpos))))
+	    `(if ,(or (first switch) (first optional-arglist))
+		 ,(body-no-optional-arg 
+		   name
+		   (append mandatory-arglist optional-arglist)
+		   (second gsl-name)
+		   (second c-arguments)
+		   key-args)
+		 ,(body-no-optional-arg
+		   name
+		   (append mandatory-arglist
+			   (when switch
+			     (remove-if
+			      (lambda (arg)
+				(member (if (listp arg) (first arg) arg) switch))
+			      optional-arglist)))
+		   (first gsl-name)
+		   (first c-arguments)
+		   key-args)))))))
diff --git a/special-functions/coulomb.lisp b/special-functions/coulomb.lisp
index 795e88e1455e705cbc40db0fee236bcc4d24e71a..e2a85b128eac9da714dd9fe0540fa0aa9047f2b4 100644
--- a/special-functions/coulomb.lisp
+++ b/special-functions/coulomb.lisp
@@ -1,6 +1,6 @@
 ;; Coulumb functions
 ;; Liam Healy, Sat Mar 18 2006 - 23:23
-;; Time-stamp: <2008-12-29 20:49:45EST coulomb.lisp>
+;; Time-stamp: <2009-01-11 10:14:54EST coulomb.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -48,7 +48,9 @@
   overflow occurs, the condition 'overflow is signalled and scaling
   exponents are stored in the modifiable parameters exp-F, exp-G.")
 
-(defmfun coulomb-wave-F-array (L-min eta x fc-array)
+(defmfun coulomb-wave-F-array
+    (L-min eta x &optional (size-or-array *default-sf-array-size*)
+	   &aux (fc-array (vdf size-or-array)))
   "gsl_sf_coulomb_wave_F_array"
   ((L-min :double) ((1- (dim0 fc-array)) :int) (eta :double) (x :double)
    ((c-pointer fc-array) :pointer) (F-exponent :double))
@@ -59,7 +61,12 @@
   L = Lmin ... Lmin + kmax, storing the results in fc-array.
   In the case of overflow the exponent is stored in the second value returned.")
 
-(defmfun coulomb-wave-FG-array (L-min eta x fc-array gc-array)
+(defmfun coulomb-wave-FG-array
+    (L-min eta x
+	   &optional (fc-size-or-array *default-sf-array-size*)
+	   gc-size-or-array
+	   &aux (fc-array (vdf fc-size-or-array))
+	   (gc-array (vdf (or gc-size-or-array (dim0 fc-array)))))
   "gsl_sf_coulomb_wave_FG_array"
   ((L-min :double) ((1- (dim0 fc-array)) :int) (eta :double) (x :double)
    ((c-pointer fc-array) :pointer) ((c-pointer gc-array) :pointer)
@@ -72,7 +79,16 @@
   results in fc_array and gc_array.  In the case of overflow the
   exponents are stored in F_exponent and G_exponent.")
 
-(defmfun coulomb-wave-FGp-array (L-min eta x fc-array fcp-array gc-array gcp-array)
+(defmfun coulomb-wave-FGp-array
+    (L-min eta x
+	   &optional (fc-size-or-array *default-sf-array-size*)
+	   fcp-size-or-array
+	   gc-size-or-array
+	   gcp-size-or-array
+	   &aux (fc-array (vdf fc-size-or-array))
+	   (fcp-array (vdf (or fcp-size-or-array (dim0 fc-array))))
+	   (gc-array (vdf (or gc-size-or-array (dim0 fc-array))))
+	   (gcp-array (vdf (or gcp-size-or-array (dim0 fc-array)))))
   "gsl_sf_coulomb_wave_FGp_array"
   ((L-min :double) ((1- (dim0 fc-array)) :int) (eta :double) (x :double)
    ((c-pointer fc-array) :pointer) ((c-pointer fcp-array) :pointer)
@@ -90,7 +106,9 @@
   In the case of overflow the exponents are stored in F_exponent
   and G_exponent.")
 
-(defmfun coulomb-wave-sphF-array (L-min eta x fc-array)
+(defmfun coulomb-wave-sphF-array
+    (L-min eta x &optional (size-or-array *default-sf-array-size*)
+	   &aux (fc-array (vdf size-or-array)))
   "gsl_sf_coulomb_wave_sphF_array"
   ((L-min :double) ((1- (dim0 fc-array)) :int) (eta :double) (x :double)
    ((c-pointer fc-array) :pointer) (F-exponent :double))
@@ -113,14 +131,16 @@
   "The Coulomb wave function normalization constant C_L(\eta)
    for L > -1.")
 
-(defmfun coulomb-CL-array (L-min eta cl)
+(defmfun coulomb-CL-array
+    (L-min eta &optional (size-or-array *default-sf-array-size*)
+	   &aux (array (vdf size-or-array)))
   "gsl_sf_coulomb_CL_array"
-  ((L-min :double) ((1- (dim0 cl)) :int) (eta :double)
-   ((c-pointer cl) :pointer))
+  ((L-min :double) ((1- (dim0 array)) :int) (eta :double)
+   ((c-pointer array) :pointer))
+  :outputs (array)
   :documentation			; FDL
   "The Coulomb wave function normalization constant C_L(\eta)
-   for L = Lmin ... Lmin + kmax, Lmin > -1."
-  :outputs (cl))
+   for L = Lmin ... Lmin + kmax, Lmin > -1.")
 
 ;;;;****************************************************************************
 ;;;; Examples and unit test
diff --git a/special-functions/gegenbauer.lisp b/special-functions/gegenbauer.lisp
index be2c6d122a74a480f4cdadf7b3568b27ab97b57a..62fc66100aa76bbefd9a4c6d086e03905cec7b0f 100644
--- a/special-functions/gegenbauer.lisp
+++ b/special-functions/gegenbauer.lisp
@@ -1,6 +1,6 @@
 ;; Gegenbauer polynomials
 ;; Liam Healy, Fri Apr 28 2006 - 20:40
-;; Time-stamp: <2008-12-26 12:25:31EST gegenbauer.lisp>
+;; Time-stamp: <2009-01-11 09:45:00EST gegenbauer.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -27,17 +27,16 @@
   "The Gegenbauer polynomial C^{(\lambda)}_n(x)} for a specific value of n,
   lambda, x subject to \lambda > -1/2, n >= 0.")
 
-(defmfun gegenbauer-array (lambda x result)
+(defmfun gegenbauer-array
+    (lambda x &optional (size-or-array *default-sf-array-size*)
+	    &aux (array (vdf size-or-array)))
   "gsl_sf_gegenpoly_array"
-  (((1- (dim0 result)) :int)
-   (lambda :double) (x :double) ((c-pointer result) :pointer))
-  :outputs (result)
+  (((1- (dim0 array)) :int)
+   (lambda :double) (x :double) ((c-pointer array) :pointer))
+  :outputs (array)
   :documentation			; FDL
   "Compute an array of Gegenbauer polynomials C^{(\lambda)}_n(X)}
-   for n = 0, 1, 2, ..., length(result)-1}, subject to \lambda > -1/2.")
-
-;;; (defparameter vec (make-data 'vector nil 3))
-;;; (gegenbauer-array 1.0d0 3.0d0 vec)
+   for n = 0, 1, 2, ..., length(array)-1}, subject to \lambda > -1/2.")
 
 ;;;;****************************************************************************
 ;;;; Examples and unit test
@@ -48,5 +47,4 @@
   (gegenbauer-2 1.0d0 3.0d0)
   (gegenbauer-3 1.0d0 3.0d0)
   (gegenbauer 4 1.0d0 3.0d0)
-  (let ((arr (make-marray 'double-float :dimensions 4)))
-      (gegenbauer-array 1.0d0 3.0d0 arr) (cl-array arr)))
+  (cl-array (gegenbauer-array 1.0d0 3.0d0 4)))
diff --git a/special-functions/legendre.lisp b/special-functions/legendre.lisp
index 9da1514611bf193951410d114bc2f4637ee73026..15bd0b21fb450a5dab1564813a5d970c63c7c6b7 100644
--- a/special-functions/legendre.lisp
+++ b/special-functions/legendre.lisp
@@ -1,6 +1,6 @@
 ;; Legendre functions
 ;; Liam Healy, Sat Apr 29 2006 - 19:16
-;; Time-stamp: <2008-12-26 11:52:33EST legendre.lisp>
+;; Time-stamp: <2009-01-11 10:21:24EST legendre.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -35,21 +35,25 @@
   "The Legendre polynomial P_l(x) for a specific value of l,
    x subject to l >= 0, |x| <= 1.")
 
-(defmfun legendre-Pl-array (x array)
+(defmfun legendre-Pl-array
+    (x &optional (size-or-array *default-sf-array-size*)
+       &aux (array (vdf size-or-array)))
   "gsl_sf_legendre_Pl_array"
   (((1- (dim0 array)) :int) (x :double) ((c-pointer array) :pointer))
+  :outputs (array)
   :documentation			; FDL
   "Compute an array of Legendre polynomials
-  P_l(x) for l = 0, ..., length(array), |x| <= 1."
-  :outputs (array))
+  P_l(x) for l = 0, ..., length(array), |x| <= 1.")
 
-(defmfun legendre-Pl-deriv-array (x array)
+(defmfun legendre-Pl-deriv-array
+    (x &optional (size-or-array *default-sf-array-size*)
+       &aux (array (vdf size-or-array)))
   "gsl_sf_legendre_Pl_deriv_array"
   (((1- (dim0 array)) :int) (x :double) ((c-pointer array) :pointer))
+  :outputs (array)
   :documentation			; FDL
   "Compute an array of Legendre polynomials derivatives
-  dP_l(x)/dx, for l = 0, ...,  length(array), |x| <= 1."
-  :outputs (array))
+  dP_l(x)/dx, for l = 0, ...,  length(array), |x| <= 1.")
 
 (defmfun legendre-Q0 (x)
   "gsl_sf_legendre_Q0_e" ((x :double) (ret sf-result))
@@ -91,25 +95,31 @@
   "The associated Legendre polynomial
    P_l^m(x) for m >= 0, l >= m, |x| <= 1.")
 
-(defmfun legendre-Plm-array (m x array)
+(defmfun legendre-Plm-array
+    (m x &optional (size-or-array *default-sf-array-size*)
+       &aux (array (vdf size-or-array)))
   "gsl_sf_legendre_Plm_array"
   (((+ (dim0 array) m -1) :int) (m :int) (x :double)
    ((c-pointer array) :pointer))
+  :outputs (array)
   :documentation			; FDL
   "An array of Legendre polynomials
     P_l^m(x), for m >= 0, 
-    l = |m|, ..., |m|+length(array)-1} and |x| <= 1."
-  :outputs (array))
+    l = |m|, ..., |m|+length(array)-1} and |x| <= 1.")
 
-(defmfun legendre-Plm-deriv-array (m x values derivatives)
+(defmfun legendre-Plm-deriv-array
+    (m x &optional (values-size-or-array *default-sf-array-size*)
+       (derivatives-size-or-array *default-sf-array-size*)
+       &aux (values (vdf values-size-or-array))
+       (derivatives (vdf derivatives-size-or-array)))
   "gsl_sf_legendre_Plm_deriv_array"
   (((+ (dim0 values) m -1) :int) (m :int) (x :double)
    ((c-pointer values) :pointer) ((c-pointer derivatives) :pointer))
+  :outputs (values derivatives)
   :documentation			; FDL
   "An array of Legendre polynomials
     values and derivatives dP_l^m(x)/dx for m >= 0, 
-    l = |m|, ..., length(values) and |x| <= 1."
-  :outputs (values derivatives))
+    l = |m|, ..., length(values) and |x| <= 1.")
 
 (defmfun legendre-sphPlm (l m x)
   "gsl_sf_legendre_sphPlm_e" ((l :int) (m :int) (x :double) (ret sf-result))
@@ -120,25 +130,31 @@
    m >= 0, l >= m, |x| <= 1.  These routines avoid the overflows
    that occur for the standard normalization of P_l^m(x).")
 
-(defmfun legendre-sphPlm-array (m x array)
+(defmfun legendre-sphPlm-array
+    (m x &optional (size-or-array *default-sf-array-size*)
+       &aux (array (vdf size-or-array)))
   "gsl_sf_legendre_sphPlm_array"
   (((+ (dim0 array) m -1) :int) (m :int) (x :double)
    ((c-pointer array) :pointer))
+  :outputs (array)
   :documentation			; FDL
   "An array of normalized associated Legendre functions
    \sqrt(2l+1)/(4\pi) \sqrt(l-m)!/(l+m)! P_l^m(x),
-   for m >= 0, l = |m|, ..., length(array)}, |x| <= 1.0."
-  :outputs (array))
+   for m >= 0, l = |m|, ..., length(array)}, |x| <= 1.0.")
 
-(defmfun legendre-sphPlm-deriv-array (m x values derivatives)
+(defmfun legendre-sphPlm-deriv-array
+    (m x &optional (values-size-or-array *default-sf-array-size*)
+       (derivatives-size-or-array *default-sf-array-size*)
+       &aux (values (vdf values-size-or-array))
+       (derivatives (vdf derivatives-size-or-array)))
   "gsl_sf_legendre_sphPlm_deriv_array"
   (((+ (dim0 values) m -1) :int) (m :int) (x :double)
    ((c-pointer values) :pointer) ((c-pointer derivatives) :pointer))
+  :outputs (values derivatives)
   :documentation			; FDL
   "An array of normalized associated Legendre functions
    values and derivatives for m >= 0,
-   l = |m|, ..., length(array)}, |x| <= 1.0."
-  :outputs (values derivatives))
+   l = |m|, ..., length(array)}, |x| <= 1.0.")
 
 (defmfun legendre-array-size (lmax m)
   "gsl_sf_legendre_array_size" ((lmax :int) (m :int))
@@ -234,7 +250,9 @@
    \eta >= 0, l >= 0.  In the flat limit this takes the form
    L^{H3d}_l(\lambda,\eta) = j_l(\lambda\eta).")
 
-(defmfun legendre-H3d-array (lambda eta array)
+(defmfun legendre-H3d-array
+    (lambda eta &optional (size-or-array *default-sf-array-size*)
+	    &aux (array (vdf size-or-array)))
   "gsl_sf_legendre_H3d_array"
   (((1- (dim0 array)) :int) (lambda :double) (eta :double)
    ((c-pointer array) :pointer))
@@ -243,10 +261,6 @@
   "An array of radial eigenfunctions
    L^{H3d}_l(\lambda, \eta) for 0 <= l <= length(array).")
 
-;;; (defparameter hleg (make-data 'vector nil 3))
-;;; (legendre-H3d-array 1.0d0 0.5d0 hleg)
-;;; #<GSL-VECTOR #(0.9200342692589383d0 0.21694026450392123d0 0.047950660488307775d0) {C07CB51}>
-
 ;;;;****************************************************************************
 ;;;; Examples and unit test
 ;;;;****************************************************************************