Skip to content
Snippets Groups Projects
Commit cf83a302 authored by liam's avatar liam
Browse files

Interpolation and splines added. New macro defvariable for access to

variables in the library.


git-svn-id: svn+ssh://pop/opt/space/mathematics/gsl/trunk@3243 a3d8a0fb-c1db-0310-ace7-a616afeb9e30
parent 3f6e03fa
No related branches found
No related tags found
No related merge requests found
...@@ -3,7 +3,7 @@ ...@@ -3,7 +3,7 @@
; description: Definition of GSLL system ; description: Definition of GSLL system
; date: ; date:
; author: Liam Healy ; author: Liam Healy
; modified: Thu Oct 18 2007 - 22:18 ; modified: Sat Nov 10 2007 - 21:15
;******************************************************** ;********************************************************
;;; $Id: $ ;;; $Id: $
...@@ -159,4 +159,11 @@ ...@@ -159,4 +159,11 @@
(:file "stepping") (:file "stepping")
(:file "control") (:file "control")
(:file "evolution") (: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")))))
...@@ -85,7 +85,8 @@ to find that the Lisp function name is #'JACOBIAN-ELLIPTIC-FUNCTIONS. ...@@ -85,7 +85,8 @@ to find that the Lisp function name is #'JACOBIAN-ELLIPTIC-FUNCTIONS.
</pre> </pre>
to get an explanation of the arguments etc. to get an explanation of the arguments etc.
</li> </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> </ul>
</div> </div>
...@@ -128,7 +129,7 @@ svn checkout svn://common-lisp.net/project/gsll/subversion/trunk</pre> ...@@ -128,7 +129,7 @@ svn checkout svn://common-lisp.net/project/gsll/subversion/trunk</pre>
<h3>Project Development</h3> <h3>Project Development</h3>
<div class="content"> <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 done, but there is still much work to do. The status is summarized
below, by reference to the Chapters in the GSL manual.</p> below, by reference to the Chapters in the GSL manual.</p>
<p>It is possible that the library API will change but probably not <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> ...@@ -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 also serve as examples. The judgement of floating point numbers as
equivalent is a hack, and needs to be replaced; I am looking for 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 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> tests.</p>
<p>A Trac web page will soon be set up for those interested in <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 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> ...@@ -272,7 +273,7 @@ svn checkout svn://common-lisp.net/project/gsll/subversion/trunk</pre>
</tr> </tr>
<tr> <tr>
<td><a href="Interpolation.html">Interpolation</a></td> <td><a href="Interpolation.html">Interpolation</a></td>
<td>Not started</td> <td>Done</td>
</tr> </tr>
<tr> <tr>
<td><a href="Numerical-Differentiation.html">Numerical Differentiation</a></td> <td><a href="Numerical-Differentiation.html">Numerical Differentiation</a></td>
......
...@@ -3,7 +3,7 @@ ...@@ -3,7 +3,7 @@
; description: Macros to interface GSL functions. ; description: Macros to interface GSL functions.
; date: Mon Mar 6 2006 - 22:35 ; date: Mon Mar 6 2006 - 22:35
; author: Liam M. Healy ; author: Liam M. Healy
; modified: Sat Sep 29 2007 - 22:43 ; modified: Sun Nov 4 2007 - 17:50
;******************************************************** ;********************************************************
(in-package :gsl) (in-package :gsl)
...@@ -215,10 +215,6 @@ ...@@ -215,10 +215,6 @@
(map-name ',(or index name) ,gsl-name) (map-name ',(or index name) ,gsl-name)
,@(when export `((export ',name)))))) ,@(when export `((export ',name))))))
;;;;****************************************************************************
;;;; Macro defun-gsl
;;;;****************************************************************************
(defmacro defun-optionals (defmacro defun-optionals
(name arglist no-optional optionals &optional documentation) (name arglist no-optional optionals &optional documentation)
"Define a function with and without optional arguments." "Define a function with and without optional arguments."
...@@ -234,3 +230,14 @@ ...@@ -234,3 +230,14 @@
(,(intern (,(intern
(concatenate 'string (string name) (string no-optional))) (concatenate 'string (string name) (string no-optional)))
,@mandatory-arglist))))) ,@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)))
;********************************************************
; 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.")
;********************************************************
; 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}.")
;********************************************************
; 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.")
;********************************************************
; 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)))))))
;********************************************************
; 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.")
...@@ -3,7 +3,7 @@ ...@@ -3,7 +3,7 @@
; description: Stepping functions for ODE systems ; description: Stepping functions for ODE systems
; date: Mon Sep 24 2007 - 21:33 ; date: Mon Sep 24 2007 - 21:33
; author: Liam Healy ; author: Liam Healy
; modified: Tue Oct 16 2007 - 22:45 ; modified: Sun Nov 4 2007 - 17:54
;******************************************************** ;********************************************************
;;; $Id: $ ;;; $Id: $
...@@ -78,39 +78,34 @@ ...@@ -78,39 +78,34 @@
those from @code{gsl_odeiv_step_apply} itself, any user-defined return those from @code{gsl_odeiv_step_apply} itself, any user-defined return
values should be distinct from the standard GSL error codes.") values should be distinct from the standard GSL error codes.")
(defmacro defstep (cl-symbol gsl-symbol documentation) (defvariable *step-rk2* "gsl_odeiv_step_rk2"
`(progn
(cffi:defcvar (,gsl-symbol ,cl-symbol) :pointer :read-only t)
(setf (documentation ',cl-symbol 'variable) ,documentation)))
(defstep *step-rk2* "gsl_odeiv_step_rk2"
"Embedded Runge-Kutta (2, 3) method.") "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.") "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 "Embedded Runge-Kutta-Fehlberg (4, 5) method. This method is a good
general-purpose integrator.") 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.") "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.") "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.") "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.") "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 "Implicit Bulirsch-Stoer method of Bader and Deuflhard. This algorithm
requires the Jacobian.") requires the Jacobian.")
(defstep *step-gear1* "gsl_odeiv_step_gear1" (defvariable *step-gear1* "gsl_odeiv_step_gear1"
"M=1 implicit Gear method.") "M=1 implicit Gear method.")
(defstep *step-gear2* "gsl_odeiv_step_gear2" (defvariable *step-gear2* "gsl_odeiv_step_gear2"
"M=2 implicit Gear method.") "M=2 implicit Gear method.")
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment