diff --git a/gsll.asd b/gsll.asd
index 2143b84b96b774edb6f8ba57a8106be8dfdbf6a0..24c63f74ca28c54a4a87d4b4f4dfecbc827bbab2 100644
--- a/gsll.asd
+++ b/gsll.asd
@@ -3,7 +3,7 @@
 ; description: Definition of GSLL system 
 ; date:        
 ; author:      Liam Healy
-; modified:    Fri Apr 28 2006 - 17:24
+; modified:    Fri Apr 28 2006 - 22:35
 ;********************************************************
 ;;; $Id: $
 
@@ -66,4 +66,5 @@
 	     (:file "exponential-integrals")
 	     (:file "fermi-dirac")
 	     (:file "gamma")
+	     (:file "gegenbauer")
 	     ))))
diff --git a/init/interface.lisp b/init/interface.lisp
index d89c0dbba7884cfc01566646e4d5d5097f77f5fb..9363e9c9f43238f8a51c2f6e3990f2eb05ae5325 100644
--- a/init/interface.lisp
+++ b/init/interface.lisp
@@ -3,7 +3,7 @@
 ; description: Macros to interface GSL functions.
 ; date:        Mon Mar  6 2006 - 22:35                   
 ; author:      Liam M. Healy
-; modified:    Fri Apr 21 2006 - 09:26
+; modified:    Fri Apr 28 2006 - 22:31
 ;********************************************************
 
 (in-package :gsl)
@@ -203,46 +203,23 @@ and a scaling exponent e10, such that the value is val*10^e10."
 
 ;;; All arrays and vectors are pass with a gsl-data object that is made
 ;;; outside the function, and pieces are spliced in the function call.
-;;; Both raw arrays and GSL structures are be passed as gsl-data objects,
-;;; but their declarations are be different:
-;;;   Raw array (xyz (:double *))
-;;;   GSL struct: (xyz gsl-vector-c)
-;;; They are expanded in the foreign-funcall argument list differently:
-;;;   Raw array to two arguments (:POINTER (gsl-array xyz) :SIZE (storage-size xyz))
-;;;   GSL struct to one: (:pointer (pointer xyz))
-
-;;; (splice-arguments '((argraw (:double *)) (argstruct gsl-vector-c) (x :double)))
+;;; If a raw array is needed, use the :function argument to defun-gsl,
+;;; then use the appropriate function(s) e.g. gsl-array, dim0
+;;; on that array to map to the GSL arguments.
 (defun splice-arguments (arguments &optional mode)
   "Convert the argument declarations to a list of declarations appropriate
    for foreign-funcall.  If mode is T, a mode argument will be added to the end,
    if it is an integer, it will be put at that position."
-  (let (sizes)				; save the dimension arg if not '*
-    (flet ((splicearg (spec)
-	     ;; No accomodation for matrices (two indices) yet
-	     (if (rst-arrayp spec)
-		 (let* ((value `(first (storage-size ,(rst-symbol spec))))
-			(use-form
-			 (if (eq (rst-dim spec) '*)
-			     value
-			     (progn
-			       (push `(,(rst-dim spec) ,value) sizes)
-			       (rst-dim spec)))))
-		   (if (rst-pointer-last-p spec)
-		       `(:size ,use-form
-			 :pointer (gsl-array ,(rst-symbol spec)))
-		       `(:pointer (gsl-array ,(rst-symbol spec))
-			 :size ,use-form)))
-		 `(,(rst-type spec)
-		   ,(wrap-arg spec)))))
-      (values
-       (mapcan #'splicearg
-	       (if mode
-		   (let ((mode (if (integerp mode) mode (length arguments))))
-		     (append (subseq arguments 0 mode)
-			     '((mode sf-mode))
-			     (subseq arguments mode)))
-		   arguments))
-       sizes))))
+  (flet ((splicearg (spec)
+	   `(,(rst-type spec)
+	     ,(wrap-arg spec))))
+    (mapcan #'splicearg
+	    (if mode
+		(let ((mode (if (integerp mode) mode (length arguments))))
+		  (append (subseq arguments 0 mode)
+			  '((mode sf-mode))
+			  (subseq arguments mode)))
+		arguments))))
 
 ;;;;****************************************************************************
 ;;;; Checking results
@@ -283,8 +260,8 @@ and a scaling exponent e10, such that the value is val*10^e10."
 ;;; Warning isn't quite right for lambdas.
 (defmacro defun-gsl
     (cl-name arguments gsl-name
-	     &key documentation return mode (c-return-value :error-code)
-	     return-input check-null-pointers method after)
+     &key documentation return mode (c-return-value :error-code)
+     return-input check-null-pointers function method after)
   "Define a CL function that provides an interface to a GSL function.
    If cl-name is :lambda, make a lambda.  Arguments:
      arguments:       a list of input arguments (symbol type) to the GSL function
@@ -302,9 +279,10 @@ and a scaling exponent e10, such that the value is val*10^e10."
      method           Make output a defmethod with the value as the arglist;
                       'arguments should then include explicit mapping of all arguments
                       to GSL form.
+     function         Arguments for CL function (like :method, but make a function)
      after            Functions to call after the GSL function has been called;
                       result is discarded."
-  (let ((clargs (or method (mapcar #'rst-symbol arguments)))
+  (let ((clargs (or function method (mapcar #'rst-symbol arguments)))
 	(return-symb-type 
 	 (unless (or (eq c-return-value :return) return-input)
 	   (return-symbol-type return))))
diff --git a/linear-algebra/qr.lisp b/linear-algebra/qr.lisp
index 4284cb6566a04fc9e2c76cf8a64a6e8a468c2e79..ea78e04e2e4494126a0685abb957a64dc8bfc666 100644
--- a/linear-algebra/qr.lisp
+++ b/linear-algebra/qr.lisp
@@ -3,7 +3,7 @@
 ; description: QR decomposition                          
 ; date:        Fri Apr 28 2006 - 16:53                   
 ; author:      Liam Healy                                
-; modified:    Fri Apr 28 2006 - 17:23
+; modified:    Fri Apr 28 2006 - 21:03
 ;********************************************************
 ;;; $Id: $
 
@@ -38,22 +38,28 @@
    as used by @sc{lapack}.
 
    The algorithm used to perform the decomposition is Householder QR (Golub
-   & Van Loan, @cite{Matrix Computations}, Algorithm 5.2.1).")
+   & Van Loan, @cite{Matrix Computations}, Algorithm 5.2.1)."
+  :after ((cl-invalidate A))
+  :return-input (A tau))
 
 (defun-gsl QR-solve
     ((QR gsl-matrix-c) (tau gsl-vector-c) (b gsl-vector-c) (x gsl-vector-c))
   "gsl_linalg_QR_solve"
   :documentation "Solve the square system @math{A x = b} using the @math{QR}
    decomposition of @math{A} into (@var{QR}, @var{tau}) given by
-   @code{gsl_linalg_QR_decomp}. The least-squares solution for rectangular systems can
-   be found using QR-lssolve.")
+   QR-decomp. The least-squares solution for rectangular systems can
+   be found using QR-lssolve."
+  :after ((cl-invalidate x))
+  :return-input (x))
 
 (defun-gsl QR-svx ((QR gsl-matrix-c) (tau gsl-vector-c) (x gsl-vector-c))
   "gsl_linalg_QR_svx"
   :documentation "Solves the square system @math{A x = b} in-place using the
   @math{QR} decomposition of @math{A} into (@var{QR},@var{tau}) given by
   QR-decomp.  On input @var{x} should contain the
-  right-hand side @math{b}, which is replaced by the solution on output.")
+  right-hand side @math{b}, which is replaced by the solution on output."
+  :after ((cl-invalidate x))
+  :return-input (x))
 
 (defun-gsl QR-lssolve
     ((QR gsl-matrix-c) (tau gsl-vector-c) (b gsl-vector-c) (x gsl-vector-c)
@@ -65,7 +71,9 @@
    residual, @math{||Ax - b||}.The routine uses the @math{QR} decomposition
    of @math{A} into (@var{QR}, @var{tau}) given by
    @code{gsl_linalg_QR_decomp}.  The solution is returned in @var{x}.  The
-   residual is computed as a by-product and stored in @var{residual}.")
+   residual is computed as a by-product and stored in @var{residual}."
+  :after ((cl-invalidate x))
+  :return-input (x))
 
 (defun-gsl QR-QTvec ((QR gsl-matrix-c) (tau gsl-vector-c) (v gsl-vector-c))
   "gsl_linalg_QR_QTvec"
@@ -73,7 +81,9 @@
   (@var{QR},@var{tau}) to the vector @var{v}, storing the result @math{Q^T
   v} in @var{v}.  The matrix multiplication is carried out directly using
   the encoding of the Householder vectors without needing to form the full
-  matrix @math{Q^T}.")
+  matrix @math{Q^T}."
+  :after ((cl-invalidate v))
+  :return-input (v))
 
 (defun-gsl QR-Qvec ((QR gsl-matrix-c) (tau gsl-vector-c) (v gsl-vector-c))
   "gsl_linalg_QR_Qvec"
@@ -81,13 +91,17 @@
    (@var{QR},@var{tau}) to the vector @var{v}, storing the result @math{Q
    v} in @var{v}.  The matrix multiplication is carried out directly using
    the encoding of the Householder vectors without needing to form the full
-   matrix @math{Q}.")
+   matrix @math{Q}."
+  :after ((cl-invalidate v))
+  :return-input (v))
 
 (defun-gsl QR-Rsolve ((QR gsl-matrix-c) (b gsl-vector-c) (x gsl-vector-c))
   "gsl_linalg_QR_Rsolve"
   :documentation "Solve the triangular system @math{R x = b} for
    @var{x}. It may be useful if the product @math{b' = Q^T b} has already
-   been computed using QR-QTvec}.")
+   been computed using QR-QTvec}."
+  :after ((cl-invalidate x))
+  :return-input (x))
 
 (defun-gsl QR-Rsvx ((QR gsl-matrix-c) (x gsl-vector-c))
   "gsl_linalg_QR_Rsvx"
@@ -95,21 +109,27 @@
   in-place. On input @var{x} should contain the right-hand side @math{b}
   and is replaced by the solution on output. This function may be useful if
   the product @math{b' = Q^T b} has already been computed using
-  QR-QTvec}.")
+  QR-QTvec}."
+  :after ((cl-invalidate x))
+  :return-input (x))
 
 (defun-gsl QR-unpack
     ((QR gsl-matrix-c) (tau gsl-vector-c) (Q gsl-matrix-c) (R gsl-matrix-c))
   "gsl_linalg_QR_unpack"
   :documentation "Unpack the encoded @math{QR} decomposition
   (@var{QR},@var{tau}) into the matrices @var{Q} and @var{R}, where
-  @var{Q} is @math{M}-by-@math{M} and @var{R} is @math{M}-by-@math{N}.")
+  @var{Q} is @math{M}-by-@math{M} and @var{R} is @math{M}-by-@math{N}."
+  :after ((cl-invalidate Q R))
+  :return-input (Q R))
 
 (defun-gsl QR-QRsolve
     ((Q gsl-matrix-c) (R gsl-matrix-c) (b gsl-vector-c) (x gsl-vector-c))
   "gsl_linalg_QR_QRsolve"
   :documentation "Solves the system @math{R x = Q^T b} for @var{x}. It can
   be used when the @math{QR} decomposition of a matrix is available in
-  unpacked form as (@var{Q}, @var{R}).")
+  unpacked form as (@var{Q}, @var{R})."
+  :after ((cl-invalidate x))
+  :return-input (x))
 
 (defun-gsl QR-update
     ((Q gsl-matrix-c) (R gsl-matrix-c) (w gsl-vector-c) (v gsl-vector-c))
@@ -118,15 +138,21 @@
   decomposition (@var{Q}, @var{R}). The update is given by @math{Q'R' = Q
   R + w v^T} where the output matrices @math{Q'} and @math{R'} are also
   orthogonal and right triangular. Note that @var{w} is destroyed by the
-  update.")
+  update."
+  :after ((cl-invalidate w Q R))
+  :return-input (Q R))
 
 (defun-gsl R-solve ((R gsl-matrix-c) (b gsl-vector-c) (x gsl-vector-c))
   "gsl_linalg_R_solve"
   :documentation "Solves the triangular system @math{R x = b} for the
-   @math{N}-by-@math{N} matrix @var{R}.")
+   @math{N}-by-@math{N} matrix @var{R}."
+  :after ((cl-invalidate x))
+  :return-input (x))
 
 (defun-gsl R-svx ((R gsl-matrix-c) (x gsl-vector-c))
   "gsl_linalg_R_svx"
   :documentation "Solve the triangular system @math{R x = b} in-place. On
   input @var{x} should contain the right-hand side @math{b}, which is
-  replaced by the solution on output.")
+  replaced by the solution on output."
+  :after ((cl-invalidate x))
+  :return-input (x))
diff --git a/polynomial.lisp b/polynomial.lisp
index effd66bb9be7f9922d819bf23f1192f473b922ed..026f92fba0d911cf75376cb733e4171f69675ccb 100644
--- a/polynomial.lisp
+++ b/polynomial.lisp
@@ -3,7 +3,7 @@
 ; description: Polynomials                               
 ; date:        Tue Mar 21 2006 - 18:33                   
 ; author:      Liam M. Healy                             
-; modified:    Sat Apr 22 2006 - 15:29
+; modified:    Fri Apr 28 2006 - 22:27
 ;********************************************************
 ;;; $Id: $
 
@@ -21,8 +21,10 @@
 ;;; (setf (data vec) #(1.0d0 2.0d0 3.0d0))
 ;;; (polynomial-eval vec -1.0d0)
 ;;; 2.0d0
-(defun-gsl polynomial-eval ((coefficients (:double *)) (x :double))
+(defun-gsl polynomial-eval
+    (((gsl-array coefficients) :pointer) ((dim0 coefficients) :size) (x :double))
   "gsl_poly_eval"
+  :function (coefficients x)
   :documentation
   "Evaluate the polyonomial with coefficients at the point x."
   :return (:double)
@@ -211,9 +213,11 @@
        ,workspace))))
 
 (defun-gsl polynomial-solve-ws
-    ((coefficients (:double n)) (workspace poly-complex-workspace))
+    (((gsl-array coefficients) :pointer) ((dim0 coefficients) :size)
+     (workspace poly-complex-workspace))
   "gsl_poly_complex_solve"
-  :return ((gsl-complex (1- n))))
+  :function (coefficients workspace)
+  :return ((gsl-complex (1- (dim0 coefficients)))))
 
 #+future
 (export '(polynomial-solve))
diff --git a/sorting.lisp b/sorting.lisp
index 7bce331eb092de83e3b825bad65b80441172fcf7..bab9243af5edc7ae343b1a48be3cc6f11fc5b973 100644
--- a/sorting.lisp
+++ b/sorting.lisp
@@ -3,7 +3,7 @@
 ; description: Sorting                                   
 ; date:        Fri Apr 14 2006 - 20:20                   
 ; author:      Liam M. Healy                             
-; modified:    Wed Apr 19 2006 - 17:43
+; modified:    Fri Apr 28 2006 - 22:21
 ;********************************************************
 ;;; $Id: $
 
@@ -74,20 +74,26 @@ place.  The first element of @var{p} gives the index of the least element
 in @var{v}, and the last element of @var{p} gives the index of the
 greatest element in @var{v}.  The vector @var{v} is not changed.")
 
-(defun-gsl sort-vector-smallest ((dest (:double *)) (v gsl-vector-c))
+(defun-gsl sort-vector-smallest
+    (((gsl-array dest) :pointer) ((dim0 dest) :size) (v gsl-vector-c))
   "gsl_sort_vector_smallest"
+  :function (dest v)
   :documentation
   "Find the smallest elements of the vector @var{v} and put them into dest,
    which must be shorter than v."
   :c-return-value :void
+  :after ((cl-invalidate dest))
   :return-input (dest))
 
-(defun-gsl sort-vector-largest ((dest (:double *)) (v gsl-vector-c))
+(defun-gsl sort-vector-largest
+    (((gsl-array dest) :pointer) ((dim0 dest) :size) (v gsl-vector-c))
   "gsl_sort_vector_largest"
+  :function (dest v)
   :documentation
   "Find the largest elements of the vector @var{v} and put them into dest,
    which must be shorter than v."
   :c-return-value :void
+  :after ((cl-invalidate dest))
   :return-input (dest))
 
 ;;;;****************************************************************************
diff --git a/special-functions/bessel.lisp b/special-functions/bessel.lisp
index 95da6ec84ccd13218589be9d32a3410d6eaf563e..742af2fb84c2098057da898efdabaa8105c1896a 100644
--- a/special-functions/bessel.lisp
+++ b/special-functions/bessel.lisp
@@ -3,7 +3,7 @@
 ; description: Bessel functions                          
 ; date:        Fri Mar 17 2006 - 18:42                   
 ; author:      Liam M. Healy
-; modified:    Wed Apr 19 2006 - 17:06
+; modified:    Fri Apr 28 2006 - 22:16
 ;********************************************************
 
 (in-package :gsl)
@@ -321,20 +321,24 @@
   "The regular cylindrical Bessel function of fractional order @math{\nu}, @math{J_\nu(x)}."
   :return (sf-result))
 
-;;; (double @var{nu}, gsl_mode_t @var{mode}, size_t @var{size}, double @var{v}[])
-;;; The mode argument comes before v, need to make that an option on
-;;; the :mode argument to defun-gsl.  Need to return v as well.
-(defun-gsl bessel-sequence-Jnu ((nu :double) (v (:double * t)))
+(defun-gsl bessel-sequence-Jnu
+    ((nu :double) ((dim0 v) :size) ((gsl-array v) :pointer))
   "gsl_sf_bessel_sequence_Jnu_e"
+  :function (nu v)
   :documentation
   "The regular cylindrical Bessel function of
 fractional order @math{\nu}, @math{J_\nu(x)}, evaluated at a series of
 @math{x} values.  The array @var{v} of length @var{size} contains the
 @math{x} values.  They are assumed to be strictly ordered and positive.
 The array is over-written with the values of @math{J_\nu(x_i)}."
+  :after ((cl-invalidate v))
   :return-input (v)
   :mode 1)
 
+;;; (defparameter vec (make-data 'vector nil 3))
+;;; (setf (data vec) #(1.0d0 2.0d0 3.0d0))
+;;; (bessel-sequence-Jnu 0.5d0 vec)
+
 ;;;;****************************************************************************
 ;;;; Irregular Bessel Function - Fractional Order
 ;;;;****************************************************************************
diff --git a/special-functions/gegenbauer.lisp b/special-functions/gegenbauer.lisp
new file mode 100644
index 0000000000000000000000000000000000000000..7344e60031d9a98cc61e88c9bcde8a28fe58218d
--- /dev/null
+++ b/special-functions/gegenbauer.lisp
@@ -0,0 +1,62 @@
+;********************************************************
+; file:        gegenbauer.lisp                           
+; description: Gegenbauer polynomials                    
+; date:        Fri Apr 28 2006 - 20:40                   
+; author:      Liam M. Healy                             
+; modified:    Fri Apr 28 2006 - 22:44
+;********************************************************
+;;; $Id: $
+
+(in-package :gsl)
+
+(defun-gsl gegenbauer-1 ((lambda :double) (x :double))
+  "gsl_sf_gegenpoly_1_e"
+  :return (sf-result)
+  :documentation "The Gegenbauer polynomial @math{C^@{(\lambda)@}_1(x)}.")
+
+(defun-gsl gegenbauer-2 ((lambda :double) (x :double))
+  "gsl_sf_gegenpoly_2_e"
+  :return (sf-result)
+  :documentation "The Gegenbauer polynomial @math{C^@{(\lambda)@}_2(x)}.")
+
+(defun-gsl gegenbauer-3 ((lambda :double) (x :double))
+  "gsl_sf_gegenpoly_3_e"
+  :return (sf-result)
+  :documentation "The Gegenbauer polynomial @math{C^@{(\lambda)@}_3(x)}.")
+
+(defun-gsl gegenbauer ((n :int) (lambda :double) (x :double))
+  "gsl_sf_gegenpoly_n_e"
+  :return (sf-result)
+  :documentation "The Gegenbauer polynomial
+  @math{C^@{(\lambda)@}_n(x)} for a specific value of @var{n},
+  @var{lambda}, @var{x} subject to @math{\lambda > -1/2},
+  @math{n >= 0}.")
+
+(defun-gsl gegenbauer-array
+    (((dim0 result) :int)
+     (lambda :double) (x :double) ((gsl-array result) :pointer))
+  "gsl_sf_gegenpoly_array"
+  :function (lambda x result)
+  :documentation "Compute an array of Gegenbauer polynomials
+  @math{C^@{(\lambda)@}_n(x)} for @math{n = 0, 1, 2, \dots, nmax}, subject
+  to @math{\lambda > -1/2}, @math{nmax >= 0}."
+  :after ((cl-invalidate result))
+  :return-input (result))
+
+;;; (defparameter vec (make-data 'vector nil 3))
+;;; (gegenbauer-array 1.0d0 3.0d0 vec)
+
+;;;;****************************************************************************
+;;;; Examples and unit test
+;;;;****************************************************************************
+
+(lisp-unit:define-test gegenbauer
+  (lisp-unit:assert-first-fp-equal
+   "0.600000000000d+01"
+   (gegenbauer-1 1.0d0 3.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.350000000000d+02"
+   (gegenbauer-2 1.0d0 3.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.204000000000d+03"
+   (gegenbauer-3 1.0d0 3.0d0)))