diff --git a/gsll.asd b/gsll.asd
index ec532806ea0da02421e947d45aef3b345afbfcbe..6c31eb69505aa9608555d0b030bbed6a24851b43 100644
--- a/gsll.asd
+++ b/gsll.asd
@@ -3,7 +3,7 @@
 ; description: Definition of GSLL system 
 ; date:        
 ; author:      Liam Healy
-; modified:    Mon May  1 2006 - 15:29
+; modified:    Mon May  1 2006 - 22:33
 ;********************************************************
 ;;; $Id: $
 
@@ -56,6 +56,8 @@
 	     (:file "legendre")
 	     (:file "logarithm")
 	     (:file "power")
+	     (:file "psi")
+	     (:file "synchrotron")
 	     ))
    (:module data
 	    :depends-on (init)
diff --git a/init/interface.lisp b/init/interface.lisp
index f9814f797ec61efb6d8faf099d1a0c544cd3f0a2..f6c7982f7b5c3beab7c8914486be90ef48a043a8 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:    Sun Apr 30 2006 - 22:35
+; modified:    Mon May  1 2006 - 22:55
 ;********************************************************
 
 (in-package :gsl)
@@ -322,7 +322,7 @@ and a scaling exponent e10, such that the value is val*10^e10."
 	    ,@(case c-return-value
 		    (:void '((declare (ignore creturn))))
 		    (:error-code
-		     (if method
+		     (if (or function method)
 			 `((check-gsl-status creturn `(,',cl-name))) ; need args
 			 `((check-gsl-status creturn `(,',cl-name ,,@clargs))))))
 	    ,@(check-null-pointers check-null-pointers)
diff --git a/special-functions/bessel.lisp b/special-functions/bessel.lisp
index 742af2fb84c2098057da898efdabaa8105c1896a..88fab4709130a76cb591a7296243bf2c0957712f 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:    Fri Apr 28 2006 - 22:16
+; modified:    Mon May  1 2006 - 23:46
 ;********************************************************
 
 (in-package :gsl)
@@ -30,11 +30,23 @@
   "The regular cylindrical Bessel function of order @var{n}, @math{J_n(x)}."
   :return (sf-result))
 
-(defun-gsl cylindrical-bessel-Jn-array ((nmin :int) (nmax :int) (x :double))
+(defun-gsl cylindrical-bessel-Jn-array
+    ((nmin :int) ((+ nmin (dim0 array)) :int) (x :double)
+     ((gsl-array array) :pointer))
   "gsl_sf_bessel_Jn_array"
+  :function (x array &optional (nmin 0))
   :documentation
-  "The values of the regular cylindrical Bessel functions @math{J_n(x)} for @math{n} from @var{nmin} to @var{nmax} inclusive.  The values are computed using recurrence relations for efficiency, and therefore may differ slightly from the exact values."
-  :return ((:double (1+ (- nmax nmin)))))
+  "The values of the regular cylindrical Bessel functions @math{J_n(x)}
+   for @math{n} from @var{nmin} to @var{nmax} inclusive.
+   The values are computed using recurrence relations for efficiency,
+   and therefore may differ slightly from the exact values."
+  :invalidate (array)
+  :return-input (array))
+
+;;; (defparameter besarr (make-data 'vector nil 4))
+;;; (cylindrical-bessel-Jn-array 2.0d0 2 besarr)
+;;; #<GSL-VECTOR #(0.35283402861563773d0 0.12894324947440206d0
+;;;               0.033995719807568436d0 0.007039629755871685d0) {C7D6B51}>
 
 ;;;;****************************************************************************
 ;;;; Irregular Cylindrical Bessel Functions
@@ -58,11 +70,18 @@
   "The irregular cylindrical Bessel function of order @var{n}, @math{Y_n(x)}."
   :return (sf-result))
 
-(defun-gsl cylindrical-bessel-Yn-array ((nmin :int) (nmax :int) (x :double))
+(defun-gsl cylindrical-bessel-Yn-array
+    ((nmin :int) ((+ nmin (dim0 array)) :int) (x :double)
+     ((gsl-array array) :pointer))
   "gsl_sf_bessel_Yn_array"
+  :function (x array &optional (nmin 0))
   :documentation
-  "The values of the irregular cylindrical Bessel functions @math{Y_n(x)} for @math{n} from @var{nmin} to @var{nmax} inclusive.  The values are computed using recurrence relations for efficiency, and therefore may differ slightly from the exact values."
-  :return ((:double (1+ (- nmax nmin)))))
+  "The values of the irregular cylindrical Bessel functions
+  @math{Y_n(x)} for @math{n} from @var{nmin} to @var{nmax} inclusive.
+  The values are computed using recurrence relations for efficiency, and
+  therefore may differ slightly from the exact values."
+  :invalidate (array)
+  :return-input (array))
 
 ;;;;****************************************************************************
 ;;;; Regular Modified Cylindrical Bessel Functions
@@ -163,12 +182,19 @@
   "The scaled irregular modified cylindrical Bessel function of order @var{n}, @math{\exp(-|x|) K_n(x)}}."
   :return (sf-result))
 
-(defun-gsl cylindrical-bessel-Kn-scaled-array ((nmin :int) (nmax :int) (x :double))
+(defun-gsl cylindrical-bessel-Kn-scaled-array
+    ((nmin :int) ((+ nmin (dim0 array)) :int) (x :double)
+     ((gsl-array array) :pointer))
   "gsl_sf_bessel_Kn_scaled_array"
+  :function (x array &optional (nmin 0))
   :documentation
   ;; Bug in original documentation?  doesn't say "modified"
-  "The values of the scaled irregular modified cylindrical Bessel functions @math{K_n(x)} for @math{n} from @var{nmin} to @var{nmax} inclusive.  The values are computed using recurrence relations for efficiency, and therefore may differ slightly from the exact values."
-  :return ((:double (1+ (- nmax nmin)))))
+  "The values of the scaled irregular modified cylindrical Bessel
+  functions @math{K_n(x)} for @math{n} from @var{nmin} to @var{nmax}
+  inclusive.  The values are computed using recurrence relations for
+  efficiency, and therefore may differ slightly from the exact values."
+  :invalidate (array)
+  :return-input (array))
 
 ;;;;****************************************************************************
 ;;;; Regular Spherical Bessel Functions
@@ -177,37 +203,61 @@
 (defun-gsl spherical-bessel-j0 ((x :double))
   "gsl_sf_bessel_j0_e"
   :documentation
-  "The regular spherical Bessel function of zeroth order, @math{j_0(x) = \sin(x)/x}."
+  "The regular spherical Bessel function of zeroth order, @math{j_0(x)
+  = \sin(x)/x}."
   :return (sf-result))
 
 (defun-gsl spherical-bessel-j1 ((x :double))
   "gsl_sf_bessel_j1_e"
   :documentation
-  "The regular spherical Bessel function of first order, @math{j_1(x) = (\sin(x)/x - \cos(x))/x}."
+  "The regular spherical Bessel function of first order, @math{j_1(x)
+   = (\sin(x)/x - \cos(x))/x}."
   :return (sf-result))
 
 (defun-gsl spherical-bessel-j2 ((x :double))
   "gsl_sf_bessel_j2_e"
   :documentation
-  "The regular spherical Bessel function of second order, @math{j_2(x) = ((3/x^2 - 1)\sin(x) - 3\cos(x)/x)/x}."
+  "The regular spherical Bessel function of second order, @math{j_2(x)
+   = ((3/x^2 - 1)\sin(x) - 3\cos(x)/x)/x}."
   :return (sf-result))
 
 (defun-gsl spherical-bessel-jl ((l :int) (x :double))
   "gsl_sf_bessel_jl_e"
-  :documentation "The regular spherical Bessel function of order @var{l}, @math{j_l(x)}, for @c{$l \geq 0$} @math{l >= 0} and @c{$x \geq 0$} @math{x >= 0}."
+  :documentation "The regular spherical Bessel function of order
+  @var{l}, @math{j_l(x)}, for @math{l >= 0} and @math{x >= 0}."
   :return (sf-result))
 
-(defun-gsl spherical-bessel-jl-array ((lmax :int) (x :double))
+(defun-gsl spherical-bessel-jl-array
+    (((dim0 array) :int) (x :double)
+     ((gsl-array array) :pointer))
   "gsl_sf_bessel_jl_array"
-  :documentation
-  "The regular spherical Bessel functions @math{j_l(x)} for @math{l} from 0 to @var{lmax} inclusive  for @c{$lmax \geq 0$} @math{lmax >= 0} and @c{$x \geq 0$} @math{x >= 0}, storing the results in the array @var{result_array}. The values are computed using recurrence relations for efficiency, and therefore may differ slightly from the exact values."
-  :return ((:double (1+ lmax))))
-
-(defun-gsl spherical-bessel-jl-steed-array ((lmax :int) (x :double))
+  :function (x array)
+  :documentation
+  "The regular spherical Bessel functions @math{j_l(x)} for @math{l}
+  from 0 to @var{lmax} inclusive for @math{lmax >= 0}
+  and @math{x >= 0}, storing the results in the array
+  @var{result_array}. The values are computed using recurrence relations
+  for efficiency, and therefore may differ slightly from the exact
+  values."
+  :invalidate (array)
+  :return-input (array))
+
+(defun-gsl spherical-bessel-jl-steed-array
+    (((dim0 array) :int) (x :double)
+     ((gsl-array array) :pointer))
   "gsl_sf_bessel_jl_steed_array"
-  :documentation
-  "Uses Steed's method to compute the values of the regular spherical Bessel functions @math{j_l(x)} for @math{l} from 0 to @var{lmax} inclusive for @c{$lmax \geq 0$} @math{lmax >= 0} and @c{$x \geq 0$} @math{x >= 0}, storing the results in the array @var{result_array}. The Steed/Barnett algorithm is described in @cite{Comp. Phys. Comm.} 21, 297 (1981).  Steed's method is more stable than the recurrence used in the other functions but is also slower."
-  :return ((:double (1+ lmax))))
+  :function (x array)
+  :documentation
+  "Uses Steed's method to compute the values of the regular spherical
+  Bessel functions @math{j_l(x)} for @math{l} from 0 to @var{lmax}
+  inclusive for @math{lmax >= 0} and @math{x >= 0},
+  storing the results in the array
+  @var{result_array}. The Steed/Barnett algorithm is described in
+  @cite{Comp. Phys. Comm.} 21, 297 (1981).  Steed's method is more
+  stable than the recurrence used in the other functions but is also
+  slower."
+  :invalidate (array)
+  :return-input (array))
 
 ;;;;****************************************************************************
 ;;;; Irregular Spherical Bessel Functions
@@ -216,32 +266,44 @@
 (defun-gsl spherical-bessel-y0 ((x :double))
   "gsl_sf_bessel_y0_e"
   :documentation
-  "The irregular spherical Bessel function of zeroth order, @math{y_0(x) = -\cos(x)/x}."
+  "The irregular spherical Bessel function of zeroth order,
+  @math{y_0(x) = -\cos(x)/x}."
   :return (sf-result))
 
 (defun-gsl spherical-bessel-y1 ((x :double))
   "gsl_sf_bessel_y1_e"
   :documentation
-  "The irregular spherical Bessel function of first order, @math{y_1(x) = -(\cos(x)/x + \sin(x))/x}."
+  "The irregular spherical Bessel function of first order,
+  @math{y_1(x) = -(\cos(x)/x + \sin(x))/x}."
   :return (sf-result))
 
 (defun-gsl spherical-bessel-y2 ((x :double))
   "gsl_sf_bessel_y2_e"
   :documentation
-  "The irregular spherical Bessel function of second order, @math{y_2(x) = (-3/x^3 + 1/x)\cos(x) - (3/x^2)\sin(x)}."
+  "The irregular spherical Bessel function of second order,
+  @math{y_2(x) = (-3/x^3 + 1/x)\cos(x) - (3/x^2)\sin(x)}."
   :return (sf-result))
 
 (defun-gsl spherical-bessel-yl ((l :int) (x :double))
   "gsl_sf_bessel_yl_e"
   :documentation
-  "The irregular spherical Bessel function of order @var{l}, @math{y_l(x)}, for @c{$l \geq 0$} @math{l >= 0}."
+  "The irregular spherical Bessel function of order @var{l},
+  @math{y_l(x)}, for @math{l >= 0}."
   :return (sf-result))
 
-(defun-gsl spherical-bessel-yl-array ((lmax :int) (x :double))
+(defun-gsl spherical-bessel-yl-array
+    (((dim0 array) :int) (x :double)
+     ((gsl-array array) :pointer))
   "gsl_sf_bessel_yl_array"
+  :function (x array)
   :documentation
-  "The irregular spherical Bessel functions @math{y_l(x)} for @math{l} from 0 to @var{lmax} inclusive  for @c{$lmax \geq 0$} @math{lmax >= 0}, storing the results in the array @var{result_array}. The values are computed using recurrence relations for efficiency, and therefore may differ slightly from the exact values."
-  :return ((:double (1+ lmax))))
+  "The irregular spherical Bessel functions @math{y_l(x)} for @math{l}
+  from 0 to @var{lmax} inclusive for @math{lmax >= 0},
+  storing the results in the array @var{result_array}. The values are
+  computed using recurrence relations for efficiency, and therefore may
+  differ slightly from the exact values."
+  :invalidate (array)
+  :return-input (array))
 
 ;;;;****************************************************************************
 ;;;; Regular Modified Spherical Bessel Functions
@@ -250,7 +312,8 @@
 (defun-gsl spherical-bessel-i0-scaled ((x :double))
   "gsl_sf_bessel_i0_scaled_e"
   :documentation
-  "The scaled regular modified spherical Bessel function of zeroth order, @math{\exp(-|x|) i_0(x)}."
+  "The scaled regular modified spherical Bessel function of zeroth
+  order, @math{\exp(-|x|) i_0(x)}."
   :return (sf-result))
 
 (defun-gsl spherical-bessel-i1-scaled ((x :double))
@@ -271,11 +334,20 @@
   "The scaled regular modified spherical Bessel function of order @var{l}, @math{ \exp(-|x|) i_l(x) }."
   :return (sf-result))
 
-(defun-gsl spherical-bessel-il-scaled-array ((lmax :int) (x :double))
+(defun-gsl spherical-bessel-il-scaled-array
+    (((dim0 array) :int) (x :double)
+     ((gsl-array array) :pointer))
   "gsl_sf_bessel_il_scaled_array"
+  :function (x array)
   :documentation
-  "The values of the scaled regular modified cylindrical Bessel functions @math{\exp(-|x|) i_l(x)} for @math{l} from 0 to @var{lmax} inclusive for @c{$lmax \geq 0$} @math{lmax >= 0}, storing the results in the array @var{result_array}.  The values are computed using recurrence relations for efficiency, and therefore may differ slightly from the exact values."
-  :return ((:double (1+ lmax))))
+  "The values of the scaled regular modified cylindrical Bessel
+  functions @math{\exp(-|x|) i_l(x)} for @math{l} from 0 to @var{lmax}
+  inclusive for @math{lmax >= 0}, storing the results
+  in the array @var{result_array}.  The values are computed using
+  recurrence relations for efficiency, and therefore may differ slightly
+  from the exact values."
+  :invalidate (array)
+  :return-input (array))
 
 ;;;;****************************************************************************
 ;;;; Irregular Modified Spherical Bessel Functions
@@ -284,7 +356,8 @@
 (defun-gsl spherical-bessel-k0-scaled ((x :double))
   "gsl_sf_bessel_k0_scaled_e"
   :documentation
-  "The scaled irregular modified spherical Bessel function of zeroth order, @math{\exp(x) k_0(x)}, for @math{x>0}."
+  "The scaled irregular modified spherical Bessel function of zeroth
+  order, @math{\exp(x) k_0(x)}, for @math{x>0}."
   :return (sf-result))
 
 (defun-gsl spherical-bessel-k1-scaled ((x :double))
@@ -305,11 +378,20 @@
   "The scaled irregular modified spherical Bessel function of order @var{l}, @math{\exp(x) k_l(x)}, for @math{x>0}."
   :return (sf-result))
 
-(defun-gsl spherical-bessel-kl-scaled-array ((lmax :int) (x :double))
+(defun-gsl spherical-bessel-kl-scaled-array
+    (((dim0 array) :int) (x :double)
+     ((gsl-array array) :pointer))
   "gsl_sf_bessel_kl_scaled_array"
+  :function (x array)
   :documentation
-  "The values of the scaled irregular modified spherical Bessel functions @math{\exp(x) k_l(x)} for @math{l} from 0 to @var{lmax} inclusive for @c{$lmax \geq 0$} @math{lmax >= 0} and @math{x>0}, storing the results in the array @var{result_array}.  The values are computed using recurrence relations for efficiency, and therefore may differ slightly from the exact values."
-  :return ((:double (1+ lmax))))
+  "The values of the scaled irregular modified spherical Bessel
+  functions @math{\exp(x) k_l(x)} for @math{l} from 0 to @var{lmax}
+  inclusive for @math{lmax >= 0} and @math{x>0},
+  storing the results in the array @var{result_array}.  The values are
+  computed using recurrence relations for efficiency, and therefore may
+  differ slightly from the exact values."
+  :invalidate (array)
+  :return-input (array))
 
 ;;;;****************************************************************************
 ;;;; Regular Bessel Function - Fractional Order
@@ -318,7 +400,8 @@
 (defun-gsl bessel-Jnu ((nu :double) (x :double))
   "gsl_sf_bessel_Jnu_e"
   :documentation
-  "The regular cylindrical Bessel function of fractional order @math{\nu}, @math{J_\nu(x)}."
+  "The regular cylindrical Bessel function of fractional order
+  @math{\nu}, @math{J_\nu(x)}."
   :return (sf-result))
 
 (defun-gsl bessel-sequence-Jnu
@@ -327,10 +410,10 @@
   :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)}."
+  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)
@@ -346,7 +429,8 @@ The array is over-written with the values of @math{J_\nu(x_i)}."
 (defun-gsl bessel-Ynu ((nu :double) (x :double))
   "gsl_sf_bessel_Ynu_e"
   :documentation
-  "The irregular cylindrical Bessel function of fractional order @math{\nu}, @math{Y_\nu(x)}."
+  "The irregular cylindrical Bessel function of fractional order
+  @math{\nu}, @math{Y_\nu(x)}."
   :return (sf-result))
 
 ;;;;****************************************************************************
@@ -356,13 +440,15 @@ The array is over-written with the values of @math{J_\nu(x_i)}."
 (defun-gsl bessel-Inu ((nu :double) (x :double))
   "gsl_sf_bessel_Inu_e"
   :documentation
-  "The regular modified Bessel function of fractional order @math{\nu}, @math{I_\nu(x)} for @math{x>0}, @math{\nu>0}."
+  "The regular modified Bessel function of fractional order
+  @math{\nu}, @math{I_\nu(x)} for @math{x>0}, @math{\nu>0}."
   :return (sf-result))
 
 (defun-gsl bessel-Inu-scaled ((nu :double) (x :double))
   "gsl_sf_bessel_Inu_scaled_e"
   :documentation
-  "The scaled regular modified Bessel function of fractional order @math{\nu}, @math{\exp(-|x|)I_\nu(x)} for @math{x>0}, @math{\nu>0}."
+  "The scaled regular modified Bessel function of fractional order
+  @math{\nu}, @math{\exp(-|x|)I_\nu(x)} for @math{x>0}, @math{\nu>0}."
   :return (sf-result))
 
 ;;;;****************************************************************************
@@ -394,19 +480,141 @@ The array is over-written with the values of @math{J_\nu(x_i)}."
 (defun-gsl bessel-zero-J0 ((s :int))
   "gsl_sf_bessel_zero_J0_e"
   :documentation
-  "The location of the @var{s}-th positive zero of the Bessel function @math{J_0(x)}."
+  "The location of the @var{s}-th positive zero of the Bessel function
+  @math{J_0(x)}."
   :return (sf-result))
 
 (defun-gsl bessel-zero-J1 ((s :int))
   "gsl_sf_bessel_zero_J1_e"
   :documentation
-  "The location of the @var{s}-th positive zero of the Bessel function @math{J_1(x)}."
+  "The location of the @var{s}-th positive zero of the Bessel function
+  @math{J_1(x)}."
   :return (sf-result))
 
-(defun-gsl bessel-zero-Jnu ((s :int))
+(defun-gsl bessel-zero-Jnu ((nu :double) (s :int))
   "gsl_sf_bessel_zero_Jnu_e"
   :documentation
-  "These routines compute the location of the @var{s}-th positive zero of
-the Bessel function @math{J_\nu(x)}.  The current implementation does not
-support negative values of @var{nu}."
+  "These routines compute the location of the @var{s}-th positive zero
+  of the Bessel function @math{J_\nu(x)}.  The current implementation
+  does not support negative values of @var{nu}."
   :return (sf-result))
+
+;;;;****************************************************************************
+;;;; Examples and unit test
+;;;;****************************************************************************
+
+(lisp-unit:define-test bessel
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL
+   "-0.397149809864d+00"
+   (CYLINDRICAL-BESSEL-J0 4.0d0))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL
+   "-0.660433280235d-01"
+   (CYLINDRICAL-BESSEL-J1 4.0d0))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL
+   "0.364128145852d+00"
+   (CYLINDRICAL-BESSEL-JN 2 4.0d0))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL
+   "-0.169407393251d-01"
+   (CYLINDRICAL-BESSEL-Y0 4.0d0))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL
+   "0.397925710557d+00"
+   (CYLINDRICAL-BESSEL-Y1 4.0d0))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL
+   "-0.182022115953d+00"
+   (CYLINDRICAL-BESSEL-YN 3 4.0d0))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL
+   "0.113019219521d+02"
+   (CYLINDRICAL-BESSEL-I0 4.0d0))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL
+   "0.975946515370d+01"
+   (CYLINDRICAL-BESSEL-I1 4.0d0))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL
+   "0.333727577842d+01"
+   (CYLINDRICAL-BESSEL-IN 3 4.0d0))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL
+   "0.207001921224d+00"
+   (CYLINDRICAL-BESSEL-I0-SCALED 4.0d0))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL
+   "0.178750839502d+00"
+   (CYLINDRICAL-BESSEL-I1-SCALED 4.0d0))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL
+   "0.611243380297d-01"
+   (CYLINDRICAL-BESSEL-IN-SCALED 3 4.0d0))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL
+   "-0.189200623827d+00"
+   (SPHERICAL-BESSEL-J0 4.0d0))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL
+   "0.116110749259d+00"
+   (SPHERICAL-BESSEL-J1 4.0d0))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL
+   "0.276283685771d+00"
+   (SPHERICAL-BESSEL-J2 4.0d0))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL
+   "0.229243857955d+00"
+   (SPHERICAL-BESSEL-JL 3 4.0d0))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL
+   "0.163410905216d+00"
+   (SPHERICAL-BESSEL-Y0 4.0d0))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL
+   "0.230053350131d+00"
+   (SPHERICAL-BESSEL-Y1 4.0d0))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL
+   "0.912910738232d-02"
+   (SPHERICAL-BESSEL-Y2 4.0d0))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL
+   "0.912910738232d-02"
+   (SPHERICAL-BESSEL-YL 2 4.0d0))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL
+   "0.124958067172d+00"
+   (SPHERICAL-BESSEL-I0-SCALED 4.0d0))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL
+   "0.938024160356d-01"
+   (SPHERICAL-BESSEL-I1-SCALED 4.0d0))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL
+   "0.546062551448d-01"
+   (SPHERICAL-BESSEL-I2-SCALED 4.0d0))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL
+   "0.255445971046d-01"
+   (SPHERICAL-BESSEL-IL-SCALED 3 4.0d0))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL
+   "0.392699081699d+00"
+   (SPHERICAL-BESSEL-K0-SCALED 4.0d0))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL
+   "0.490873852123d+00"
+   (SPHERICAL-BESSEL-K1-SCALED 4.0d0))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL
+   "0.760854470791d+00"
+   (SPHERICAL-BESSEL-K2-SCALED 4.0d0))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL
+   "0.255445971046d-01"
+   (SPHERICAL-BESSEL-KL-SCALED 3 4.0d0))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL
+   "0.430171473876d+00"
+   (BESSEL-JNU 3.0d0 4.0d0))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL
+   "-0.182022115953d+00"
+   (BESSEL-YNU 3.0d0 4.0d0))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL
+   "0.333727577842d+01"
+   (BESSEL-INU 3.0d0 4.0d0))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL
+   "0.611243380297d-01"
+   (BESSEL-INU-SCALED 3.0d0 4.0d0))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL
+   "0.298849244168d-01"
+   (BESSEL-KNU 3.0d0 4.0d0))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL
+   "-0.351040112585d+01"
+   (BESSEL-LNKNU 3.0d0 4.0d0))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL
+   "0.163166158704d+01"
+   (BESSEL-KNU-SCALED 3.0d0 4.0d0))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL
+   "0.149309177085d+02"
+   (BESSEL-ZERO-J0 5))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL
+   "0.164706300509d+02"
+   (BESSEL-ZERO-J1 5))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL
+   "0.179598194950d+02"
+   (BESSEL-ZERO-JNU 2.0d0 5)))
diff --git a/special-functions/psi.lisp b/special-functions/psi.lisp
new file mode 100644
index 0000000000000000000000000000000000000000..cef02e7e011e8cb164d055d0610a33f7659f5651
--- /dev/null
+++ b/special-functions/psi.lisp
@@ -0,0 +1,75 @@
+;********************************************************
+; file:        psi.lisp                                  
+; description: Psi (digamma) functions                   
+; date:        Mon May  1 2006 - 22:11                   
+; author:      Liam M. Healy                             
+; modified:    Mon May  1 2006 - 22:25
+;********************************************************
+;;; $Id: $
+
+(in-package :gsl)
+
+;;;;****************************************************************************
+;;;; Digamma Function
+;;;;****************************************************************************
+
+(defgeneric psi (x)
+  (:documentation "The psi, or digamma, function."))
+
+(defun-gsl psi ((n :int))
+  "gsl_sf_psi_int_e"
+  :method ((n fixnum))
+  :documentation "Domain: n integer, n > 0."
+  :return (sf-result))
+
+(defun-gsl psi ((x :double))
+  "gsl_sf_psi_e"
+  :method ((x double-float))
+  :documentation "Domain: x /= 0.0, -1.0, -2.0, ..."
+  :return (sf-result))
+
+(defun-gsl psi-1piy ((x :double))
+  "gsl_sf_psi_1piy_e"
+  :documentation "The real part of the digamma function
+  on the line @math{1+i y}, @math{\Re[\psi(1 + i y)]}."
+  :return (sf-result))
+
+;;;;****************************************************************************
+;;;; Trigamma Function
+;;;;****************************************************************************
+
+(defgeneric psi-1 (x)
+  (:documentation "The Trigamma function."))
+
+(defun-gsl psi-1 ((n :int))
+  "gsl_sf_psi_1_int_e"
+  :method ((n fixnum))
+  :documentation "Domain: n integer, n > 0."
+  :return (sf-result))
+
+(defun-gsl psi-1 ((x :double))
+  "gsl_sf_psi_1_e"
+  :method ((x double-float))
+  :documentation "Domain: x /= 0.0, -1.0, -2.0, ..."
+  :return (sf-result))
+
+;;;;****************************************************************************
+;;;; Polygamma
+;;;;****************************************************************************
+
+(defun-gsl psi-n ((m :int) (x :double))
+  "gsl_sf_psi_n_e"
+  :documentation "The polygamma function @math{\psi^@{(m)@}(x)} for
+   @math{m >= 0}, @math{x > 0}."
+  :return (sf-result))
+
+;;;;****************************************************************************
+;;;; Examples and unit test
+;;;;****************************************************************************
+
+(lisp-unit:define-test psi
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL "0.125611766843d+01" (PSI 4))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL "0.125611766843d+01" (PSI 4.0d0))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL "0.283822955737d+00" (PSI-1 4))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL "0.283822955737d+00" (PSI-1 4.0d0))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL "-0.800397322451d-01" (PSI-N 2 4.0d0)))
diff --git a/special-functions/synchrotron.lisp b/special-functions/synchrotron.lisp
new file mode 100644
index 0000000000000000000000000000000000000000..5c6ee1442614f7d6db4c7c992534cadc5ada5640
--- /dev/null
+++ b/special-functions/synchrotron.lisp
@@ -0,0 +1,30 @@
+;********************************************************
+; file:        synchrotron.lisp                                
+; description: Synchrotron functions
+; date:        Mon May  1 2006 - 22:29
+; author:      Liam M. Healy                             
+; modified:    Mon May  1 2006 - 22:30
+;********************************************************
+;;; $Id: $
+
+(in-package :gsl)
+
+(defun-gsl synchrotron-1 ((x :double))
+  "gsl_sf_synchrotron_1_e"
+  :documentation "The first synchrotron function 
+  @math{x \int_x^\infty dt K_@{5/3@}(t)} for @math{x >= 0}."
+  :return (sf-result))
+
+(defun-gsl synchrotron-2 ((x :double))
+  "gsl_sf_synchrotron_2_e"
+  :documentation "The second synchrotron function 
+  @math{x K_@{2/3@}(x)} for @math{x >= 0}."
+  :return (sf-result))
+
+;;;;****************************************************************************
+;;;; Examples and unit test
+;;;;****************************************************************************
+
+(lisp-unit:define-test synchrotron
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL "0.528273966979d-01" (SYNCHROTRON-1 4.0d0))
+  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL "0.469232058261d-01" (SYNCHROTRON-2 4.0d0)))