diff --git a/gsll.asd b/gsll.asd
index d04918298fd66ea23f367c028d4a46a7689d32f0..a6c814c02aee82828b0aeed87c3b6317c4315fd7 100644
--- a/gsll.asd
+++ b/gsll.asd
@@ -3,7 +3,7 @@
 ; description: Definition of GSLL system 
 ; date:        
 ; author:      Liam Healy
-; modified:    Fri Mar 17 2006 - 18:52
+; modified:    Sun Mar 19 2006 - 00:41
 ;********************************************************
 ;;; $Id: $
 
@@ -27,5 +27,8 @@
 	    :depends-on (init interface)
 	    :components
 	    ((:file "airy")
-	     (:file "bessel")
-	     (:file "dilogarithm")))))
+	     (:file "bessel")		; one left 
+	     (:file "clausen")
+	     (:file "coulomb")
+	     (:file "dilogarithm")
+	     ))))
diff --git a/special-functions/bessel.lisp b/special-functions/bessel.lisp
index d3bacd80ea379ad38ecc9154c3d8d6857b5f3cb4..4a270d6256af1517a849713288c982229bfdc8b8 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:    Sat Mar 18 2006 - 00:15
+; modified:    Sat Mar 18 2006 - 22:57
 ;********************************************************
 
 (in-package :gsl)
@@ -324,15 +324,46 @@
   "The regular cylindrical Bessel function of fractional order @math{\nu}, @math{J_\nu(x)}."
   :return (sf-result))
 
-;;; Needs work:
-;;; has mode, same array for input and output
 ;;; (double @var{nu}, gsl_mode_t @var{mode}, size_t @var{size}, double @var{v}[])
-;;; "gsl_sf_bessel_sequence_Jnu_e"
-"The regular cylindrical Bessel function of
+#+development
+(defun-sf bessel-sequence-Jnu ((nu :double) (v :pointer))
+  "gsl_sf_bessel_sequence_Jnu_e"
+  :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)}."
+  :mode t
+  :return (v :pointer))
+
+#+development
+(DEFUNX BESSEL-SEQUENCE-JNU
+    (NU V &OPTIONAL (MODE :DOUBLE-PREC))
+  "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)}."
+  (LET ((STATUS
+	 (FOREIGN-FUNCALL
+	  "gsl_sf_bessel_sequence_Jnu_e"
+	  :DOUBLE NU
+	  SF-MODE MODE
+	  :POINTER V
+	  :INT)))
+    (UNLESS
+	(EQL :SUCCESS
+	     (FOREIGN-ENUM-KEYWORD 'GSL-ERRORNO
+				   STATUS))
+      (WARN 'GSL-WARNING
+	    :GSL-ERRNO
+	    STATUS
+	    :GSL-CONTEXT
+	    `(BESSEL-SEQUENCE-JNU ,NU ,V)))
+    (VALUES
+     (pick-result '((v (double (length v))))))))
+
 
 ;;;;****************************************************************************
 ;;;; Irregular Bessel Function - Fractional Order
diff --git a/special-functions/clausen.lisp b/special-functions/clausen.lisp
new file mode 100644
index 0000000000000000000000000000000000000000..c9a98cc2046e616df1778629d51d7977842e2f93
--- /dev/null
+++ b/special-functions/clausen.lisp
@@ -0,0 +1,16 @@
+;********************************************************
+; file:        clausen.lisp                              
+; description: Clausen function                          
+; date:        Sat Mar 18 2006 - 23:18                   
+; author:      Liam M. Healy                             
+; modified:    Sat Mar 18 2006 - 23:22
+;********************************************************
+;;; $Id:$
+
+(in-package :gsl)
+
+(defun-sf clausen ((x :double))
+  "gsl_sf_clausen_e"
+  :documentation
+  "The Clausen integral @math{Cl_2(x)}."
+  :return (sf-result))
diff --git a/special-functions/coulomb.lisp b/special-functions/coulomb.lisp
new file mode 100644
index 0000000000000000000000000000000000000000..11dc2d4abf641f8c9d0bf013ffb8c8b1846736fa
--- /dev/null
+++ b/special-functions/coulomb.lisp
@@ -0,0 +1,117 @@
+;********************************************************
+; file:        coulumb.lisp                              
+; description: Coulumb functions                         
+; date:        Sat Mar 18 2006 - 23:23                   
+; author:      Liam M. Healy                             
+; modified:    Sun Mar 19 2006 - 00:41
+;********************************************************
+;;; $Id:$
+
+(in-package :gsl)
+
+;;;;****************************************************************************
+;;;; Normalized Hydrogenic Bound States
+;;;;****************************************************************************
+
+(defun-sf hydrogenicR-1 ((x :double) (r :double))
+  "gsl_sf_hydrogenicR_1_e"
+  :documentation
+  "The lowest-order normalized hydrogenic bound state radial wavefunction @c{$R_1 := 2Z \sqrt{Z} \exp(-Z r)$} @math{R_1 := 2Z \sqrt@{Z@} \exp(-Z r)}."
+  :return (sf-result))
+
+(defun-sf hydrogenicR ((n :int) (l :int) (x :double) (r :double))
+  "gsl_sf_hydrogenicR_e"
+  :documentation
+  "The @var{n}-th normalized hydrogenic bound state radial wavefunction,
+$$R_n := {2 Z^{3/2} \over n^2}  \left({2Z \over n}\right)^l  \sqrt{(n-l-1)! \over (n+l)!} \exp(-Z r/n) L^{2l+1}_{n-l-1}(2Z/n r).
+$$
+The normalization is chosen such that the wavefunction @math{\psi} is given by 
+@c{$\psi(n,l,r) = R_n Y_{lm}$}
+@math{\psi(n,l,r) = R_n Y_@{lm@}}."
+  :return (sf-result))
+
+;;;;****************************************************************************
+;;;; Coulomb Wave Functions
+;;;;****************************************************************************
+
+;;; Comments are direct from GSL and aren't lispized yet.
+
+;;; Returns should be reorganized sensibly.
+(defun-sf coulomb-wave-FG ((eta :double) (x :double) (L-F :double) (k :int))
+  "gsl_sf_coulomb_wave_FG_e"
+  :documentation
+  "The Coulomb wave functions @math{F_L(\eta,x)},
+@c{$G_{L-k}(\eta,x)$} 
+@math{G_@{L-k@}(\eta,x)} and their derivatives 
+@math{F'_L(\eta,x)}, 
+@c{$G'_{L-k}(\eta,x)$}
+@math{G'_@{L-k@}(\eta,x)}
+with respect to @math{x}.  The parameters are restricted to @math{L,
+L-k > -1/2}, @math{x > 0} and integer @math{k}.  Note that @math{L}
+itself is not restricted to being an integer. The results are stored in
+the parameters @var{F}, @var{G} for the function values and @var{Fp},
+@var{Gp} for the derivative values.  If an overflow occurs,
+@code{GSL_EOVRFLW} is returned and scaling exponents are stored in
+the modifiable parameters @var{exp_F}, @var{exp_G}."
+  :return
+  (sf-result sf-result sf-result sf-result :double :double))
+
+(defun-sf coulomb-wave-F-array
+    ((L-min :double) (kmax :int) (eta :double) (x :double))
+  "gsl_sf_coulomb_wave_F_array"
+  :documentation
+  "The Coulomb wave function @math{F_L(\eta,x)} for
+@math{L = Lmin \dots Lmin + kmax}, storing the results in @var{fc_array}.
+In the case of overflow the exponent is stored in @var{F_exponent}."
+  :return ((:double (1+ kmax)) :double))
+
+(defun-sf coulomb-wave-FG-array
+    ((L-min :double) (kmax :int) (eta :double) (x :double))
+  "gsl_sf_coulomb_wave_F_array"
+  :documentation
+  "The functions @math{F_L(\eta,x)},
+@math{G_L(\eta,x)} for @math{L = Lmin \dots Lmin + kmax} storing the
+results in @var{fc_array} and @var{gc_array}.  In the case of overflow the
+exponents are stored in @var{F_exponent} and @var{G_exponent}."
+  :return ((:double (1+ kmax)) (:double (1+ kmax)) :double :double))
+
+(defun-sf coulomb-wave-FGp-array 
+    ((L-min :double) (kmax :int) (eta :double) (x :double))
+  "gsl_sf_coulomb_wave_FGp_array"
+  :documentation
+  "The functions @math{F_L(\eta,x)},
+@math{G_L(\eta,x)} and their derivatives @math{F'_L(\eta,x)},
+@math{G'_L(\eta,x)} for @math{L = Lmin \dots Lmin + kmax} storing the
+results in @var{fc_array}, @var{gc_array}, @var{fcp_array} and @var{gcp_array}.
+In the case of overflow the exponents are stored in @var{F_exponent} 
+and @var{G_exponent}."
+  :return ((:double (1+ kmax)) (:double (1+ kmax))
+	   (:double (1+ kmax)) (:double (1+ kmax))
+	   :double :double))
+
+(defun-sf coulomb-wave-sphF-array
+    ((L-min :double) (kmax :int) (eta :double) (x :double))
+  "gsl_sf_coulomb_wave_sphF_array"
+  :documentation
+  "The Coulomb wave function divided by the argument
+@math{F_L(\eta, x)/x} for @math{L = Lmin \dots Lmin + kmax}, storing the
+results in @var{fc_array}.  In the case of overflow the exponent is
+stored in @var{F_exponent}. This function reduces to spherical Bessel
+functions in the limit @math{\eta \to 0}."
+  :return ((:double (1+ kmax)) :double :double))
+
+;;;;****************************************************************************
+;;;; Coulomb Wave Function Normalization Constant
+;;;;****************************************************************************
+
+(defun-sf coulomb-CL ((L :double) (eta :double))
+  "gsl_sf_coulomb_CL_e"
+  :documentation
+  "The Coulomb wave function normalization constant @math{C_L(\eta)} for @math{L > -1}."
+  :return (sf-result))
+
+(defun-sf coulomb-CL-array ((Lmin :double) (kmax :int) (eta :double))
+  "gsl_sf_coulomb_CL_array"
+  :documentation
+  "The Coulomb wave function normalization constant @math{C_L(\eta)} for @math{L = Lmin \dots Lmin + kmax}, @math{Lmin > -1}."
+  :return ((:double (1+ kmax))))