diff --git a/gsll.asd b/gsll.asd
index 313e172b685751e392a760827f9b5b5160f80263..76131314ff113e48910ee105e21f89add6b1b667 100644
--- a/gsll.asd
+++ b/gsll.asd
@@ -3,7 +3,7 @@
 ; description: Definition of GSLL system 
 ; date:        
 ; author:      Liam Healy
-; modified:    Thu Oct 18 2007 - 22:18
+; modified:    Sat Nov 10 2007 - 21:15
 ;********************************************************
 ;;; $Id: $
 
@@ -159,4 +159,11 @@
 	     (:file "stepping")
 	     (:file "control")
 	     (:file "evolution")
-	     (:file "ode-example" :depends-on (ode-system stepping))))))
+	     (:file "ode-example" :depends-on (ode-system stepping))))
+   (:module interpolation
+	    :depends-on (init)
+	    :components
+	    ((:file "interpolation")
+	     (:file "types")
+	     (:file "lookup")
+	     (:file "evaluation")))))
diff --git a/index.html b/index.html
index f38da0519cab83a557d580db413e5abb9edf6fa0..148b7bcc28c8b5d17b87be3d17a4ec0cadf49b30 100644
--- a/index.html
+++ b/index.html
@@ -85,7 +85,8 @@ to find that the Lisp function name is #'JACOBIAN-ELLIPTIC-FUNCTIONS.
    </pre>
    to get an explanation of the arguments etc.
  </li>
- <li>Look at the end of the Lisp source file for examples.</li>
+ <li>Look at the end of the Lisp source file, or in an "example" file
+ in the same directory, for examples.</li>
  </ul>
  </div>
 
@@ -128,7 +129,7 @@ svn checkout svn://common-lisp.net/project/gsll/subversion/trunk</pre>
  
  <h3>Project Development</h3>
   <div class="content">
-  <p>This project has been in development for a year.  A majority is
+  <p>This project has been in development since March 2006.  A majority is
   done, but there is still much work to do.  The status is summarized
   below, by reference to the Chapters in the GSL manual.</p>
    <p>It is possible that the library API will change but probably not
@@ -148,7 +149,7 @@ svn checkout svn://common-lisp.net/project/gsll/subversion/trunk</pre>
    also serve as examples.  The judgement of floating point numbers as
    equivalent is a hack, and needs to be replaced; I am looking for
    good ideas.  Note that most of the result have not been verified
-   with an indpendent source; the values given are from previous
+   with an independent source; the values given are from previous
    tests.</p>
   <p>A Trac web page will soon be set up for those interested in
   development, and for reporting bugs.  In the meantime, send email to
@@ -272,7 +273,7 @@ svn checkout svn://common-lisp.net/project/gsll/subversion/trunk</pre>
 	</tr>
 	<tr>
 	<td><a href="Interpolation.html">Interpolation</a></td>
-	<td>Not started</td>
+	<td>Done</td>
 	</tr>
 	<tr>
 	<td><a href="Numerical-Differentiation.html">Numerical Differentiation</a></td>
diff --git a/init/interface.lisp b/init/interface.lisp
index 13a088df4b5c594d0c050154ad9a1e26f03aa66d..61654aa44c1689f6c634c10e4fb42e756b55e683 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:    Sat Sep 29 2007 - 22:43
+; modified:    Sun Nov  4 2007 - 17:50
 ;********************************************************
 
 (in-package :gsl)
@@ -215,10 +215,6 @@
        (map-name ',(or index name) ,gsl-name)
        ,@(when export `((export ',name))))))
 
-;;;;****************************************************************************
-;;;; Macro defun-gsl 
-;;;;****************************************************************************
-
 (defmacro defun-optionals
     (name arglist no-optional optionals &optional documentation)
   "Define a function with and without optional arguments."
@@ -234,3 +230,14 @@
 	   (,(intern
 	      (concatenate 'string (string name) (string no-optional)))
 	     ,@mandatory-arglist)))))
+
+;;;;****************************************************************************
+;;;; Variables in library
+;;;;****************************************************************************
+
+(defmacro defvariable (cl-symbol gsl-symbol documentation)
+  "Define a library variable pointer."
+  `(progn
+    (cffi:defcvar (,gsl-symbol ,cl-symbol) :pointer :read-only t)
+    (map-name ',cl-symbol ,gsl-symbol)
+    (setf (documentation ',cl-symbol 'variable) ,documentation)))
diff --git a/interpolation/evaluation.lisp b/interpolation/evaluation.lisp
new file mode 100644
index 0000000000000000000000000000000000000000..fd6864ffb2fe70ce3bd548b2630abe4935f6baa6
--- /dev/null
+++ b/interpolation/evaluation.lisp
@@ -0,0 +1,88 @@
+;********************************************************
+; file:        evaluation.lisp                           
+; description: Evaluation of interpolation functions.    
+; date:        Sun Nov  4 2007 - 18:40                   
+; author:      Liam Healy                                
+; modified:    Sun Nov  4 2007 - 19:04
+;********************************************************
+;;; $Id: $
+
+(in-package :gsl)
+
+(defun-gsl evaluate-interpolation (interpolation xa ya x acceleration)
+  "gsl_interp_eval"
+  ((interpolation :pointer) (xa :pointer) (ya :pointer) (x :double)
+   (acceleration :pointer))
+  :c-return :double
+  :documentation
+  "Find the interpolated value of @var{y} for a given
+   point @var{x}, using the interpolation object 'interpolation, data arrays
+   @var{xa} and @var{ya} and the accelerator 'acceleration.")
+
+(defun-gsl derivative-evaluate-interpolation (interpolation xa ya x acceleration)
+  "gsl_interp_eval_deriv"
+  ((interpolation :pointer) (xa :pointer) (ya :pointer) (x :double)
+   (acceleration :pointer))
+  :c-return :double
+  :documentation
+  "Find the derivative of an interpolated function for a given point
+   @var{x}, using the interpolation object 'interpolation, data arrays
+   @var{xa} and @var{ya} and the accelerator 'acceleration.")
+
+(defun-gsl second-derivative-evaluate-interpolation (interpolation xa ya x acceleration)
+  "gsl_interp_eval_deriv2"
+  ((interpolation :pointer) (xa :pointer) (ya :pointer) (x :double)
+   (acceleration :pointer))
+  :c-return :double
+  :documentation
+  "Find the second derivative of an interpolated function for a given point
+   @var{x}, using the interpolation object 'interpolation, data arrays
+   @var{xa} and @var{ya} and the accelerator 'acceleration.")
+
+(defun-gsl integral-evaluate-interpolation (interpolation xa ya x low high acceleration)
+  "gsl_interp_eval_integ"
+  ((interpolation :pointer) (xa :pointer) (ya :pointer) (x :double)
+   (low :double) (high :double) (acceleration :pointer))
+  :c-return :double
+  :documentation
+  "Find the numerical integral of an interpolated function over the
+   range [low, high], using the interpolation object 'interpolation,
+   data arrays @var{xa} and @var{ya} and the accelerator 'acceleration.")
+
+;;; Spline
+(defun-gsl evaluate-spline (spline x acceleration)
+  "gsl_spline_eval"
+  ((spline :pointer) (x :double) (acceleration :pointer))
+  :c-return :double
+  :documentation
+  "Find the interpolated value of @var{y} for a given
+   point @var{x}, using the spline object 'spline, data arrays
+   @var{xa} and @var{ya} and the accelerator 'acceleration.")
+
+(defun-gsl derivative-evaluate-spline (spline x acceleration)
+  "gsl_spline_eval_deriv"
+  ((spline :pointer) (x :double) (acceleration :pointer))
+  :c-return :double
+  :documentation
+  "Find the derivative of an interpolated function for a given point
+   @var{x}, using the spline object 'spline, data arrays
+   @var{xa} and @var{ya} and the accelerator 'acceleration.")
+
+(defun-gsl second-derivative-evaluate-spline (spline x acceleration)
+  "gsl_spline_eval_deriv2"
+  ((spline :pointer) (x :double) (acceleration :pointer))
+  :c-return :double
+  :documentation
+  "Find the second derivative of an interpolated function for a given point
+   @var{x}, using the spline object 'spline, data arrays
+   @var{xa} and @var{ya} and the accelerator 'acceleration.")
+
+(defun-gsl integral-evaluate-spline (spline x low high acceleration)
+  "gsl_spline_eval_integ"
+  ((spline :pointer) (x :double) (low :double) (high :double)
+   (acceleration :pointer))
+  :c-return :double
+  :documentation
+  "Find the numerical integral of an interpolated function over the
+   range [low, high], using the spline object 'spline,
+   data arrays @var{xa} and @var{ya} and the accelerator 'acceleration.")
diff --git a/interpolation/interpolation.lisp b/interpolation/interpolation.lisp
new file mode 100644
index 0000000000000000000000000000000000000000..9147c0d94565f44ecb4614b38ab5a7ee0b6b55b6
--- /dev/null
+++ b/interpolation/interpolation.lisp
@@ -0,0 +1,69 @@
+;********************************************************
+; file:        interpolation.lisp                        
+; description: Interpolation allocation, initialization, 
+;              and freeing.                              
+; date:        Sun Nov  4 2007 - 17:24                   
+; author:      Liam Healy                                
+; modified:    Sun Nov  4 2007 - 19:00
+;********************************************************
+;;; $Id: $
+
+(in-package :gsl)
+
+;;; A spline is an interpolation that also stores the arrays xa and ya,
+;;; so they need not be supplied on each call.
+
+
+;;; Interpolation
+(defun-gsl allocate-interpolation (type size)
+  "gsl_interp_alloc"
+  ((type :pointer) (size :size))
+  :c-return :pointer
+  :documentation
+  "Allocate an interpolation object of type for size data-points,
+   and return the pointer to it.")
+
+(defun-gsl initialize-interpolation (interpolation xa ya size)
+  "gsl_interp_init"
+  ((interpolation :pointer) (xa :pointer) (ya :pointer) (size :size))
+  :documentation
+  "Initialize the interpolation object @var{interp} for the
+  data (@var{xa},@var{ya}) where @var{xa} and @var{ya} are arrays of size
+  @var{size}.  The interpolation object (@code{gsl_interp}) does not save
+  the data arrays @var{xa} and @var{ya} and only stores the static state
+  computed from the data.  The @var{xa} data array is always assumed to be
+  strictly ordered; the behavior for other arrangements is not defined.")
+
+(defun-gsl free-interpolation (interpolation)
+  "gsl_interp_free"
+  ((interpolation :pointer))
+  :c-return :void
+  :documentation
+  "Frees the interpolation object @var{interp}.")
+
+;;; Spline
+(defun-gsl allocate-spline (type size)
+  "gsl_spline_alloc"
+  ((type :pointer) (size :size))
+  :c-return :pointer
+  :documentation
+  "Allocate an interpolation object of type for size data-points,
+   and return the pointer to it.")
+
+(defun-gsl initialize-spline (interpolation xa ya size)
+  "gsl_spline_init"
+  ((interpolation :pointer) (xa :pointer) (ya :pointer) (size :size))
+  :documentation
+  "Initialize the interpolation object @var{interp} for the
+  data (@var{xa},@var{ya}) where @var{xa} and @var{ya} are arrays of size
+  @var{size}.  The interpolation object (@code{gsl_spline}) does not save
+  the data arrays @var{xa} and @var{ya} and only stores the static state
+  computed from the data.  The @var{xa} data array is always assumed to be
+  strictly ordered; the behavior for other arrangements is not defined.")
+
+(defun-gsl free-spline (interpolation)
+  "gsl_spline_free"
+  ((interpolation :pointer))
+  :c-return :void
+  :documentation
+  "Frees the interpolation object @var{interp}.")
diff --git a/interpolation/lookup.lisp b/interpolation/lookup.lisp
new file mode 100644
index 0000000000000000000000000000000000000000..a3aadd3d1d6fc51629044ca58d7b2120e297e895
--- /dev/null
+++ b/interpolation/lookup.lisp
@@ -0,0 +1,49 @@
+;********************************************************
+; file:        lookup.lisp                               
+; description: Index lookup and acceleration             
+; date:        Sun Nov  4 2007 - 18:09                   
+; author:      Liam Healy                                
+; modified:    Sat Nov 10 2007 - 21:16
+;********************************************************
+;;; $Id: $
+
+(in-package :gsl)
+
+;;; Lacks with-acceleration macro.
+
+(defun-gsl interpolation-search (x-array x low-index high-index)
+  "gsl_interp_bsearch"
+  ((x-array :pointer) (x :double) (low-index :size) (high-index :size))
+  :c-return :size
+  :documentation
+  "Find the index @math{i} of the array @var{x-array} such
+   that @code{x-array[i] <= x < x-array[i+1]}.  The index is searched for
+   in the range [low-index, high-index].")
+
+(defun-gsl allocate-acceleration ()
+  "gsl_interp_accel_alloc"
+  ()
+  :c-return :pointer
+  :documentation
+  "Allocate an accelerator object, which is a
+   kind of iterator for interpolation lookups.  It tracks the state of
+   lookups, thus allowing for application of various acceleration
+   strategies.")
+
+(defun-gsl accelerated-interpolation-search (x-array x acceleration)
+  "gsl_interp_accel_find"
+  ((acceleration :pointer) (x-array :pointer) (x :double))
+  :c-return :size
+  :documentation
+  "Search the data array @var{x_array}
+   of size @var{size}, using the given accelerator @var{a}.  This is how
+   lookups are performed during evaluation of an interpolation.  The
+   function returns an index @math{i} such that @code{x_array[i] <= x <
+   x_array[i+1]}.")
+
+(defun-gsl free-acceleration (acceleration)
+  "gsl_interp_accel_free"
+  ((acceleration :pointer))
+  :c-return :void
+  :documentation
+  "Frees the accelerator object.")
diff --git a/interpolation/spline-example.lisp b/interpolation/spline-example.lisp
new file mode 100644
index 0000000000000000000000000000000000000000..2ba7a9dd41db4f9ff0898ae50f840bd39ff18e6e
--- /dev/null
+++ b/interpolation/spline-example.lisp
@@ -0,0 +1,38 @@
+;********************************************************
+; file:        spline-example.lisp                       
+; description: Example spline                            
+; date:        Sat Nov 10 2007 - 21:18                   
+; author:      Liam Healy                                
+; modified:    Sat Nov 10 2007 - 22:10
+;********************************************************
+;;; $Id: $
+
+(in-package :gsl)
+
+;;; Possible future improvements:
+;;; Direct double-float vectors.
+;;; Use generic alloc/free defined in data.lisp
+
+(defun spline-example-arrays ()
+  (let ((xarr (make-array '(10) :element-type 'double-float))
+	(yarr (make-array '(10) :element-type 'double-float)))
+    (loop for i from 0 below 10
+	  do (setf (aref xarr i) (+ i (* 0.5d0 (sin (coerce i 'double-float))))
+		   (aref yarr i) (+ i (cos (expt (coerce i 'double-float) 2)))))
+    (values xarr yarr)))
+
+(defun spline-example ()
+  "The first example in Sec. 26.7 of the GSL manual."
+  (let ((acc (allocate-acceleration))
+	(spline (allocate-spline *cubic-spline-interpolation* 10)))
+    (multiple-value-bind (xarr yarr)
+	(spline-example-arrays)
+      (with-data (cxarr vector-double 10)
+	(with-data (cyarr vector-double 10)
+	  (setf (data cxarr) xarr (data cyarr) yarr)
+	  (initialize-spline spline (gsl-array cxarr) (gsl-array cyarr) 10)
+	  (prog1
+	      (loop for xi from (aref xarr 0) below (aref xarr 9) by 0.01d0
+		    collect (list xi (evaluate-spline spline xi acc)))
+	    (free-spline spline)
+	    (free-acceleration acc)))))))
diff --git a/interpolation/types.lisp b/interpolation/types.lisp
new file mode 100644
index 0000000000000000000000000000000000000000..38231554796531518ee160635f29b18297bd1ec5
--- /dev/null
+++ b/interpolation/types.lisp
@@ -0,0 +1,76 @@
+;********************************************************
+; file:        types.lisp                                
+; description: Interpolation types                       
+; date:        Sun Nov  4 2007 - 17:41                   
+; author:      Liam Healy                                
+; modified:    Sun Nov  4 2007 - 19:02
+;********************************************************
+;;; $Id: $
+
+(in-package :gsl)
+
+(defvariable *linear-interpolation* "gsl_interp_linear"
+  "Linear interpolation.  This interpolation method does not require any
+   additional memory.")
+
+(defvariable *polynomial-interpolation* "gsl_interp_polynomial"
+  "Polynomial interpolation.  This method should only be used for
+   interpolating small numbers of points because polynomial interpolation
+   introduces large oscillations, even for well-behaved datasets.  The
+   number of terms in the interpolating polynomial is equal to the number
+   of points.")
+
+(defvariable *cubic-spline-interpolation* "gsl_interp_cspline"
+  "Cubic spline with natural boundary conditions.  The resulting curve is
+   piecewise cubic on each interval, with matching first and second
+   derivatives at the supplied data-points.  The second derivative is
+   chosen to be zero at the first point and last point.")
+
+(defvariable *periodic-cubic-spline-interpolation* "gsl_interp_cspline_periodic"
+  "Cubic spline with periodic boundary conditions.  The resulting curve
+   is piecewise cubic on each interval, with matching first and second
+   derivatives at the supplied data-points.  The derivatives at the first
+   and last points are also matched.  Note that the last point in the
+   data must have the same y-value as the first point, otherwise the
+   resulting periodic interpolation will have a discontinuity at the
+   boundary.")
+
+(defvariable *akima-interpolation* "gsl_interp_akima"
+  "Non-rounded Akima spline with natural boundary conditions.  This method
+   uses the non-rounded corner algorithm of Wodicka.")
+
+(defvariable *periodic-akima-interpolation* "gsl_interp_akima_periodic"
+  "Non-rounded Akima spline with periodic boundary conditions.  This method
+   uses the non-rounded corner algorithm of Wodicka.")
+
+(defun-gsl interpolation-name (interpolation)
+  "gsl_interp_name"
+  ((interpolation :pointer))
+  :c-return :string
+  :documentation
+  "This function returns the name of the interpolation type.")
+
+(defun-gsl interpolation-minimum-size (interpolation)
+  "gsl_interp_min_size"
+  ((interpolation :pointer))
+  :c-return :uint
+  :documentation
+  "The minimum number of points required by the
+   interpolation.  For example, Akima spline interpolation
+   requires a minimum of 5 points.")
+
+(defun-gsl spline-name (interpolation)
+  "gsl_spline_name"
+  ((interpolation :pointer))
+  :c-return :string
+  :documentation
+  "This function returns the name of the interpolation type.")
+
+(defun-gsl spline-minimum-size (interpolation)
+  "gsl_spline_min_size"
+  ((interpolation :pointer))
+  :c-return :uint
+  :documentation
+  "The minimum number of points required by the
+   interpolation.  For example, Akima spline interpolation
+   requires a minimum of 5 points.")
diff --git a/ordinary-differential-equations/stepping.lisp b/ordinary-differential-equations/stepping.lisp
index 8aa01d92f6211f0bec1493be58356cf14ece8a55..60082a1bbd11f0a49c5b1302508b3f64e90787de 100644
--- a/ordinary-differential-equations/stepping.lisp
+++ b/ordinary-differential-equations/stepping.lisp
@@ -3,7 +3,7 @@
 ; description: Stepping functions for ODE systems        
 ; date:        Mon Sep 24 2007 - 21:33                   
 ; author:      Liam Healy                                
-; modified:    Tue Oct 16 2007 - 22:45
+; modified:    Sun Nov  4 2007 - 17:54
 ;********************************************************
 ;;; $Id: $
 
@@ -78,39 +78,34 @@
    those from @code{gsl_odeiv_step_apply} itself, any user-defined return
    values should be distinct from the standard GSL error codes.")
 
-(defmacro defstep (cl-symbol gsl-symbol documentation)
-  `(progn
-    (cffi:defcvar (,gsl-symbol ,cl-symbol) :pointer :read-only t)
-    (setf (documentation ',cl-symbol 'variable) ,documentation)))
-
-(defstep *step-rk2* "gsl_odeiv_step_rk2"
+(defvariable *step-rk2* "gsl_odeiv_step_rk2"
   "Embedded Runge-Kutta (2, 3) method.")
 
-(defstep *step-rk4* "gsl_odeiv_step_rk4"
+(defvariable *step-rk4* "gsl_odeiv_step_rk4"
   "4th order (classical) Runge-Kutta.")
 
-(defstep *step-rkf45* "gsl_odeiv_step_rkf45"
+(defvariable *step-rkf45* "gsl_odeiv_step_rkf45"
   "Embedded Runge-Kutta-Fehlberg (4, 5) method.  This method is a good
    general-purpose integrator.")
 
-(defstep *step-rkck* "gsl_odeiv_step_rkck"
+(defvariable *step-rkck* "gsl_odeiv_step_rkck"
   "Embedded Runge-Kutta Cash-Karp (4, 5) method.")
 
-(defstep *step-rk8pd* "gsl_odeiv_step_rk8pd"
+(defvariable *step-rk8pd* "gsl_odeiv_step_rk8pd"
   "Embedded Runge-Kutta Prince-Dormand (8,9) method.")
 
-(defstep *step-rk2imp* "gsl_odeiv_step_rk2imp"
+(defvariable *step-rk2imp* "gsl_odeiv_step_rk2imp"
   "Implicit 2nd order Runge-Kutta at Gaussian points.")
 
-(defstep *step-rk4imp* "gsl_odeiv_step_rk4imp"
+(defvariable *step-rk4imp* "gsl_odeiv_step_rk4imp"
   "Implicit 4th order Runge-Kutta at Gaussian points.")
 
-(defstep *step-bsimp* "gsl_odeiv_step_bsimp"
+(defvariable *step-bsimp* "gsl_odeiv_step_bsimp"
   "Implicit Bulirsch-Stoer method of Bader and Deuflhard.  This algorithm
    requires the Jacobian.")
 
-(defstep *step-gear1* "gsl_odeiv_step_gear1"
+(defvariable *step-gear1* "gsl_odeiv_step_gear1"
   "M=1 implicit Gear method.")
 
-(defstep *step-gear2* "gsl_odeiv_step_gear2"
+(defvariable *step-gear2* "gsl_odeiv_step_gear2"
   "M=2 implicit Gear method.")